用數值模擬行星運動，很難得到穩定解？

01-14

我在看《費曼物理學講義（第1卷）》（新千年版），第9章牛頓動力學定律，用彈簧運動和行星運動作為例子，還給出了數值解法。感覺費曼好超前，在那個還沒有PC的年代，就跟學生講用計算機解決問題。我自己用python試了一下，本來跟書上一樣的，後來我想增加計算步數試一試，結果一試發現隨著步數的增加，行星的運動軌跡是一個比一個大的橢圓，相差很明顯。請問這是書中給定的初始條件的原因還是因為這是數值模擬？像真實的行星運動那樣的持續四十多億年的條件是不是很難找？謝謝。
程序字數超了，放鏈接吧： https://github.com/zwdnet/zhihuask/blob/master/9-7.py

首先，感謝大家的對這個答案的支持。

但是很慚愧，我沒有弄清穩態的本質，也就是回答偏了。

題主問題的根本原因，應該是 @浪客所說的「長期模擬行星運動要用專門設計的保能量的辛演算法。常規的演算法會導致能量累積或衰減，高階龍格庫塔或者隱式方法只是誤差累積的慢一些。」

我專門查了一點相關文獻，的確如此。

題主問題中的穩定性，應該是指的「不出現隨著時間一直累積的數值計算誤差」，這個跟微分方程數值求解過程本身的穩定性，是兩個不同的概念。

這裡也非常感謝@ Sun AO，@KL Hahns等知友對我答案的指正。

下面是原答案，請大家隨便看看就好——因為答偏了。

##################### 原答案 #######################

個人也認為主要是問題出在所選用的數值演算法上，比如，時間步長不夠小。

說來也奇怪，我們國家的數值分析類的課程，對微分方程的隱式數值演算法介紹的相當不夠，相關內容老舊，N年不換，導致遇到實際應用的時候，例如，比較剛性的微分方程(組)，經常很難或者無法成功數值求解。

隱式演算法對數值求解穩定性方面很有優勢，也可以控制求解精度，同時，可以進行變步長計算(在要求的精度下)，以節約求解時間。

如上面的知友所述，可以考慮隱式歐拉法。隱式歐拉法是最簡單的隱式常微分方程(組)演算法，但需要利用牛頓法等迭代求解非線性方程(組)。題主有興趣可以搜搜相關文獻和程序。

數值計算方法的重要性，其實遠超大部分科學和工程技術人員的想像。

##################### 原答案 #######################

我在這裡貼一個我之前遇到的相關聯的問題。。

為什麼這個演算法誤差的看起來這麼小？ - 知乎

簡而言之，確實需要專門做理論推導才能設計出來合適的差分方程來進行計算機模擬。

你的計算精度太低了。誤差積少成多。把delta取小一些會更好些。拿你的代碼，我取delta = 0.00001，然後run（1000000）,大約運行了十幾圈，沒有問題。但你能看到軌道遠日點處比較粗，這就是已經發生的誤差。這種誤差都是先從遠日點開始的。所以說單純減小步長圈數多了仍然會有問題。我記得我剛學編程的時候寫的模擬太陽系差不多也就是轉100年的樣子，之後好像是火星軌道就開始發生明顯可見的問題了。

不過其實減小步長的辦法是最直接的解決辦法。從原理上講，它可以將誤差無限減小。但這個辦法的缺點很明顯。計算量增加得太大了。有很多減小和修正誤差的辦法。實際處理中會平衡運算量和誤差，比如評論里提到的龍格庫塔和蛙跳演算法。

得到高穩定性且高精度的系統大約至少需要先修個數值計算。高精度（且相對低計算量）的行星運動模擬應該是專門的一門學問，不要小看這門科目。至少攝動什麼的肯定是要學的。

——————————————————————————————————————

更新：剛看了一下其他答主的答案。他們比我更專業，我這個答案就算拋磚引玉吧。

