三种方法解决voterra方程.doc

《三种方法解决voterra方程.doc》由会员分享，可在线阅读，更多相关《三种方法解决voterra方程.doc（11页珍藏版）》请在课桌文档上搜索。

1、 一、 实验目的：通过MATLAB实现书上含初值一阶Volterra方程组的解析，在解析过程之中注意加深对三种方法：向前欧拉公式、改良的欧拉公式和四阶龙格-库塔公式的理解，注意各个方法的细节和精度比拟。进一步熟悉MATLAB编程。二、 实验原理：在数学建模理论课上，我们在书上6.4节中介绍了Volterra模型，并用相平面分析法对问题作了简要分析。随后我们又在6.7节中学习到了常微分方程数值解，学习了欧拉方法和龙格-库塔方法，所以我们尝试借助6.7节中的数值方法来考虑如下Volterra方程：为了跟书上的解析图相对应，我选取区间；步长=0.1。1) 向前欧拉公式(1阶精度)：我们先根据步长将区

2、间分为150等分，在其中任一区间上取对应左端点的和作为递推和的初值，将其带入到向前欧拉公式组得到如下方程组：这样，我们就可以根据初值层层递推出对应后面150个t值的x值。具体在MATLAB中的实现过程：首先将各自初值赋予初始变量，并将步长输入。根据步长，在取定循环值t=0.1:0.1:15形成一个150次的循环，循环过程为：根据上面方程组的前两行带入初值求解下一次数值，然后将下次数值作为初值不停迭代下去，即得到150个要求的值。因为总数据太过庞大，我数据输出只输出了t在0,2上的数据，为使结果更为形象直观，我在把全部结果其标注在图形上，并随后进展30阶拟合成Volterra解析曲线。2) 改良

3、欧拉公式(2阶精度)：为在提高精度的同时不使计算复杂，改良的欧拉公式将向前欧拉公式的简单与梯形公式的高精度结合在一起。它把迭代过程简化成两步：先由向前欧拉公式算出x(k)的预测值，再把预测值代入梯形公式的右端进展校正，得到如下方程组：分组递推迭代过程跟向前欧拉公式一样。具体在MATLAB中实现过程也跟向前欧拉公式类似，只不过随着方程的改变而每个循环中间多做了计算而已。3) 4阶龙格-库塔公式(4阶精度)：由上面两个公式得到启发，取四个点的斜率值作加权平均来作为斜率的近似值，这样我们就提高了精度。所以我们在每一小区间上取四个点构造得到下面4阶龙格-库塔公式：分组递推迭代过程同上。而在MATLAB

4、中实现过程依然和上面类似，只不过相比改良欧拉公式每次循环中计算过程更为复杂而已。三、 实验源代码：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)*(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

5、)=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,-r,tt,xx2,-r);%30阶拟合曲线xlabel(t);ylabel(x1,x2);title(向前欧拉公式拟合Volterra曲线30阶);clear%绘图，清理变量；2改良

6、欧拉公式函数文件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: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)

7、-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上的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阶拟合曲线；

8、xlabel(t);ylabel(x1,x2);title(改良欧拉公式拟合Volterra曲线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(

9、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*(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)+

10、h/6*(k21+2*k22+2*k23+k24);%代入龙格-库塔公式； 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,LineWidth,2); %30阶拟合曲线；

11、xlabel(t);ylabel(x1,x2);title(改良欧拉公式拟合Volterra曲线30阶);clear%标注图形，清理变量。四、 实验结果：因为本实验的数据实在太多，我只取了前21对数值进展显示，而为验证实验的正确性，我描出了所有的点值，并采用30阶拟合逼近得出各种方法下所算出的解析曲线从而可以与书上的理论曲线比对来验证。1向前欧拉公式Volterra1运行结果:步长为：xx1 = Columns 1 through 11 Columns 12 through 21xx2 = Columns 1 through 11 Columns 12 through 212改良欧拉公式Vol

12、terra2运行结果:步长为：xx1 = Columns 1 through 11 Columns 12 through 21xx2 = Columns 1 through 11 Columns 12 through 213龙格-库塔公式Volterra4运行结果:步长为：xx1 = Columns 1 through 11 Columns 12 through 21xx2 = Columns 1 through 11 Columns 12 through 21五、 实验结果分析与实验小结：1结果分析：为了验证程序的正确性，我们将上面的三个图分别与书上图6.16比对发现：大体轮廓根本类似，但仔细分辨:向前欧拉公式的精度明显比拟低，而改良欧拉公式已经逼近书上曲线，而书上曲线即是四阶龙格-库塔的解曲线，所以最后一个图与书上几乎一样，证明了实验的正确性。2实验小结：通过本实验的操作，我又把教师上课所讲的三种典型方法的细节和操作过程复习了解了一遍，把上课并未熟悉的4阶龙格-库塔公式演练了一遍，当然这都得益于MATLAB强大的计算功能。通过这次实验，我除了以上的收获之外，还提高了自己编程画图的能力，在控制图形的样式方面有了改良，且回顾加强了拟合曲线的能力。11 / 11