小行星运动轨迹的Runge-Kutta法模拟

1、小行星运动的Runge-Kutta法模拟一、 背景介绍由于两个恒星作用下行星运动问题没有解析解，只能用数值方法求解微分方程。但是在用一阶近似求解微分方程的时候存在严重的误差累积。当只考虑一个恒星引力影响时的模型如下：.(1)当初始值是时，行星做圆周运动。此时，微分方程的解是。在后面的讨论中，用这个初始条件的方程作为测试方程。如果采用一阶近似，就会有严重的误差累积。如下图所示当行星偏离理想轨道很小的量以后，之后的偏差就会越来越大，直至脱离恒星的束缚。在离散化以后，原来临界稳定的系统变得发散了。二、 用高阶系统去求解单恒星问题当用高于一阶的方法近似求解以上方程时，会取得较好一些的近似。把二阶常微分

2、方程组(1)转化为一阶常微分方程组： ,初始条件是一阶常微分方程组的经典4阶RK法的公式是当时，迭代100000次，模拟行星绕行星圈的轨迹图如下：从上图中可以看出，当模拟绕中心159圈后，轨道的偏移依然很小。为了定量衡量偏差的大小，可以用行星的总能量E=。初始状态时的，经过100000次迭代后总能量变为。可见用4阶KR方法的解精度很高。总能量的偏差量随迭代次数改变的曲线三、 用高阶系统解双恒星问题考虑一种简单情况，即行星初始速度在三个天体所在平面内。行星在两个恒星作用下的微分方程是，其中两个恒星位置是.用经典4阶RK法求解以上微分方程，并且在求解过程中根据行星的速度自适应调整迭代的步长(变动范

3、围是0.0005到0.005之间)。当初值条件为时，迭代5000步后的轨迹图如下：在两个恒星作用下，初始值选的不好时，行星在迭代有限次数后会撞到恒星上去。如以上的初始条件在迭代5780次就会出现行星和恒星的距离小于0.01。当选取初始值为，迭代50000次时的运动轨迹如下：在以上初始值下行星的运动接近周期运动，在上图中行星运行了31周。对以上初始值稍作改动，。迭代35185次时行星与恒星的距离小于0.01。运动轨迹图如下：当初始值改为。迭代34297次时行星与恒星的距离小于0.01。运动轨迹图如下：当初始值改为。迭代50000次的运动轨迹图如下：以上各组测试数据表明，行星在双恒星的引力作用下运

4、动轨迹对初始值很敏感。四、 参考文献吴勃英, 王德明等. 数值分析原理. 北京:科学出版社. 2003：309-310Matlab程序1clear all;clc;close all;%J=-1;L=1f1=(x,vx,y,vy) vx;f2=(x,vx,y,vy) -x/sqrt(x*x+y*y)3);%-(x-L)/sqrt(x-L)*(x-L)+y*y)3)f3=(x,vx,y,vy) vy;f4=(x,vx,y,vy) -y/sqrt(x*x+y*y)3);%-y/sqrt(x-L)*(x-L)+y*y)3)h=0.001;N=10000;X=zeros(1,N);X(1)=1;Vx=

