scipy.odeint 返回二阶非线性微分方程的错误值

Question

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

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

对于卫星在一维空间中接近地球（在本例中为质点）的运动

对一系列一阶微分方程使用 numpy.odeint，但与 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 位置值的数组。例如，人们会期望物体从 6371000m 下降到 6370000m 的表面 1000 米大约需要 10 秒（因为 1000 = 1/2(9.8)(t^2) 的解几乎正好是 10），微分方程的 Mathematica 解也使该值略高于 10s。

然而，根据 odeint 解和 sol 数组，该值接近 14.4。

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

是odeint解决方案有明显错误，还是我function/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

相对于 ~10 的 t，这只会给我们提供

的距离

>>> 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])

scipy.odeint 返回二阶非线性微分方程的错误值

scipy.odeint returning incorrect values for second order non-linear differential equation

python

numpy

matplotlib

scipy

differential-equations