小行星運(yùn)動(dòng)軌跡的Runge-Kutta法模擬概述

1、小行星運(yùn)動(dòng)的Runge-Kutta法模擬一、背景介紹山于兩個(gè)恒星作用下行星運(yùn)動(dòng)問題沒有解析解，只能用數(shù)值方法求解微分方 程。但是在用一階近似求解微分方程的時(shí)候存在嚴(yán)重的誤差累積。當(dāng)只考慮一個(gè)恒星引力影響時(shí)的模型如下：宀（x2 + y2y（x2 + r）7o=xo,xo=xo丿0 =兒，兒=兒當(dāng)初始值是心=5=0用。=0,兒=1時(shí)，行星做圓周運(yùn)動(dòng)。此時(shí)，微分方 程的解是r = COS(0o在后面的討論中，用這個(gè)初始條件的方程作為測(cè)試方程。 y = sin(/)如果采用-階近似,嚴(yán)+1)"何+衣何,5+1)"何+衣'(“)，就會(huì)有 y(n + l) = y(n) + h

2、y n y n +1) = y S) + hy n)嚴(yán)重的誤差累積。如下圖所示當(dāng)行星偏離理想軌道很小的量以后，之后的偏差就會(huì)越來越大，直至脫離恒 星的束縛。在離散化以后，原來臨界穩(wěn)定的系統(tǒng)變得發(fā)散了。二、用高階系統(tǒng)去求解單恒星問題當(dāng)用高于一階的方法近似求解以上方程時(shí)，會(huì)取得較好一些的近似。 把二階常微分方程組(1)轉(zhuǎn)化為一階常微分方程組：無=1>0=0Vx0=0% = 1（V2 + V2'丿初始條件是V =-V3一階常微分方程組Y'= F（xX）的經(jīng)典4階RK法的公式是丫”+產(chǎn)打+£（仏+瓦+瓦+瓦）K嚴(yán)F（和打）K第+第）2眄+笄+軌）K4=F（xn+h,Yn

3、+hK3）2穴當(dāng)力= 0.01時(shí)，迭代100000次，模擬行星繞行星型里T°°l 159圈的軌跡從上圖中可以看岀，當(dāng)模擬繞中心159圈后，軌道的偏移依然很小。為了定量衡量偏差的大小,可以用行星的總能量E二-一+丄（氣2+叫2）。（2+廠）2初始狀態(tài)時(shí)的£ = -0.5,經(jīng)過100000次迭代后總能量變?yōu)镋 = -0.5-1.4x10-9o 可見用4階KR方法的解精度很高?？偰芰康钠盍侩S迭代次數(shù)改變的曲線三、用高階系統(tǒng)解雙恒星問題考慮一種簡(jiǎn)單情況，即行星初始速度在三個(gè)天體所在平面內(nèi)。行星在兩個(gè)恒 星作用下的微分方程是（x l）（x + l）（A-i）2 + /r

4、（x+i）2+y2pyv< y, =:+:，其中兩個(gè)恒星位置是（-1.0）和（1,0）.（x-i）2 + rf Cv+i）2+rF用經(jīng)典4階RK法求解以上微分方程，并且在求解過程中根據(jù)行星的速度自 適應(yīng)調(diào)整迭代的步長(zhǎng)力（變動(dòng)范圍是0.0005到0.005之間）。當(dāng)初值條件為心=0,兒=0275用。=-03571,幾=1時(shí)，迭代5000步后的軌 跡圖如下:在兩個(gè)恒星作用下，初始值選的不好時(shí)，行星在迭代有限次數(shù)后會(huì)撞到恒星 上去。如以上的初始條件在迭代5780次就會(huì)出現(xiàn)行星和恒星的距離小于0.01 o當(dāng)選取初始值為x0 = 0。=0,x；= -0.349,兒=11,迭代50000次時(shí)的運(yùn)動(dòng)在

