天體運動的簡單數值計算

03-12

閱讀原文

預備知識　萬有引力，彈簧振子受迫運動的簡單數值計算

　　直角坐標系中，設中心天體質量為 $M$ ，固定在原點不動．根據牛頓萬有引力定律，質量為 $m$ 的行星受到中心天體的力為

$mathbf F = -Gfrac{Mm}{r^2}hat {mathbf r} = - G frac{Mm}{r^3} mathbf r$ 　　(1)

其中 $mathbf r$ 為行星的位矢（設行星在 xy 平面上運動）．根據牛頓第二定律，加速度為

$mathbf a = frac{mathbf F}{m} = -G frac M{r^3}mathbf r$ 　　(2)

用 $r = sqrt{x^2 + y^2}$ 以及 $mathbf r = x hat {mathbf x} + yhat {mathbf y}$ ，其中 $x,y$ 看成 $t$ 的函數．考慮到 $a_x = mathrm d^2 x/mathrm d t^2$ ， $a_y = mathrm d^2 y/mathrm d t^2$ 可以列出二階微分方程組

$left{egin{aligned} &frac{mathrm d^2 x}{mathrm d t^2} = - frac{GM}{(x^2 + y^2)^{3/2}} x\ &frac{mathrm d^2 y}{mathrm d t^2} = - frac{GM}{(x^2 + y^2)^{3/2}} y end{aligned} ight.$ 　　(3)

假設已知初值條件 $x(0) = x_0, , y(0) = y_0, , dot x(0) = v_{x0}, , dot y(0) = v_{y0}$ ．下面用「簧振子受迫振動的簡單數值計算」中類似的方法求接下來行星的運動軌跡．

（剩下部分見頂部的「閱讀原文」）