1、1 / 11 小行星運動的runge-kutta 法模擬一、背景介紹由于兩個恒星作用下行星運動問題沒有解析解，只能用數值方法求解微分方程。但是在用一階近似求解微分方程的時候存在嚴重的誤差累積。當只考慮一個恒星引力影響時的模型如下： .(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、積。如下圖所示當行星偏離理想軌道很小的量以后，之后的偏差就會越來越大， 直至脫離恒星的束縛。在離散化以后，原來臨界穩定的系統變得發散了。二、用高階系統去求解單恒星問題當用高于一階的方法近似求解以上方程時，會取得較好一些的近似。把二階常微分方程組 (1)轉化為一階常微分方程組：3222322200000000()(),xxxyyyxyxxxxyyyy2 / 11 32223222()()xxyyxvxvxyyvyvxy,初始條件是00001001xyxyvv一階常微分方程組00( ,)()xxyfyyy的經典 4 階 rk 法的公式是112341213243()6(,)(,)22(,)22(,)

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

4、微分方程是332222223322222200000000(1)(1)(1)(1)(1)(1), ,xxxxyxyyyyxyxyxx xxyyyy，其中兩個恒星位置是( 1,0)1 0和（，）.用經典 4 階 rk 法求解以上微分方程， 并且在求解過程中根據行星的速度自適應調整迭代的步長 h(變動圍是 0.0005 到 0.005 之間 )。當初值條件為00000,0.275,0.3571,1xyxy時，迭代 5000 步后的軌跡圖如下：4 / 11 在兩個恒星作用下， 初始值選的不好時， 行星在迭代有限次數后會撞到恒星上去。如以上的初始條件在迭代5780 次就會出現行星和恒星的距離小于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 跡圖如下：以上各組測試數據表明， 行星在雙恒星的引力作用下

6、運動軌跡對初始值很敏感。四、參考文獻吳勃英 , 王德明等 .數值分析原理 .:科學.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)+