1、東南大學(xué)電氣工程學(xué)院MATLAB數(shù)學(xué)建模 實(shí)驗(yàn)報(bào)告三種方法解析Volterra方程劉海東（16006213）2008/12/13一、 實(shí)驗(yàn)?zāi)康模和ㄟ^MATLAB實(shí)現(xiàn)書上含初值一階Volterra方程組的解析，在解析過程之中注意加深對三種方法：向前歐拉公式、改進(jìn)的歐拉公式和四階龍格-庫塔公式的理解，注意各個(gè)方法的細(xì)節(jié)和精度比較。進(jìn)一步熟悉MATLAB編程。二、 實(shí)驗(yàn)原理：在數(shù)學(xué)建模理論課上，我們在書上6.4節(jié)中介紹了Volterra模型，并用相平面分析法對問題作了簡要分析。隨后我們又在6.7節(jié)中學(xué)習(xí)到了常微分方程數(shù)值解，學(xué)習(xí)了歐拉方法和龍格-庫塔方法，所以我們嘗試借助6.7節(jié)中的數(shù)值方法來考慮如

2、下Volterra方程：為了跟書上的解析圖相對應(yīng)，我選取區(qū)間；步長=0.1。1) 向前歐拉公式(1階精度)：我們先根據(jù)步長將區(qū)間分為150等分，在其中任一區(qū)間上（）取對應(yīng)左端點(diǎn)的和作為遞推和的初值，將其帶入到向前歐拉公式組得到如下方程組：這樣，我們就可以根據(jù)初值層層遞推出對應(yīng)后面150個(gè)t值的x值。具體在MATLAB中的實(shí)現(xiàn)過程：首先將各自初值賦予初始變量，并將步長輸入。根據(jù)步長，在取定循環(huán)值t=0.1:0.1:15形成一個(gè)150次的循環(huán)，循環(huán)過程為：根據(jù)上面方程組的前兩行帶入初值求解下一次數(shù)值，然后將下次數(shù)值作為初值不停迭代下去，即得到150個(gè)要求的值。因?yàn)榭倲?shù)據(jù)太過龐大，我數(shù)據(jù)輸出只輸出了

3、t在0,2上的數(shù)據(jù)，為使結(jié)果更為形象直觀，我在把全部結(jié)果其標(biāo)注在圖形上，并隨后進(jìn)行30階擬合成Volterra解析曲線。2) 改進(jìn)歐拉公式(2階精度)：為在提高精度的同時(shí)不使計(jì)算復(fù)雜，改進(jìn)的歐拉公式將向前歐拉公式的簡單與梯形公式的高精度結(jié)合在一起。它把迭代過程簡化成兩步：先由向前歐拉公式算出x(k)的預(yù)測值，再把預(yù)測值代入梯形公式的右端進(jìn)行校正，得到如下方程組：分組遞推迭代過程跟向前歐拉公式一樣。具體在MATLAB中實(shí)現(xiàn)過程也跟向前歐拉公式類似，只不過隨著方程的改變而每個(gè)循環(huán)中間多做了計(jì)算而已。3) 4階龍格-庫塔公式(4階精度)：由上面兩個(gè)公式得到啟發(fā)，取四個(gè)點(diǎn)的斜率值作加權(quán)平均來作為斜率的

