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

大作业3 对于常微方程数值解问题 检查各种数值算法的长期行为 观察步长对于收敛效果的影响 给定方程组 证明方程组的解是xOy平面上的一个椭圆； 利用①改进的欧拉折线法，②4阶标准龙格－库塔法，选几个不同的步长h，计算上述方程组的轨道，看看哪种方法和步长能够保持椭圆轨道不变。（计算的时间步要足够多――至少10000步） 证明： 因为 所以 所给方程的特征方程为 所以是一对共轭复根， 所以所求通解为 因为所以 同理可得： 因为 即 所以， 所以 所以 方程组的解是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.m function [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); R3=-b*(x(i)+0.5*h*R2); R4=-b*(x(i)+h*R3); y(i+1)=y(i)+h/6*(R1+2*R2+2*R3+R4); end