scipy.odeint返回二階非線性微分方程不正確的值

我一直在試圖解決牛頓萬有引力定律（平方反比定律）二階非線性微分方程：scipy.odeint返回二階非線性微分方程不正確的值

x(t)'' = -GM/(x**2)

的衛星接近地球的在一名維

使用numpy.odeint用一系列第一ORD的運動（在這種情況下，點質量） er微分方程，但與Mathematica或簡化形式的法則（Δx=（1/2）在^ 2）相比，該操作產生了不正確的結果。

這是程序中的代碼：

import numpy as np 
from scipy.integrate import odeint 

def deriv(x, t):   #derivative function, where x[0] is x, x[1] is x' or v, and x2 = x'' or a 
    x2 = -mu/(x[0]**2) 
    return x[1], x2 

init = 6371000, 0   #initial values for x and x' 

a_t = np.linspace(0, 20, 100) #time scale 

mu = 398600000000000  #gravitational constant 

x, _ = odeint(deriv, init, a_t).T 

sol = np.column_stack([a_t, x])

其產生與耦合A_T和x位置值作爲衛星的陣列從6371000米（地球的平均半徑的初始距離接近地球）。例如，人們會期望物體從地面6371000m降到6370000m需要大約10秒鐘的時間（因爲1000 = 1/2（9.8）（t^2）的解幾乎是10）並且微分方程的數學解決方案也將該值稍微高於10s。

然而根據odeint溶液和溶膠陣列值是幾乎14.4。

t     x 
    [ 1.41414141e+01, 6.37001801e+06], 
    [ 1.43434343e+01, 6.36998975e+06],

有沒有在odeint解決方案顯著的錯誤，或者是有一個主要的問題在我的函數/ odeint使用？謝謝！

來源

2015-08-24 JAustin

一切看起來基本上是正確的。我會試着用更小的時間步驟來看看你會得到什麼。 –

（因爲該解決方案，以1000 = 1/2（9.8）（T^2）幾乎正好10），

這是一個正確的完整性檢查，但有事情和算術的了。使用這種近似，我們得到一個t的

>>> (1000/(1/2 * 9.8))**0.5 
14.285714285714285

截至〜10反對，這會給我們只有

>>> 1/2 * 9.8 * 10**2 
490.00000000000006

的這種期望〜14.29的距離非常接近的結果你觀察到：

>>> sol[abs((sol[:,1] - sol[0,1]) - -1000).argmin()] 
array([ 1.42705427e+01, 6.37000001e+06])

來源

2015-08-24 21:35:33 DSM

非常感謝。我現在覺得自己像個白癡，因爲我花了好幾個小時試圖弄清楚Python爲什麼不同意我的錯誤算法。我必須弄清楚爲什麼Wolfram情節顯然不正確，但這可能是我的錯誤，因爲我對Mathematica缺乏經驗。 – JAustin

scipy.odeint返回二階非線性微分方程不正確的值

回答

相關問題