5、zeros(1,N);Vx(1)=0.1;Y=zeros(1,N);Y(1)=1;Vy=zeros(1,N);Vy(1)=0.7;d=0.09;for n=1:N-1 Kx1=f1(X(n),Vx(n),Y(n),Vy(n); Kvx1=f2(X(n),Vx(n),Y(n),Vy(n); Ky1=f3(X(n),Vx(n),Y(n),Vy(n); Kvy1=f4(X(n),Vx(n),Y(n),Vy(n); Kx2=f1(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kvx2=f2(X(n)+h/2*Kx1,Vx(n)+h/

6、2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Ky2=f3(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kvy2=f4(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kx3=f1(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kvx3=f2(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); K

7、y3=f3(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kvy3=f4(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kx4=f1(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3); Kvx4=f2(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3); Ky4=f3(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3)

8、; Kvy4=f4(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3); X(n+1)=X(n)+h/6*(Kx1+2*Kx2+2*Kx3+Kx4); Vx(n+1)=Vx(n)+h/6*(Kvx1+2*Kvx2+2*Kvx3+Kvx4); Y(n+1)=Y(n)+h/6*(Ky1+2*Ky2+2*Ky3+Ky4); Vy(n+1)=Vy(n)+h/6*(Kvy1+2*Kvy2+2*Kvy3+Kvy4); %if X(n+1)*X(n+1)+Y(n+1)*Y(n+1)<d2 | (X(n+1)-L)*(X(n+1)-L)+Y(n+1)*Y(n

9、+1)<d2 %error('d<0.05'); % break; %endendfigure,plot(X,Y);grid,axis equalfigure,plot(X);E=-1./(sqrt(X.*X+Y.*Y)+0.5*(Vx.*Vx+Vy.*Vy);%-2./(sqrt(X-d).*(X-d)+Y.*Y)E0=E(1)figure,plot(E-E(1);程序2clear all;clc;%close all; f1=(x,vx,y,vy) vx;%f2=(x,vx,y,vy,t) -(x-O1x(t)/sqrt(x-O1x(t)*(x-O1x(t)+(

10、y-O1y(t)*(y-O1y(t)3)-(x-O2x(t)/sqrt(x-O2x(t)*(x-O2x(t)+(y-O2y(t)*(y-O2y(t)3);f2=(x,vx,y,vy,t) -(x+1)/sqrt(x+1)*(x+1)+y*y)3)-(x-1)/sqrt(x-1)*(x-1)+y*y)3);f3=(x,vx,y,vy) vy;%f4=(x,vx,y,vy,t) -y/sqrt(x-O1x(t)*(x-O1x(t)+(y-O1y(t)*(y-O1y(t)3)-(y-O2y(t)/sqrt(x-O2x(t)*(x-O2x(t)+(y-O2y(t)*(y-O2y(t)3);f4=(x,

11、vx,y,vy,t) -y/sqrt(x+1)*(x+1)+y*y)3)-y/sqrt(x-1)*(x-1)+y*y)3);h=0.005;t=0;N=50000;X=zeros(1,N);X(1)=0;Vx=zeros(1,N);Vx(1)=-0.349;Y=zeros(1,N);Y(1)=0.1;Vy=zeros(1,N);Vy(1)=1.1;d=0.01;for n=1:N-1 d1=(X(n)-1)*(X(n)-1)+Y(n)*Y(n); d2=(X(n)+1)*(X(n)+1)+Y(n)*Y(n); if d1<d2 | d2<d2 %error('d<0.

12、05'); break; elseif d1<9*d2 | d2<9*d2 h=0.0005; elseif d1<64*d2 | d2<64*d2 h=0.001; else h=0.005; end Kx1=f1(X(n),Vx(n),Y(n),Vy(n); Kvx1=f2(X(n),Vx(n),Y(n),Vy(n),t); Ky1=f3(X(n),Vx(n),Y(n),Vy(n); Kvy1=f4(X(n),Vx(n),Y(n),Vy(n),t); Kx2=f1(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+

13、h/2*Kvy1); Kvx2=f2(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1,t+h/2); Ky2=f3(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kvy2=f4(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1,t+h/2); Kx3=f1(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kvx3=f2(X(n)+h/

14、2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2,t+h/2); Ky3=f3(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kvy3=f4(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2,t+h/2); Kx4=f1(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3); Kvx4=f2(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+

15、h*Kvy3,t+h); Ky4=f3(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3); Kvy4=f4(X(n)+h*Kx3,Vx(n)+h*Kvx3,Y(n)+h*Ky3,Vy(n)+h*Kvy3,t+h); X(n+1) =X(n) +h/6*(Kx1+2*Kx2+2*Kx3+Kx4); Vx(n+1)=Vx(n)+h/6*(Kvx1+2*Kvx2+2*Kvx3+Kvx4); Y(n+1) =Y(n) +h/6*(Ky1+2*Ky2+2*Ky3+Ky4); Vy(n+1)=Vy(n)+h/6*(Kvy1+2*Kvy2+2*Kvy3+Kvy

16、4); t=t+h;endfigure,plot(X(1:n);%E=0.5*(Vx.*Vx+Vy.*Vy);%-2./(sqrt(X-d).*(X-d)+Y.*Y)+1./(sqrt(X.*X+Y.*Y)%E0=E(1)%figure,plot(E);%-E(1)figure,plot(X(1:n),Y(1:n);grid,axis equal程序3clear all;clc;%close all;w=1/sqrt(8);f1=(x,vx,y,vy) vx;%f2=(x,vx,y,vy,t) -(x-O1x(t)/sqrt(x-O1x(t)*(x-O1x(t)+(y-O1y(t)*(y-O1

17、y(t)3)-(x-O2x(t)/sqrt(x-O2x(t)*(x-O2x(t)+(y-O2y(t)*(y-O2y(t)3);f2=(x,vx,y,vy,t) -(x+1)/sqrt(x+1)*(x+1)+y*y)3)-(x-1)/sqrt(x-1)*(x-1)+y*y)3)+(cos(w*t)*x-sin(w*t)*y)/8;f3=(x,vx,y,vy) vy;%f4=(x,vx,y,vy,t) -y/sqrt(x-O1x(t)*(x-O1x(t)+(y-O1y(t)*(y-O1y(t)3)-(y-O2y(t)/sqrt(x-O2x(t)*(x-O2x(t)+(y-O2y(t)*(y-O2y

18、(t)3);f4=(x,vx,y,vy,t) -y/sqrt(x+1)*(x+1)+y*y)3)-y/sqrt(x-1)*(x-1)+y*y)3)+(sin(w*t)*x+cos(w*t)*y)/8;h=0.005;t=0;N=50000;X=zeros(1,N);X(1)=0;Vx=zeros(1,N);Vx(1)=1.1;Y=zeros(1,N);Y(1)=-0.015;Vy=zeros(1,N);Vy(1)=-0.45;d=0.01;for n=1:N-1 d1=(X(n)-1)*(X(n)-1)+Y(n)*Y(n); d2=(X(n)+1)*(X(n)+1)+Y(n)*Y(n); if

19、 d1<d2 | d2<d2 | d1>1000 | d2>1000 %error('d<0.05'); break; elseif d1<9*d2 | d2<9*d2 h=0.0005; elseif d1<64*d2 | d2<64*d2 h=0.001; else h=0.005; end Kx1=f1(X(n),Vx(n),Y(n),Vy(n); Kvx1=f2(X(n),Vx(n),Y(n),Vy(n),t); Ky1=f3(X(n),Vx(n),Y(n),Vy(n); Kvy1=f4(X(n),Vx(n),Y(n

20、),Vy(n),t); Kx2=f1(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kvx2=f2(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1,t+h/2); Ky2=f3(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1); Kvy2=f4(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,Vy(n)+h/2*Kvy1,t+h/2); Kx3=f1(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2); Kvx3=f2(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2,Y(n)+h/2*Ky2,Vy(n)+h/2*Kvy2,t+h/2); Ky3=f3(X(n)+h/2*Kx2,Vx(n)+h/2*Kvx2