兩體問題一般是二次曲線，在遠離拋物線、圓、直線這些臨界情形時不涉及初始狀態。是橢圓就是橢圓，是雙曲線就是雙曲線。

你用的是顯式歐拉折線法。

隱式歐拉折線法無條件穩定。

顯式歐拉折線法條件穩定，且需要步長充分小。

不過這個穩定性指的是小的初始誤差帶來的誤差的，而你口中的穩定性似乎只是舍入誤差和截斷誤差。

不過，這一步的舍入誤差和截斷誤差放到下一步就是初始的誤差，所以套用穩定性的分析也沒錯。

改成極坐標，並用龍格庫塔四階方法去計算，會好些。

推薦你一本書

Computational Physics

by Mark Newman

google可看前四章

可能離題了。樓主問的是數值方法，目的是減少誤差，但是實際上算軌道力學更簡便的近似法；二體問題就用開普勒軌道元素算就能得到非常好的接近解。天然衛星和人造衛星有各自的多餘因素需要考慮。

https://en.m.wikipedia.org/wiki/Orbital_elements

開普勒元素對於行星來說很長時間都是比較穩定的，以下這篇兩頁的論文一步步教你從軌道元素算出位置和速度向量

https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf

考慮中心天體的形狀不是完全的球體，要把J2項帶入計算從而得到升交點和近日（地）點的進動。太陽同步衛星就是這麼來的。

對於偏心率大的天體，要考慮相對論效應而產生的近日點進動。例如水星。

對於近地軌道上的衛星，要考慮大氣的drag。

我記得解微分方程，第一講就是如果搞不定就么上計算機。

直接模擬有系統誤差，比如你就算模擬一個二次函數，如果直接按定義來模擬，由於步長不是0，一定會偏向外側/內側。

估計用容格庫塔會好很多，一味減小步長不是正路。

不能忍受忽略了所有物理常數的物理公式！！

應該是這樣：

import numpy as np import matplotlib.pyplot as plt


def run(N):
    t = np.arange(0.0, N, 0.1)
    x = np.zeros(N*10)

    vx = np.zeros(N*10)

    ax = np.zeros(N*10)
    y = np.zeros(N*10)

    vy = np.zeros(N*10)

    ay = np.zeros(N*10)
    r = np.zeros(N*10)
    G = 1

    m = 1

    d = 1
    def Acceleration(d, G, m):

        return -G*m/(d**2)
    x[0] = 0.5*d

    vx[0] = -0.2*d

    y[0] = 0.0*d

    vy[0] = 1.63*d

    r[0] = np.sqrt(x[0]**2 + y[0]**2)

    ax[0] = Acceleration(r[0], G, m)*x[0]/r[0]

    ay[0] = Acceleration(r[0], G, m)*y[0]/r[0]
    delta = 0.01/(np.sqrt(ax[0]**2 + ay[0]**2))

    for i in range(1, N*10):

        vx[i] = vx[i-1] + ax[i-1]*delta

        vy[i] = vy[i-1] + ay[i-1]*delta

        x[i] = x[i-1] + vx[i-1]*delta

        y[i] = y[i-1] + vy[i-1]*delta

        r[i] = np.sqrt(x[i]**2 + y[i]**2)

        ax[i] = Acceleration(r[i], G, m)*x[i]/r[i]

        ay[i] = Acceleration(r[i], G, m)*y[i]/r[i]
    plt.title("Simulation Times: %d" % int(N*10))

    plt.plot(x,y,"o")

    plt.show()

if __name__ == "__main__": run(2500)

作個圖發現軌跡跑遠了：

說明步長設定為0.01太大了。改小點試試！

delta = 0.001/(np.sqrt(ax[0]**2 + ay[0]**2))

果然好多了：

收工。

我建議你試一下 Verlet 演算法。The Verlet algorithm

看你想模擬什麼問題了。混沌的系統理論上是不可模擬的。當然你嘗試的例子也可能能夠穩定模擬，但是這是特例不是普遍規律。