计算科学 - 为什么是1 /r21/r2力定律给出螺旋轨迹？ - 吾爱随笔录

我已经编写了一个程序来解决牛顿第二运动定律在二维极坐标中给定的力定律。众所周知，如果力定律的形式为 $k/r^2$ ，我们得到圆锥截面作为解决方案（https://en.wikipedia.org/wiki/Kepler_problem）。

但是当我运行程序时，我得到一个螺旋轨迹（它不是圆锥曲线之一）。

分量方向，牛顿第二定律的中心力 $F(r)=k/r^2$ （和单位质量）是-

\frac{d^{2} r}{d t^{2}} = k / r^{2}

$\frac{d^2r}{dt^2}=k/r^2$

\frac{d^{2} θ}{d t^{2}} = 0

$\frac{d^2\theta}{dt^2}=0$ 其中第一个方程是径向部分，第二个方程是角部分（由于角动量守恒，中心力问题消失）。

程序中要使用的方程是一组 4 个耦合的一阶 ODE，从上述二阶 ODE 导出。（因为 odeint 不能直接处理二阶 ODE）。我在程序中使用的方程是——

d r / d t = v_{r}

$dr/dt=v_r$

d v_{r} / d t = k / r^{2}

$dv_r/dt=k/r^2$

d θ / d t = v_{θ}

$d\theta/dt=v_\theta$

d v_{θ} / d t = 0

$dv_\theta/dt=0$ 和

k = 1

$k=1$ 为简单起见。

我尝试了许多可能的初始条件，并尝试运行不同的值 $k$ 确保。为什么会这样？这是因为初始条件的选择还是一些错误的实现或一些数值问题？还是他们的物理问题？

下面是我的最小代码

def vec(w, t):
    r, vr, theta, vtheta = w
    return [vr, r**(-2.0), vtheta, 0]

def newton(vec, initial, t):
    wsol = odeint(vec, initial, t)
    return [wsol[:, 0], wsol[:, 2]]

T = np.linspace(0, 50, 1000)    
initial = [2, 10, radians(0),2] #The order is radius, radial velocity, angle, angular velocity
R = newton(vec, initial, T)[0]
Theta = newton(vec, initial, T)[1]
plt.polar(Theta, R, "r", lw="1")
plt.show()