（建議閱讀全文）

預備知識　萬有引力， 彈簧振子受迫運動的簡單數值計算 　　 下面我們來用一種極其簡單的算法對單個天體在中心天體的萬有引力作用下的運動進行數值計算． 事實上該問題存在解析解（見開普勒三定律）， 所以以下的算法只是用于演示數值解常微分方程的大致原理． 這種方法可以輕易地拓展到多個天體的情況， 而多體情況沒有一般的解析解． 一種效率更高的常見常微分方程數值算法可參考 “四階龍格庫塔法”． 　　 直角坐標系中， 設中心天體質量為

， 固定在原點不動．根據牛頓萬有引力定律，質量為 的行星受到中心天體的力為

其中

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

用

以及 ， 其中 看成 的函數． 考慮到 ， ， 可以列出二階微分方程組（變量上方兩點表示關于時間的二階導數）

假設已知初值條件

， ， ， ． 下面用 “彈簧振子受迫振動的簡單數值計算” 中類似的方法求接下來行星的運動軌跡． 　　 1．將初始條件代入式 3 ，得到初始加速度

　　

2．設經過一段極微小的時間步長

（例如 ， 數值越小誤差越小）， 根據微分近似（微分近似在這里的物理意義是在 內速度和加速度都近似為常數）

　　

3．把

再次代入式 3 ， 得到 ， 再次利用微分近似求出 如此循環下去就可以得到每隔 的數值解． 代碼 1：kepler.m % 參數設定 GM = 1; % 萬有引力常數乘以中心天體質量 x0 = 1; y0 = 0; % 初始位置 vx0 = 0; vy0 = 0.7; % 初始速度 T = 4; Nstep = 4000; % 總時間和步數 dt = T/Nstep; % 步長% 矩陣預賦值 x = nan(Nstep,1); y = x; x1 = x; y1 = x; x2 = x; y2 = x;% 初始位置，速度，加速度 x(1) = x0; y(1) = y0; % 初位置 x1(1) = vx0; y1(1) = vy0; % 初速度 x2(1) = -GM*x(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 x''(0) y2(1) = -GM*y(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 y''(0)% 迭代循環 for ii = 2:Nstepx(ii) = x(ii-1)+x1(ii-1)*dt; % x的微分y(ii) = y(ii-1)+y1(ii-1)*dt; % y的微分x1(ii) = x1(ii-1)+x2(ii-1)*dt; % x' 的微分y1(ii) = y1(ii-1)+y2(ii-1)*dt; % y' 的微分x2(ii) = -GM*x(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 x''y2(ii) = -GM*y(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 y'' end% 畫圖 plot(x,y); % 畫行星軌道 axis equal; % xy坐標長度一致 hold on; % 繼續畫圖 scatter(0,0); % 標出中心天體程序運行結果如圖 1 所示． 注意行星軌道并不是一個閉合的橢圓， 這是由于這種算法誤差較大， 為了減小誤差， 可以增加程序中 Nstep 的值．

圖 1：運行結果

總結

以上是生活随笔為你收集整理的pstate0 vid数值意义_天体运动的简单数值计算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。