4、近似值，這樣我們就提高了精度。所以我們在每一小區(qū)間上取四個(gè)點(diǎn)構(gòu)造得到下面4階龍格-庫塔公式：分組遞推迭代過程同上。而在MATLAB中實(shí)現(xiàn)過程依然和上面類似，只不過相比改進(jìn)歐拉公式每次循環(huán)中計(jì)算過程更為復(fù)雜而已。三、 實(shí)驗(yàn)源代碼：1）向前歐拉公式函數(shù)文件Volterra1:function Volterra1%定義函數(shù)；display(步長為：);%顯示步長；h=0.1,k=1;t=0;xx1(1)=25;x1(1)=25;xx2(1)=2;x2(1)=2;%賦予初值；for t=0.1:h:15;%進(jìn)入150次循環(huán)； x1(k+1)=x1(k)+h*x1(k)*(1-0.1*x2(k); x2

5、(k+1)=x2(k)+h*x2(k)*(0.02*x1(k)-0.5);%代入向前歐拉公式； k=k+1; if t=2 xx1(k)=x1(k); xx2(k)=x2(k); end %方便前21對數(shù)值的輸出endt=0:0.1:15;xx1xx2%顯示t在0,2上的x值；plot(t,x1,x,t,x2,x);%描出全部的點(diǎn)值；hold onp=polyfit(t,x1,30);tt=0:0.01:15;xx1=polyval(p,tt);p1=polyfit(t,x2,30);xx2=polyval(p1,tt);plot(tt,xx1,-r,tt,xx2,-r);%30階擬合曲線xl

6、abel(t);ylabel(x1,x2);title(向前歐拉公式擬合Volterra曲線（30階）);clear%繪圖，清理變量；2）改進(jìn)歐拉公式函數(shù)文件Volterra2:function Volterra2%定義函數(shù)；display(步長為：);%顯示步長；h=0.1,k=1;t=0;xx1(1)=25;x1(1)=25;xx2(1)=2;x2(1)=2;%賦予初值；for t=0.1:h:15;%進(jìn)入150次循環(huán)； x11=x1(k)+h*x1(k)*(1-0.1*x2(k); x21=x2(k)+h*x2(k)*(0.02*x1(k)-0.5); x1(k+1)=x1(k)+h*(

7、x1(k)*(1-0.1*x2(k)+x11*(1-0.1*x21)/2;x2(k+1)=x2(k)+h*(x2(k)*(0.02*x1(k)-0.5)+x21*(0.02*x11-0.5)/2;%代入改進(jìn)歐拉公式； k=k+1; if t=2 xx1(k)=x1(k); xx2(k)=x2(k); end %方便21對數(shù)值的輸出；endt=0:0.1:15;xx1xx2%顯示t在0,2上的x值；plot(t,x1,x,t,x2,x);%描出全部的點(diǎn)值；hold onp=polyfit(t,x1,30);tt=0:0.01:15;xx1=polyval(p,tt);p1=polyfit(t,x

8、2,30);xx2=polyval(p1,tt);plot(tt,xx1,-k,tt,xx2,-k,LineWidth,2); %30階擬合曲線；xlabel(t);ylabel(x1,x2);title(改進(jìn)歐拉公式擬合Volterra曲線（30階）);clear%標(biāo)注圖形，清理變量。3）龍格-庫塔公式函數(shù)文件Volterra4:function Volterra4%定義函數(shù)；display(步長為：);%顯示步長；h=0.1,k=1;t=0;xx1(1)=25;x1(1)=25;xx2(1)=2;x2(1)=2;%賦予初值；for t=0.1:h:15;%進(jìn)入150次循環(huán)； k11=x1(

9、k)*(1-0.1*x2(k); k21=x2(k)*(0.02*x1(k)-0.5); k12=(x1(k)+h/2*k11)*(1-0.1*(x2(k)+h/2*k21); k22=(x2(k)+h/2*k21)*(0.02*(x1(k)+h/2*k11)-0.5); k13=(x1(k)+h/2*k12)*(1-0.1*(x2(k)+h/2*k22); k23=(x2(k)+h/2*k22)*(0.02*(x1(k)+h/2*k12)-0.5); k14=(x1(k)+h/2*k13)*(1-0.1*(x2(k)+h/2*k23); k24=(x2(k)+h/2*k23)*(0.02*(

10、x1(k)+h/2*k13)-0.5); x1(k+1)=x1(k)+h/6*(k11+2*k12+2*k13+k14); x2(k+1)=x2(k)+h/6*(k21+2*k22+2*k23+k24);%代入龍格-庫塔公式； k=k+1; if t=2 xx1(k)=x1(k); xx2(k)=x2(k); end %方便前21對數(shù)值的輸出；endt=0:0.1:15;xx1xx2%顯示t在0,2上的x值；plot(t,x1,x,t,x2,x);%描出全部點(diǎn)值；hold onp=polyfit(t,x1,30);tt=0:0.01:15;xx1=polyval(p,tt);p1=polyfi

11、t(t,x2,30);xx2=polyval(p1,tt);plot(tt,xx1,-g,tt,xx2,-g,LineWidth,2); %30階擬合曲線；xlabel(t);ylabel(x1,x2);title(改進(jìn)歐拉公式擬合Volterra曲線（30階）);clear%標(biāo)注圖形，清理變量。四、 實(shí)驗(yàn)結(jié)果：因?yàn)楸緦?shí)驗(yàn)的數(shù)據(jù)實(shí)在太多，我只取了前21對數(shù)值進(jìn)行顯示，而為驗(yàn)證實(shí)驗(yàn)的正確性，我描出了所有的點(diǎn)值，并采用30階擬合逼近得出各種方法下所算出的解析曲線從而可以與書上的理論曲線比對來驗(yàn)證。1）向前歐拉公式Volterra1運(yùn)行結(jié)果:步長為：h =0.1000xx1 = Columns 1

12、through 11 25.0000 27.0000 29.1600 31.4905 34.0019 36.7047 39.6089 42.7234 46.0561 49.6130 53.3969 Columns 12 through 21 57.4071 61.6372 66.0740 70.6944 75.4633 80.3295 85.2216 90.0428 94.6658 98.9270xx2 = Columns 1 through 11 2.0000 2.0000 2.0080 2.0247 2.0510 2.0879 2.1368 2.1992 2.2772 2.3731 2.4

13、899 Columns 12 through 21 2.6313 2.8018 3.0072 3.2542 3.5516 3.9100 4.3427 4.8658 5.4987 6.26492）改進(jìn)歐拉公式Volterra2運(yùn)行結(jié)果:步長為：h =0.1000xx1 = Columns 1 through 11 25.0000 27.0800 29.3307 31.7629 34.3875 37.2146 40.2534 43.5118 46.9950 50.7051 54.6396 Columns 12 through 21 58.7894 63.1371 67.6538 72.2959 7

14、7.0011 81.6833 86.2278 90.4858 94.2714 97.3607xx2 = Columns 1 through 11 2.0000 2.0040 2.0167 2.0390 2.0720 2.1170 2.1757 2.2502 2.3429 2.4570 2.5965 Columns 12 through 21 2.7661 2.9719 3.2214 3.5240 3.8912 4.3371 4.8790 5.5372 6.3352 7.2986 3）龍格-庫塔公式Volterra4運(yùn)行結(jié)果:步長為：h =0.1000xx1 = Columns 1 throug

15、h 11 25.0000 27.0680 29.3051 31.7221 34.3297 37.1381 40.1568 43.3935 46.8543 50.5416 54.4539 Columns 12 through 21 58.5834 62.9145 67.4208 72.0622 76.7806 81.4956 86.0991 90.4507 94.3728 97.6500xx2 = Columns 1 through 11 2.0000 2.0034 2.0154 2.0369 2.0689 2.1127 2.1698 2.2424 2.3328 2.4440 2.5798 Columns 12 through 212.7448 2.9449 3.1871 3.4805 3.8361 4.2675 4.7912 5.4270 6.1978 7.1290五、 實(shí)驗(yàn)結(jié)果分析與實(shí)驗(yàn)小結(jié)：1）結(jié)果分析：為了驗(yàn)證程序的正確性，我們將上面的三個(gè)圖分別與書上圖6.16比對發(fā)現(xiàn)