三种方法解决voterra方程参考模板

1、东南大学电气工程学院MATLAB数学建模 实验报告三种方法解析Volterra方程刘海东（16006213）2008/12/131 / 12一、 实验目的：通过MATLAB实现书上含初值一阶Volterra方程组的解析，在解析过程之中注意加深对三种方法：向前欧拉公式、改进的欧拉公式和四阶龙格-库塔公式的理解，注意各个方法的细节和精度比较。进一步熟悉MATLAB编程。二、 实验原理：在数学建模理论课上，我们在书上6.4节中介绍了Volterra模型，并用相平面分析法对问题作了简要分析。随后我们又在6.7节中学习到了常微分方程数值解，学习了欧拉方法和龙格-库塔方法，所以我们尝试借助6.7节中的数值

2、方法来考虑如下Volterra方程：为了跟书上的解析图相对应，我选取区间；步长=0.1。1) 向前欧拉公式(1阶精度)：我们先根据步长将区间分为150等分，在其中任一区间上（）取对应左端点的和作为递推和的初值，将其带入到向前欧拉公式组得到如下方程组：这样，我们就可以根据初值层层递推出对应后面150个t值的x值。具体在MATLAB中的实现过程：首先将各自初值赋予初始变量，并将步长输入。根据步长，在取定循环值t=0.1:0.1:15形成一个150次的循环，循环过程为：根据上面方程组的前两行带入初值求解下一次数值，然后将下次数值作为初值不停迭代下去，即得到150个要求的值。因为总数据太过庞大，我数据

3、输出只输出了t在0,2上的数据，为使结果更为形象直观，我在把全部结果其标注在图形上，并随后进行30阶拟合成Volterra解析曲线。2) 改进欧拉公式(2阶精度)：为在提高精度的同时不使计算复杂，改进的欧拉公式将向前欧拉公式的简单与梯形公式的高精度结合在一起。它把迭代过程简化成两步：先由向前欧拉公式算出x(k)的预测值，再把预测值代入梯形公式的右端进行校正，得到如下方程组：分组递推迭代过程跟向前欧拉公式一样。具体在MATLAB中实现过程也跟向前欧拉公式类似，只不过随着方程的改变而每个循环中间多做了计算而已。3) 4阶龙格-库塔公式(4阶精度)：由上面两个公式得到启发，取四个点的斜率值作加权平均

4、来作为斜率的近似值，这样我们就提高了精度。所以我们在每一小区间上取四个点构造得到下面4阶龙格-库塔公式：分组递推迭代过程同上。而在MATLAB中实现过程依然和上面类似，只不过相比改进欧拉公式每次循环中计算过程更为复杂而已。三、 实验源代码：1）向前欧拉公式函数文件Volterra1:function Volterra1%定义函数；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次循环； 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对数值的输出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）改进欧拉公式函数文件Volterra2:function Volterra2%定义函数；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次循环； 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对数值的输出；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）龙格-库塔公式函数文件Volterra4:function Volterra4%定义函数；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次循环； 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对数值的输出；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%标注图形，清理变量。四、 实验结果：因为本实验的数据实在太多，我只取了前21对数值进行显示，而为验证实验的正确性，我描出了所有的点值，并采用30阶拟合逼近得出各种方法下所算出的解析曲线从而可以与书上的理论曲线比对来验证。1）向前欧拉公式Volterra1运行结果:步长为：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运行结果:步长为：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运行结果:步长为：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五、 实验结果分析与实验小结：1）结果分析：为了验证程序的正确性，我们将上面的三个图