5、以上初始值下行星的運(yùn)動(dòng)接近周期運(yùn)動(dòng)，在上圖中行星運(yùn)行了 31周。對(duì)以上初始值稍作改動(dòng),如=0,兒=0*；=-0349 = 111。迭代35185次時(shí)行星與恒星的距離小于0.01 o運(yùn)動(dòng)軌跡圖如下：當(dāng)初始值改為兀=0,兒=0/；=-034&y； = ll o迭代34297次時(shí)行星與恒星的距離小于0.01。運(yùn)動(dòng)軌跡圖如下：-1-0.500.51當(dāng)初始值改為心=0,兒=01用。=-0349,幾=11。迭代50000次的運(yùn)動(dòng)軌跡圖如下:1-0.500.51以上各組測(cè)試數(shù)據(jù)表明，行星在雙恒星的引力作用下運(yùn)動(dòng)軌跡對(duì)初始值很敬 感。四、參考文獻(xiàn)吳勃英，王徳明等.數(shù)值分析原理.北京:科學(xué)出版社.200

6、3： 309-310Matlab 程序 1clear all;clc;close all;%J=-1;L=1f1=0(x,vx,y,vy) vx;(x,vx,y,vy)-x/sqrt(x*x+y*y)A3)(x-L)/sqrc(x-L)*(x-L)+y*y)A3)f3=(x,vx,y,vy)(x,vxzyzvy) h=0001;vy;-y/sqrt(x*x+y*y)A3);%-y/sqrt(x-L)*(x-L)+y*y) A3)N=10000;X=zeros(1ZN);X(1)=1;V>:=zeros (1,N) ;Vx (1) = 0l;Y=zeros(1,N);Y(1)=1;Vy=z

