交大数值分析大作业.docx

数值分析大作业-第三次校区：徐汇 班级：B学号： 姓名：杨小亮大作业3 对于常微方程数值解问题 检查各种数值算法的长期行为 观察步长对于收敛效果的影响给定方程组1 证明方程组的解是xOy平面上的一个椭圆；2 利用改进的欧拉折线法，4阶标准龙格库塔法，选几个不同的步长h，计算上述方程组的轨道，看看哪种方法和步长能够保持椭圆轨道不变。（计算的时间步要足够多至少10000步）1. 证明：由 求导得：所给方程的特征方程为故是一对共轭复根，所以所求通解为因为所以同理可得：因为即所以，所以整理得：方程组的解是xOy平面上的一个椭圆。2. 研究不同步长下改进的欧拉折现法和四阶标准龙格库塔法轨道变化情况。取a=6,b=3，此时椭圆方程为：。取步数n=15000,分别取不同的步长h=0.0001, 0.001, 0.01, 0.1，利用matlab程序进行计算和绘图。取步长h=0.0001，计算得到的轨道如图 1 所示。此时，改进的欧拉折线法和四阶标准龙格-库塔法轨道仍是椭圆，并且计算结果同精确解几乎完全一样。=图1图2取步长h=0.001，计算得到的轨道如图 3 所示，局部放大如图4所示。此时，改进的欧拉折线法轨道仍是椭圆，而且和精确解几乎一致。但是四阶标准龙格-库塔法轨道出现了明显的向精确解外部发散。=图3图4取步长h=0.01，计算得到的轨道如图 5 所示，局部放大如图6，图7所示。此时，改进的欧拉折线法轨道仍是椭圆，通过进一步局部放大图7可以看出改进的欧拉折线法轨道也出现了微小的向外发散。四阶标准龙格-库塔法轨道出现了严重的向外发散。=图5图6图7取步长h=0.1，计算得到的轨道如图 8 所示。此时，改进的欧拉折线法轨道仍是椭圆，但向外发散程度变大。四阶标准龙格-库塔法不能保持椭圆轨道。=图8结论：通过以上实验可以看出，对于上述常微方程数值解问题的改进的欧拉折线法数值算法和四阶标准龙格-库塔法，在计算的时间步足够多（至少10000步）的情况下：（1） 当步长很小时，两种方法都能够保持椭圆轨道不变，而且都能得到较高的精度（2） 步长增加时，改进的欧拉折线法仍能保持椭圆轨道，只是会出现轻微的向外发散。而四阶标准龙格-库塔法发散情况较严重，不能得到精确解，步长继续增加时，四阶标准龙格-库塔法不能保持椭圆轨道（3） 对于本实验所研究的常微方程数值解问题，欧拉折线法数值算法效果好于四阶标准龙格-库塔法本实验用到的程序：clc; clear; a=6;b=3; n=15000; h=0.1; x1,y1=ImpEuler(a,b,h,n); x2,y2=RungeKutta4(a,b,h,n); hold on;syms u v;p=ezplot(b*u2+a*v2-a*b2); set(p,'color','k','linewidth',4);plot(x1,y1,'r-','linewidth',2) plot(x2,y2,'b-','linewidth',2) legend('精确解','改进Euler', '四阶Runge-Kutta');xlabel('x');ylabel('y');title('h=0.1');hold off;（4） 附：欧拉折线法数值算法和四阶标准龙格-库塔法程序function x,y=ImpEuler(a,b,h,n) x(1)=0; y(1)=b; for i=1:n xp=x(i)+h*a*y(i); yp=y(i)+h*(-b)*x(i); xc=x(i)+h*a*yp; yc=y(i)+h*(-b)*xp; x(i+1)=(xp+xc)/2; y(i+1)=(yp+yc)/2; endfunction x,y=RungeKutta4(a,b,h,n) % x,y为微分方程的数值解% a,b为方程组的系数 % h为步长 % n为迭代步数 x(1)=0; y(1)=b; for i=1:n Kx1=a*y(i); Kx2=a*(y(i)+0.5*h*Kx1); Kx3=a*(y(i)+0.5*h*Kx2); Kx4=a*(y(i)+h*Kx3); x(i+1)=x(i)+h/6*(Kx1+2*Kx2+2*Kx3+Kx4); Ky1=-b*x(i); Ky2=-b*(x(i)+0.5*h*Ky1); Ky3=-b*(x(i)+0.5*h*Ky2); Ky4=-b*(x(i)+h*Ky3); y(i+1)=y(i)+h/6*(Ky1+2*Ky2+2*Ky3+Ky4); end