计算方法大作业(第三次).doc

大作业3 对于常微方程数值解问题检查各种数值算法的长期行为观察步长对于收敛效果的影响给定方程组1 证明方程组的解是xOy平面上的一个椭圆；2 利用改进的欧拉折线法，4阶标准龙格库塔法，选几个不同的步长h，计算上述方程组的轨道，看看哪种方法和步长能够保持椭圆轨道不变。（计算的时间步要足够多至少10000步）1. 证明：因为所以所给方程的特征方程为所以是一对共轭复根，所以所求通解为因为所以同理可得：因为即所以，所以所以方程组的解是xOy平面上的一个椭圆。2.设a=4,b=1;（1）用改进的欧拉折线法计算，程序及结果如下：a=4;b=1;n=10000; H=6*pi; x,y=eulermend(a,b,H,n);plot(x,y)图1 h=6*pi/10000时的椭圆轨道图像 H=7*pi; x,y=eulermend(a,b,H,n); plot(x,y,g-)图2 h=7*pi/10000时的椭圆轨道图像 H=10*pi; x,y=eulermend(a,b,H,n); plot(x,y,r-)图3 h=10*pi/10000时的椭圆轨道图像 H=20*pi; x,y=eulermend(a,b,H,n); plot(x,y)图4 h=20*pi/10000时的椭圆轨道图像 H=30*pi; x,y=eulermend(a,b,H,n); figure(2),plot(x,y,r-)图5 h=30*pi/10000时的椭圆轨道图像从上述结果可知，改进的欧拉折线法对于不同步长能够保持椭圆轨道不变。以下是欧拉改进方法函数eulermend.mfunction x,y=eulermend(a,b,H,n)h=H/n;x(1)=0;y(1)=b;for i=1:n x1=x(i)+h*a*y(i); y1=y(i)+h*(-b)*x(i); x2=x(i)+h*a*y1; y2=y(i)+h*(-b)*x1; x(i+1)=(x1+x2)/2; y(i+1)=(y1+y2)/2;end（2）用4阶标准龙格库塔法计算，程序及结果如下： a=4;b=1; n=10000; H=pi; x,y=ronger4(a,b,H,n); plot(x,y)图6 h=1*pi/10000时的椭圆轨道图像 H=2*pi; x,y=ronger4(a,b,H,n); plot(x,y,r-)图7 h=2*pi/10000时的椭圆轨道图像 H=6*pi; x,y=ronger4(a,b,H,n); plot(x,y)图8 h=6*pi/10000时的椭圆轨道图像 H=7*pi; x,y=ronger4(a,b,H,n); plot(x,y,g-)图9 h=7*pi/10000时的椭圆轨道图像 H=10*pi; x,y=ronger4(a,b,H,n); plot(x,y,r-)图10 h=10*pi/10000时的椭圆轨道图像从以上结果可知，随着步长越来越大，用4阶标准龙格库塔法计算的椭圆轨道变化也越来越大，只有当步长比较小时，如图7，图8所示中步长取得相对较小，椭圆轨道才保持不变。以下是4阶标准龙格库塔法程序ronger4.m：function x,y=ronger4(a,b,H,n)h=H/n;x(1)=0;y(1)=b;for i=1:n K1=a*y(i); K2=a*(y(i)+0.5*h*K1); K3=a*(y(i)+0.5*h*K2); K4=a*(y(i)+h*K3); x(i+1)=x(i)+h/6*(K1+2*K2+2*K3+K4); R1=-b*x(i); R2=-b*(x(i)+0.5*h*R1);