小行星運動軌跡的Runge-Kutta法模擬

1、1 / 11 小行星運動的runge-kutta 法模擬一、背景介紹由于兩個恒星作用下行星運動問題沒有解析解，只能用數(shù)值方法求解微分方程。但是在用一階近似求解微分方程的時候存在嚴重的誤差累積。當只考慮一個恒星引力影響時的模型如下： .(1)當初始值是00001,0,0,1xyxy時，行星做圓周運動。此時，微分方程的解是cos( )sin( )xtyt。在后面的討論中，用這個初始條件的方程作為測試方程。如果采用一階近似，(1)( )( ), (1)( )( )(1)( )( ),(1)( )( )x nx nhx nx nx nhx ny ny nhy ny ny nhy n，就會有嚴重的誤差累

2、積。如下圖所示當行星偏離理想軌道很小的量以后，之后的偏差就會越來越大， 直至脫離恒星的束縛。在離散化以后，原來臨界穩(wěn)定的系統(tǒng)變得發(fā)散了。二、用高階系統(tǒng)去求解單恒星問題當用高于一階的方法近似求解以上方程時，會取得較好一些的近似。把二階常微分方程組 (1)轉(zhuǎn)化為一階常微分方程組：3222322200000000()(),xxxyyyxyxxxxyyyy2 / 11 32223222()()xxyyxvxvxyyvyvxy,初始條件是00001001xyxyvv一階常微分方程組00( ,)()xxyfyyy的經(jīng)典 4 階 rk 法的公式是112341213243()6(,)(,)22(,)22(,)

3、nnnnnnnnnnhxhhxhhxxhhyykkkkkfykfykkfykkfyk當0.01h時，迭代 100000次，模擬行星繞行星100000*0.011592圈的軌跡圖如下：從上圖中可以看出，當模擬繞中心159 圈后，軌道的偏移依然很小。為了定量衡量偏差的大小， 可以用行星的總能量e=222211()()2xyvvxy。3 / 11 初始狀態(tài)時的0.5e，經(jīng)過 100000次迭代后總能量變?yōu)?0.51.410e。可見用 4 階 kr 方法的解精度很高。總能量的偏差量隨迭代次數(shù)改變的曲線三、用高階系統(tǒng)解雙恒星問題考慮一種簡單情況， 即行星初始速度在三個天體所在平面。行星在兩個恒星作用下的

4、微分方程是332222223322222200000000(1)(1)(1)(1)(1)(1), ,xxxxyxyyyyxyxyxx xxyyyy，其中兩個恒星位置是( 1,0)1 0和（，）.用經(jīng)典 4 階 rk 法求解以上微分方程， 并且在求解過程中根據(jù)行星的速度自適應調(diào)整迭代的步長 h(變動圍是 0.0005 到 0.005 之間 )。當初值條件為00000,0.275,0.3571,1xyxy時，迭代 5000 步后的軌跡圖如下：4 / 11 在兩個恒星作用下， 初始值選的不好時， 行星在迭代有限次數(shù)后會撞到恒星上去。如以上的初始條件在迭代5780 次就會出現(xiàn)行星和恒星的距離小于0.0

5、1。當選取初始值為00000,0,0.349,1.1xyxy，迭代 50000 次時的運動軌跡如下：在以上初始值下行星的運動接近周期運動，在上圖中行星運行了31 周。5 / 11 對以上初始值稍作改動，00000,0,0.349,1.11xyxy。迭代35185次時行星與恒星的距離小于0.01。運動軌跡圖如下：當初始值改為00000,0,0.348,1.1xyxy。迭代 34297 次時行星與恒星的距離小于 0.01。運動軌跡圖如下：當初始值改為00000,0.1, 0.349,1.1xyxy。迭代 50000 次的運動軌6 / 11 跡圖如下：以上各組測試數(shù)據(jù)表明， 行星在雙恒星的引力作用下

6、運動軌跡對初始值很敏感。四、參考文獻吳勃英 , 王德明等 .數(shù)值分析原理 .:科學.2003：309-310 matlab 程序 1 clear 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=

7、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-17 / 11 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

8、(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*kv

9、y2); 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

10、*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);%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(

11、n+1)d2%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);程序 2 clear all;clc;%close all;f1=(x,vx,y,vy) vx;%f2=(x,vx,y,vy,t) 8 / 11 -(x-o1x(t)/sqrt(x-o1x(t)*(x-o1x(t)+(y-o1y(t)*(y-

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

13、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;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*d

14、2 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,y(n)+h/2*k

15、y1,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);9 / 11 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

16、/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*ky

17、3,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;endfigure,plot(x(1:n);%e=0.5*(vx.*vx+vy.*vy

18、);%-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程序 3 clear 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)*

19、(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;10 / 11 %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

20、*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;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;e

21、lseif d19*d2 | 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)

22、+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);11 / 11 kvx3=f2(x(n)+h/2*kx2,vx(n)+h/2*kvx2,y(n)+h/2*ky2,vy(n)+