7、eros(1,N);Vy(1)=0 7;for n=l:N-lKxl=f1(X(n),Vx(n),Y(n)zVy(n);Kvxl=f2(X(n)zVx(n),Y(n),Vy(n);Kyl=f3(X(n),Vx(n),Y(n),Vy(n);Kvyl=f4(X(n),Vx(n),Y(n)zVy(n);Kx2=fl(X(n)+h/2*KxlzVx(n)+h/2*Kvxl,Y(n)+h/2*KylzVy(n)+h/2*Kvyl); Kvx2=f2(X(n)+h/2*KxlzVx(n)+h/2*KvxlzY(n)+h/2*Kyl,Vy(n)+h/2*Kvyl);Ky2=f3(X(n)+h/2*Kxl,V

8、x(n)+h/2*KvxlzY(n)+h/2*Kyl,Vy(n)+h/2*Kvyl);Kvy2=f4(X(n)+h/2*Kxl,Vx(n)+h/2*Kvxl,Y(n)+h/2*Kyl,Vy(n)+h/2*Kvyl);Kx3=fl (X (n) +h/2*K>:2, Vx (n) +h/2*Kv>:2, Y (n) +h/2*Ky2z Vy (n) +h/2*Kvy2); Kvx3=f2 (X (n) +h/2*Kx2z Vx (n) +h/2*Kv>:2, Y (n) +h/2*Ky2, Vy (n) +h/2*Kvy2);Ky3=f3 (X (n) +h/2*K>:2

9、z Vx (n) +h/2*Kv>:2, Y (n) +h/2*Ky2z Vy (n) +h/2*Kvy2);Kvy3=f4(X(n)+h/2*Kx2zVx(n)+h/2*Kvx2zY(n)+h/2*Ky2,Vy(n)+h/2*Kvy2);Kx4 = f 1 (X (n) +h*Kx3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3);Kvx4 = f 2 (X (n) +h*Kx3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3);Ky4 = f 3 (X (n) +h*Kx3, V

10、x (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3);Kvy4 = f4 (X (n) +h*K>:3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3z Vy (n) +h*Kvy3);X (n+1) =X (n) +h/6* (Kxl+2*Kx2+2*Kx3+K>:4);Vx(n+1) =Vx(n)+h/6*(Kvxl+2*Kvx2+2*Kvx3+Kvx4);Y(n+1)=Y(n)+h/6*(Kyl+2*Ky2+2*Ky3+Ky4);Vy(n+1)=Vy(n)+h/6*(Kvy1+2 * Kvy2+2 * Kvy3+

11、Kvy4);%if X(n+1)*X(n+l)+Y(n+l)*Y(n+l)<dA2 |(X(n+1)-L)*(X(n+1)-L)+Y(n+1)*Y(n+1)<dA2%error('d<0051);% break;%endendfigure,plot(X,Y);grid,axis equalfigure,plot(X);E=-l/(sqrt(X *X+Y Y)+0 5六(Vxy *Vy);%-2 /(sqrt(X-d)(X-d)+Y.*Y)E0=E (1)figure,plot(E-E (1);程序2clear all;clc;%close all;f1=0(x,vx,

12、y,vy) vx;%f2=(x,vx,y,vyz t)-(x-Olx (t) ) /sqrt ( ( (x-Olx (t) ) * (>：-Olx (t) ) + (y-Oly (t) ) * (y-Oly (t) ) ) A3)-(>：- O2x(t)/sqrt(x-02x(t)*(x-02x(t)+(y-O2y(t)*(y-O2y(t)A3);f2=(x,vxz y,vy,t)_(x+l)/sqrt(x+1)*(x+1)+y*y)A3)_(x-1)/sqrt(x-1)*(x-1)+y*y)A3); f3=(xzvxz y,vy) vy;%f4=3(x,vxz y,vy,t)-y

13、/sqrt(x-Olx(t)*(x-Olx(t) + (y-Oly(t)*(y-Oly(t)A3)-(y-O2y(t) )/s qrt ( ( (x-02x(t)*(x-02x(t) + (y-O2y(t)*(y-O2y(t) ) A3);f4=(x, vxz y,vyz t)-y/sqrt(x+1)*(x+l)+y*y)A3)-y/sqrt(x-1)*(x-l)+y*y) A3);h=0005;t=0;N=50000;X=zeros(1ZN);X(l)=0;Vx=zeros(1,N);Vx(1)=一0349;Y=zeros(1ZN);Y(1)=0l;Vy=zeros(1z N);Vy(1)=

14、11;d=001;for n=l:N-ldl=(X(n) -1) * (X(n)-1)+Y(n) *Y(n);d2=(X(n)+l)*(X(n)+1)+Y(n)*Y(n);if dl<dA2 | d2<dA2%error (P<005J ;break;elseif dl<9*dA2 | d2<9*dA2h=00005;elseif dl<64*dA2 | d2<64*dA2h=0001;elseh=0005;endKxl=f1(X(n),Vx(n),Y(n),Vy(n);Kvxl=f2(X(n),Vx(n),Y(n),Vy(n),t);Kyl=f3(X

15、(n)zVx(n),Y(n),Vy(n);Kvyl = f 4 (X(n) zVx(n) , Y(n) zVy (n) z t);Kx2 = fl (X (n) +h/2*Kxlz Vx (n) +h/2 *Kvxl z Y (n) +h/2 *Kyl z Vy (n) +h/2*Kvyl);Kvx2 = f2 (X (n) +h/2*Kxl, Vx (n) +h/2*Kvxl, Y (n) +h/2*Kylz Vy (n) +h/2*Kvyl, t+h/2 )；Ky2=f3(X(n)+h/2*KxlzVx(n)+h/2*Kvxl,Y(n)+h/2*Kyl,Vy(n)+h/2*Kvyl);Kv

16、y2 = f4(X(n)+h/2 *Kxl,Vx(n)+h/2*Kvxl,Y(n)+h/2*Kyl,Vy(n)+h/2*Kvyl,t+h/2 )；Kx3=fl (X (n) +h/2*K>:2z Vx (n) +h/2*Kv>:2, Y (n) +h/2*Ky2z Vy (n) +h/2*Kvy2);Kvx3=f2 (X (n) +h/2*Kx2, Vx (n) +h/2*Kv>:2z Y (n) +h/2*Ky2, Vy (n) +h/2*Kvy2, t+h/2 )；Ky3=f3 (X (n) +h/2*K>:2z Vx (n) +h/2*Kv>:2, Y (n

17、) +h/2*Ky2, Vy (n) +h/2*Kvy2);Kvy3=f4 (X (n) +h/2*K>:2z Vx (n) +h/2*Kv>:2z Y (n) +h/2*Ky2, Vy (n) +h/2*Kvy2, t+h/2 )；Kx4 = f 1 (X (n) +h*K>:3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3);Kvx4 = f2 (X (n) +h*K>:3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3z t+h);Ky4 = f 3 (X (

18、n) +h*Kx3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3);Kvy4 = f 4 (X (n) +h*Kx3, Vx (n) +h*Kv>:3, Y (n) +h*Ky3, Vy (n) +h*Kvy3z t+h);X (n+1) =X(n) +h/6*(Kxl+2*Kx2+2*Kx3+Kx4);Vx(n+1)=Vx(n)+h/6*(Kvxl+2*Kvx2+2*Kvx3+Kvx4);Y (n+1) =Y (n) +h/6*(Kyl+2*Ky2+2*Ky3+Ky4);Vy(n+1)=Vy(n)+h/6*(Kvyl+2 *Kvy2

19、+2 *Kvy3+Kvy4);+h;endfigure,plot(X(1:n);%E=05(Vx水Vk+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=l/sqrt(8);f1=0(x,vx,y,vy) vx;務(wù)f2=(x,vx,y,vy,t)_ (x-Olx (t) ) /sqrt ( ( (x-01x (t) ) * (x-01>

20、: (t) ) + (y-Oly (t) ) * (y-Oly (t)人3) _ (x_ 02x(t)/sqrt(x-02x(t)*(x-02x(t)+(y-O2y(t)*(y-O2y(t)A3);f2=(x,vxz y,vy,t)_ (x+1) /sqrt ( ( (x+1) * (>:+l) +y*y) A3) _ (x-1) /sqrt ( ( (x-1) * (>：-l) +y*y) A3) + (cos (w*t)*x-sin(w*t)*y) /8;f3= (>：/vxz y, vy) vy;(>：/vxz y, vyz t)-y/sqrt(x-01x(t)*

21、(x-Olx(t)+(y-Oly(t)*(y-Oly(t)A3)-(y-O2y(t)/s qrt ( ( (x-02x(t)*(x-02x(t) + (y-O2y(t)*(y-O2y(t) ) ) A3);(x, vxz y,vy,t)-y/sqrt ( ( (x+1) * (x+1) +y*y) A3) -y/sqrt ( ( (x-1) * (x-1) +y*y) A3) + (sin (w*t) *>:+ cos(w*t)*y)/8;h=0005;t=0;N=50000;X=zeros(1,N);X(l)=0;Vx=zeros(1,N);Vx(1)=11;Y=zeros(1ZN);

22、Y(1)=-0015;Vy=zeros(1,N);Vy(1)=-0 45;d=001;for n=l:N-ldl=(X(n) -1) * (X(n)-1)+Y(n) *Y(n);d2=(X(n)+l)*(X(n)+1)+Y(n)*Y(n);if dl<dA2 | d2<dA2 | dl>1000 | d2>1000%error(1d<0051);break;elseif dl<9*dA2 | d2<9*dA2h=00005;elseif dl<64*dA2 | d2<64*dA2h=0001;elseh=0005;endKxl=f1(X(n

23、),Vx(n),Y(n),Vyn);Kvxl=f2(X(n)zVx(n),Y(n)zVy(n),t);Kyl=f3(X(n)fVx(n)zY(n),Vy(n);Kvyl = f4 (X(n) zVx(n) z Y(n) ,Vy (n) z t);Kx2=fl (X (n) +h/2*Kxl, Vx (n) +h/2 *Kvxl z Y (n) +h/2 *Kyl, Vy (n) +h/2*Kvyl);Kvx2 = f2 (X (n) +h/2*Kxl, Vx (n) +h/2*Kvxlz Y (n) +h/2*Kylz Vy (n) +h/2*Kvyl, t+h/2 )；Ky2=f3 (X (

24、n) +h/2*Kxlz Vx (n) +h/2 *Kvxl z Y (n) +h/2 *Kyl z Vy (n) +h/2*Kvyl);Kvy2 = f4(X(n)+h/2*Kxl,Vx(n)+h/2 *Kvxl,Y(n)+h/2*KylzVy(n)+h/2*Kvyl,t+h/2 )；Kx3=fl (X (n) +h/2*Kx2z Vx (n) +h/2*Kv>:2z Y (n) +h/2*Ky2, Vy (n) +h/2*Kvy2);Kvx3=f2 (X(n) +h/2*Kx2z Vx (n) +h/2*Kv>:2z Y (n) +h/2 *Ky2, Vy (n) +h/2 *Kvy2, t+h/2 )；Ky3=f3 (X (n) +h/2*K>:2z Vx (n)