小行星運動軌跡的Runge-Kutta法模擬[相關(guān)材料]

1、小行星運動的Runge-Kutta法模擬一、 背景介紹由于兩個恒星作用下行星運動問題沒有解析解，只能用數(shù)值方法求解微分方程。但是在用一階近似求解微分方程的時候存在嚴(yán)重的誤差累積。當(dāng)只考慮一個恒星引力影響時的模型如下：.(1)當(dāng)初始值是時，行星做圓周運動。此時，微分方程的解是。在后面的討論中，用這個初始條件的方程作為測試方程。如果采用一階近似，就會有嚴(yán)重的誤差累積。如下圖所示當(dāng)行星偏離理想軌道很小的量以后，之后的偏差就會越來越大，直至脫離恒星的束縛。在離散化以后，原來臨界穩(wěn)定的系統(tǒng)變得發(fā)散了。二、 用高階系統(tǒng)去求解單恒星問題當(dāng)用高于一階的方法近似求解以上方程時，會取得較好一些的近似。把二階常微分

2、方程組(1)轉(zhuǎn)化為一階常微分方程組： ,初始條件是一階常微分方程組的經(jīng)典4階RK法的公式是當(dāng)時，迭代100000次，模擬行星繞行星圈的軌跡圖如下：從上圖中可以看出，當(dāng)模擬繞中心159圈后，軌道的偏移依然很小。為了定量衡量偏差的大小，可以用行星的總能量E=。初始狀態(tài)時的，經(jīng)過100000次迭代后總能量變?yōu)椤？梢娪?階KR方法的解精度很高?？偰芰康钠盍侩S迭代次數(shù)改變的曲線三、 用高階系統(tǒng)解雙恒星問題考慮一種簡單情況，即行星初始速度在三個天體所在平面內(nèi)。行星在兩個恒星作用下的微分方程是，其中兩個恒星位置是.用經(jīng)典4階RK法求解以上微分方程，并且在求解過程中根據(jù)行星的速度自適應(yīng)調(diào)整迭代的步長(變動范

3、圍是0.0005到0.005之間)。當(dāng)初值條件為時，迭代5000步后的軌跡圖如下：在兩個恒星作用下，初始值選的不好時，行星在迭代有限次數(shù)后會撞到恒星上去。如以上的初始條件在迭代5780次就會出現(xiàn)行星和恒星的距離小于0.01。當(dāng)選取初始值為，迭代50000次時的運動軌跡如下：在以上初始值下行星的運動接近周期運動，在上圖中行星運行了31周。對以上初始值稍作改動，。迭代35185次時行星與恒星的距離小于0.01。運動軌跡圖如下：當(dāng)初始值改為。迭代34297次時行星與恒星的距離小于0.01。運動軌跡圖如下：當(dāng)初始值改為。迭代50000次的運動軌跡圖如下：以上各組測試數(shù)據(jù)表明，行星在雙恒星的引力作用下運

4、動軌跡對初始值很敏感。四、 參考文獻吳勃英, 王德明等. 數(shù)值分析原理. 北京:科學(xué)出版社. 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+1)d

9、2 %error(d0.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)+(y-O1y(t)*(y-O1y(t)3)-(

10、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

11、)*(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 d1d2 | d2d2 %error(d0.05); break; elseif d19*d2 | d29*d2 h=0.

12、0005; elseif d164*d2 | d264*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)+h/2*Kvy1); Kvx2=f2(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,Y(n)+h/2*Ky1,

13、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

14、=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)

15、+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;endfigure,plot(X(1:n);%E=0.5*(Vx.*Vx+Vy.*Vy);%-2./

16、(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-O1y(t)3)-(x-O2x(t)/sqrt(x-O2x(t)*(x-O2x(t)+(y-O2y(t)*(y-O2y(t)

17、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)

18、*(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 d1d2 | d21000 | d21000 %error(d0.05); break; elseif d19*d2

19、| d29*d2 h=0.0005; elseif d164*d2 | d264*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)+h/2*Kvy1); Kvx2=f2(X(n)+h/2*Kx1,Vx(n)+h/2*Kvx1,

20、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