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