三種方法解決voterra方程參考模板

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

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

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

4、來作為斜率的近似值，這樣我們就提高了精度。所以我們在每一小區(qū)間上取四個點構(gòu)造得到下面4階龍格-庫塔公式：分組遞推迭代過程同上。而在MATLAB中實現(xiàn)過程依然和上面類似，只不過相比改進歐拉公式每次循環(huá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;%進入150次循環(huán)； x1(k+1)=x1(k)+h*x1(k)*

5、(1-0.1*x2(k); x2(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');%描出全部的點值；hold onp=polyfit(t,x1,30);tt=0:0.01:15;xx1=polyval(p,tt);p1=polyfit(t,x2,30);xx2=polyval(p1,

6、tt);plot(tt,xx1,'-r',tt,xx2,'-r');%30階擬合曲線xlabel('t');ylabel('x1,x2');title('向前歐拉公式擬合Volterra曲線（30階）');clear%繪圖，清理變量；2）改進歐拉公式函數(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

7、:15;%進入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*(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;%代入改進歐拉公式； 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上

8、的x值；plot(t,x1,'x',t,x2,'x');%描出全部的點值；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,'-k',tt,xx2,'-k','LineWidth',2); %30階擬合曲線；xlabel('t');ylabel('x1,x2');title('改進歐拉公式擬合Volterra曲線

9、（30階）');clear%標注圖形，清理變量。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;%進入150次循環(huán)； k11=x1(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*k2

10、1)*(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*(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);%代入龍

11、格-庫塔公式； 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');%描出全部點值；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,'-g',tt,xx2,'-g','Li

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

13、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.4899 Columns 12 through

14、 21 2.6313 2.8018 3.0072 3.2542 3.5516 3.9100 4.3427 4.8658 5.4987 6.26492）改進歐拉公式Volterra2運行結(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 77.0011 81.6833 86.2278

15、 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運行結(jié)果:步長為：h =0.1000xx1 = Columns 1 through 11 25.0000 27.0680 2

16、9.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五、 實驗結(jié)果分析與實驗小結(jié)：1）結(jié)果分析：為了驗證程序的正確性，我們將上面的三個圖