小行星运动轨迹的RungeKutta法模拟Word文档格式.docx
《小行星运动轨迹的RungeKutta法模拟Word文档格式.docx》由会员分享，可在线阅读，更多相关《小行星运动轨迹的RungeKutta法模拟Word文档格式.docx（13页珍藏版）》请在冰豆网上搜索。
二、用高阶系统去求解单恒星问题
当用高于一阶的方法近似求解以上方程时，会取得较好一些的近似。
把二阶常微分方程组
（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）<
d^2||（X（n+1）-L）*（X（n+1）-L）+Y（n+1）*Y（n+1）<
d^2
%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
%closeall;
%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）;
%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
（1）=0;
Vx
（1）=-0.349;
Y
（1）=0.1;
Vy
（1）=1.1;
d=0.01;
d1=（X（n）-1）*（X（n）-1）+Y（n）*Y（n）;
d2=（X（n）+1）*（X（n）+1）+Y（n）*Y（n）;
ifd1<
d^2||d2<
break;
elseifd1<
9*d^2||d2<
9*d^2
h=0.0005;
64*d^2||d2<
64*d^2
h=0.001;
else
h=0.005;
end
Kvx1=f2（X（n）,Vx（n）,Y（n）,Vy（n）,t）;
Kvy1=f4（X（n）,Vx（n）,Y（n）,Vy（n）,t）;
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）;
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）;
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）;
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）;
Kvx4=f2（X（n）+h*Kx3,Vx（n）+h*Kvx3,Y（n）+h*Ky3,Vy（n）+h*Kvy3,t+h）;
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）;
Y（n+1）=Y（n）+h/6*（Ky1+2*Ky2+2*Ky3+Ky4）;
t=t+h;
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）
n）,Y（1:
程序3
w=1/sqrt（8）;
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;
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;
Vx
（1）=1.1;
Y
（1）=-0.015;
Vy
（1）=-0.45;
d^2||d1>
1000||d2>
1000
%figure,plot（X（1:
升级会员 