有限差分法计算天体轨道.docx

一几种差分格式的区别1以初速度或者末速度代替平均速度进行向后差分在使用差分法进行迭代计算时，每一次循环称做一轮，循环内部的一次计算称为一步对于天体运动，给定了初始条件之后，就要计算第一轮的速度，v(i + l) = v(i) + ZlVo接 着计算第一轮的位移而这时就会有两种选择，选初始时刻的速度v(i), B|Jx(i + l)=x(i) + 或者选末时刻速度,EPx(i + 1) = %(0 + v(i + l)21to这两种方法的程序结构完全相同取KT的情况(即半径为1的正圆轨道)，时间间隔和总步长相同对于两种方法的结果 对比如下Ar=0 001 N=1E9"方向蓝色线即选初时刻速度，红色线即选末时刻速度所得结果可以看出，在同等时间步长的精度下，红色轨迹己经收敛到一个圆内但是蓝色轨迹仍呈 发散的趋势意即，采用末时刻的速度会加快收敛速度实事上如果把时间间隔取得很小，蓝 色的线也可以收敛到一个比较小的圆当中这个结果是比较好理解的，因为上一步计算得到的速度马上就可以用于下一步位置的计算 如果采用初时刻的速度作为平均速度，那就等于本轮计算得到的速度要到下一轮才被采用所 以釆用末时刻的速度能充分利用每一步的计算结果，提高效率。

所以下面的计算釆用末时刻速度的向后差分进行2中心差分与向后差分的对比下面对比中心差分(蛙跳格式)与向后差分还是之前的问题，在每一轮计算中，无论采 用初时刻的速度还是末时刻的速度作为一个计算步长的平均速度都是明显有问题的，一个可取 的办法便是取中点位置的速度作为平均速度，即蛙跳格式蛙跳格式将时间域和空间域错开 在各自域内其实仍是向后差分，但是整体达到了中心差分的效果在计算上，蛙跳格式只是在 初始条件上多计算了一步，其它的都与向后差分无异1.5-1.5-1.5A t=0.510.50-0.5-1-1-0.50 0.51x方向1.5取一个很夸张的步长匸0・5來突出两者收敛性的不同红色线是蛙跳格式，蓝色线是向后 差分明显看出蛙跳格式的精度更高e方向 方向对于中心差分，虽然步长跨了三个数量级，但是三条线仍基本重合说明当匸0.1时计算 结果就已经比较稳定而向后差分在低精度时已经明显偏离了半径为1的圆看局部放大图•0.C2 -0.010 0.01 0.02 0.03x方向1 0C61.00251 0021.000510.99S60.999x方向x103AT=O.lAr=0.01At=0.00011.0015 叵H 1.001要注意两个图的坐标尺度不同，中心差分的要小很多。

明显可以看出，无论从哪方面对比， 中心差分的精度都要明显高于向后差分付出的代价仅仅是多计算了一步初始条件，但是达到 相同的精度节省了大量的时间和内存消耗因此后面的计算当中都采用蛙跳格式进行二K参数对于行星运行轨道的影响不同K值得到计算得到的行星运行轨道从这里可以看到，当31时，得到的轨迹是一个完整的正圆当KV0.5时，星体逃逸 下面探讨一下其物理过程在有心力场中，考察VXL对时间的变化率d dv dv_(vXL) = -xL=m-X(rXv)GMm(r X v)=GMm一 [r(v - r) - v(r - r)][rrvr — vr2] = —GMmr dr ldri r dt r dt\GMm r3亦即令^(vXL-GMnJ) = OrB = v X £ — GMm — r称为龙格-楞茨矢量(Rung—Lenz vector) [1]要得到天体的运行轨道，只需要计算标积厶2旷・ B =旷• (v XL) — GMmr = £ ・(r X v) — GMmr = GMmrm而厶2厂・ B = rBcosd = GMmrm得到pr = 1 + scosd其中厶2 BP GMm2 ,£ GMm这是一个标准的圆锥曲线方程。

£是偏心率，P是半正焦弦£ = 0时为圆，0V 椭圆，时为抛物线，£>1时为双曲线下面我们來证明K与天体轨道的关系当 K=1 时,即有GM = vxq ,令r = xqo运用两次右手法则，矢量VXL与矢量r是同方向的于是龙格-楞次矢 量的模值可以简化为\B\ = B = X L - GMmf] = mrv2 — GMm = mrv2 — mrv2 = 0于是£ = 0得到方程是一个标准的圆周代入半焦弦长公式，易得p = x0即为半径 K = 时，v2Xq = tv1 = 2GM(1.1)(1.2)(13)(1.4)(1.5)(1.6)(1.7)£ VI时为(1-8)B mrv2 — GMmGMm GMm2GMm — GMmGMm(1-9)1(1.10)对应的轨迹是双曲线 更一般地，可以得到由此总结出下表K£轨迹半焦弦长10圆X0.5~10~1椭圆Xo/K0.51抛物线2x00^0.5>1双曲线xo/K>1<0椭圆Xo/K。