小行星运动轨迹的RungeKutta法模拟Word文档格式.docx
《小行星运动轨迹的RungeKutta法模拟Word文档格式.docx》由会员分享,可在线阅读,更多相关《小行星运动轨迹的RungeKutta法模拟Word文档格式.docx(13页珍藏版)》请在冰豆网上搜索。
![小行星运动轨迹的RungeKutta法模拟Word文档格式.docx](https://file1.bdocx.com/fileroot1/2023-1/25/aa14a7ae-fcd7-4814-b9f3-8483e1365f76/aa14a7ae-fcd7-4814-b9f3-8483e1365f761.gif)
二、用高阶系统去求解单恒星问题
当用高于一阶的方法近似求解以上方程时,会取得较好一些的近似。
把二阶常微分方程组
(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: