小行星运动轨迹的Runge-Kutta法模拟..doc

小行星运动的Runge-Kutta法模拟 背景介绍 由于两个恒星作用下行星运动问题没有解析解，只能用数值方法求解微分方程。但是在用一阶近似求解微分方程的时候存在严重的误差累积。 当只考虑一个恒星引力影响时的模型如下： ……..(1) 当初始值是时，行星做圆周运动。此时，微分方程的解是。在后面的讨论中，用这个初始条件的方程作为测试方程。 如果采用一阶近似，，就会有严重的误差累积。如下图所示 当行星偏离理想轨道很小的量以后，之后的偏差就会越来越大，直至脱离恒星的束缚。在离散化以后，原来临界稳定的系统变得发散了。 用高阶系统去求解单恒星问题 当用高于一阶的方法近似求解以上方程时，会取得较好一些的近似。 把二阶常微分方程组(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次的运动轨迹图如下： 以上各组测试数据表明，行星在双恒星的引力作用下运动轨迹对初始值很敏感。 参考文献 吴勃英, 王德明等. 数值分析原理. 北京:科学出版社. 2003309-310 Matlab程序1 clear all;clc;close all; ;%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; 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/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/