西皮·奥德林二阶非线性微分方程的返回错误值

Question

我一直在试图解牛顿万有引力定律（平方反比定律）的二阶非线性微分方程：

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

卫星在一维接近地球的运动（在这里是点质量）

使用纽比.奥迪恩对于一系列一阶微分方程，但与Mathematica或简化形式的定律（∆x=（1/2）at^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])

当卫星从6371000米的初始距离（地球的平均半径）接近地球时，产生一个阵列，其中a_t和x位置值耦合。例如，人们会认为，物体从6371000米下降到637万米，需要大约10秒的时间（因为1000=1/2（9.8）（t^2）的解几乎正好是10秒），而微分方程的mathematica解也使该值略高于10秒。在

然而，根据odeint解和sol阵列，这个值接近14.4。在

^{pr2}$

odeint解决方案中是否存在重大错误，或者我的函数/odeint使用中是否存在重大问题？谢谢！在

Answer 1

(because the solution to 1000 = 1/2(9.8)(t^2) is almost exactly 10),

这是正确的检查，但你的算术出了问题。使用这个近似，我们得到t

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

与t为~10相反，它只给我们一个距离

^{pr2}$

这个~14.29的预期值非常接近您观察到的结果：

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