小行星运动轨迹的RungeKutta法模拟.docx
《小行星运动轨迹的RungeKutta法模拟.docx》由会员分享，可在线阅读，更多相关《小行星运动轨迹的RungeKutta法模拟.docx（17页珍藏版）》请在冰点文库上搜索。
小行星运动轨迹的RungeKutta法模拟
小行星运动的Runge-Kutta法模拟
一、背景介绍
由于两个恒星作用下行星运动问题没有解析解，只能用数值方法求解微分方程。
但是在用一阶近似求解微分方程的时候存在严重的误差累积。
当只考虑一个恒星引力影响时的模型如下：
……..
（1）
当初始值是
时，行星做圆周运动。
此时，微分方程的解是
。
在后面的讨论中，用这个初始条件的方程作为测试方程。
如果采用一阶近似，
，就会有严重的误差累积。
如下图所示
当行星偏离理想轨道很小的量以后，之后的偏差就会越来越大，直至脱离恒星的束缚。
在离散化以后，原来临界稳定的系统变得发散了。
二、用高阶系统去求解单恒星问题
当用高于一阶的方法近似求解以上方程时，会取得较好一些的近似。
把二阶常微分方程组
（1）转化为一阶常微分方程组：
初始条件是
一阶常微分方程组
的经典4阶RK法的公式是
当
时，迭代100000次，模拟行星绕行星
圈的轨迹图如下：
从上图中可以看出，当模拟绕中心159圈后，轨道的偏移依然很小。
为了定量衡量偏差的大小，可以用行星的总能量E=
。
初始状态时的
，经过100000次迭代后总能量变为
。
可见用4阶KR方法的解精度很高。
总能量的偏差量随迭代次数改变的曲线
三、用高阶系统解双恒星问题
考虑一种简单情况，即行星初始速度在三个天体所在平面。
行星在两个恒星作用下的微分方程是
，其中两个恒星位置是
.
用经典4阶RK法求解以上微分方程，并且在求解过程中根据行星的速度自适应调整迭代的步长
（变动围是0.0005到0.005之间）。
当初值条件为
时，迭代5000步后的轨迹图如下：
在两个恒星作用下，初始值选的不好时，行星在迭代有限次数后会撞到恒星上去。
如以上的初始条件在迭代5780次就会出现行星和恒星的距离小于0.01。
当选取初始值为
，迭代50000次时的运动轨迹如下：
在以上初始值下行星的运动接近周期运动，在上图中行星运行了31周。
对以上初始值稍作改动，
。
迭代35185次时行星与恒星的距离小于0.01。
运动轨迹图如下：
当初始值改为
。
迭代34297次时行星与恒星的距离小于0.01。
运动轨迹图如下：
当初始值改为
。
迭代50000次的运动轨迹图如下：
以上各组测试数据表明，行星在双恒星的引力作用下运动轨迹对初始值很敏感。
四、参考文献
吴勃英,王德明等.数值分析原理.:
科学.2003：
309-310
Matlab程序1
clearall;clc;closeall;
;%J=-1;L=1
f1=（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=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;
forn=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/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）;
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）;
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）;
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）;
%ifX（n+1）*X（n+1）+Y（n+1）*Y（n+1）%error（'d<0.05'）;
%break;
%end
end
figure,plot（X,Y）;grid,axisequal
figure,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））;
程序2
clearall;clc;%closeall;
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-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,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;
forn=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）;
ifd1%error（'d<0.05'）;
break;
elseifd1<9*d^2||d2<9*d^2
h=0.0005;
elseifd1<64*d^2||d2<64*d^2
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）+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,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）+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+Kvy4）;
t=t+h;
end
figure,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,axisequal
程序3
clearall;clc;%closeall;
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-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）+（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（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;
forn=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）;
ifd11000||d2>1000
%error（'d<0.05'）;
break;
elseifd1<9*d^2||d2<9*d^2
h=0.0005;
elseifd1<64*d^2||d2<64*d^2
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）+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,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）+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+Kvy4）;
t=t+h;
end
%figure,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,axisequal
升级会员 