python - 如何设置初始条件以获得轨迹作为运动方程的解？

Question

我想看看运动方程与静电力的解。下面的脚本有什么问题？是初始状态的问题吗？谢谢

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def dr_dt(y, t):
    """Integration of the governing vector differential equation.
    d2r_dt2 = -(e**2/(4*pi*eps_0*me*r**2)) with d2r_dt2 and r as vecotrs.
    Initial position and velocity are given.
    y[0:2] = position components
    y[3:] = velocity components"""

    e = 1.602e-19 
    me = 9.1e-31
    eps_0 = 8.8541878128e-12
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2)

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -e**2/(4 * np.pi * eps_0 * me * (y[0])**2)
    dy4 = -e**2/(4 * np.pi * eps_0 * me * (y[1])**2)
    dy5 = -e**2/(4 * np.pi * eps_0 * me * (y[2])**2)
    return [dy0, dy1, dy2, dy3, dy4, dy5]

t = np.arange(0, 100000, 0.1)
y0 = [10, 0., 0., 0., 1.e3, 0.]
y = odeint(dr_dt, y0, t)
plt.plot(y[:,0], y[:,1])
plt.show()

这是轨迹形状的预期结果：

我应用以下初始条件：

t = np.arange(0, 2000, 0.1)   
y01 = [-200, 400., 0., 1., 0, 0.]                   
y02 = [-200, -400., 0., 1., 0, 0.]

得到了这个：

为什么轨迹的形状不同？

标签： pythonnumpyscipyodeodeint

Answer 1

中心力在半径方向上F=-c/r^2的大小确实为 。c=e**2/(4*pi*eps_0*me)要获得向量值的力，需要将其与远离中心的方向向量相乘。F=-c*x/r^3这给出了where的向量力r=|x|。

def dr_dt(y, t):
    """Integration of the governing vector differential equation.
    d2x_dt2 = -(e**2*x/(4*pi*eps_0*me*r**3)) with d2x_dt2 and x as vectors, 
    r is the euclidean length of x.
    Initial position and velocity are given.
    y[:3] = position components
    y[3:] = velocity components"""

    e = 1.602e-19 
    me = 9.1e-31
    eps_0 = 8.8541878128e-12
    x, v = y[:3], y[3:]
    r = sum(x**2)**0.5
    a = -e**2*x/(4 * np.pi * eps_0 * me * r**3)
    return [*v, *a]

t = np.arange(0, 50, 0.1)
y0 = [-120, 40., 0., 4., 0., 0.]
y = odeint(dr_dt, y0, t)
plt.plot(y[:,0], y[:,1])
plt.show()

初始条件[10, 0., 0., 0., 1.e3, 0.]对应于与固定在原点的质子相互作用的电子，从 10m 开始并以 1000 m/s 垂直于半径方向移动。在旧物理学中（您需要一个适当的微米级粒子过滤器），这直观地意味着电子几乎不受阻碍地飞走，因此在 1e5 秒后它将有 1e8 m 的距离，此代码的结果证实了这一点。

---

至于添加的双曲线摆动图，请注意开普勒定律的圆锥截面参数化为

r(φ)= R/(1+E*cos(φ-φ_0))

其最小半径r=R/(1+E)为φ=φ_0。所需的渐近线大约在顺时针方向移动φ=π。φ=-π/4然后给出φ_0=3π/8对称轴的角度和E=-1/cos(π-φ_0)=1/cos(φ_0)偏心度。要将水平渐近线的高度纳入计算，请将y坐标计算为

y(φ) = R/E * sin(φ)/(cos(φ-φ_0)+cos(φ0))
     = R/E * sin(φ/2)/cos(φ/2-φ_0)

in the limit

y(π)=R/(E*sin(φ_0))

随着图形集y(π)=b，我们得到R=E*sin(φ_0)*b。

远离原点的速度为r'(-oo)=sqrt(E^2-1)*c/K，其中K为常数或面积定律r(t)^2*φ'(t)=K，与K^2=R*c。取其为初点水平方向的近似初速度[ x0,y0] = [-3*b,b]。

因此，b = 200这导致了初始条件向量

y0 = [-600.0, 200.0, 0.0, 1.749183465630342, 0.0, 0.0]

并整合T=4*b给情节

上面代码中的初始条件是b=40生成相似图像。

更多的镜头使用b=200，将初始速度乘以0.4,0.5,...,1.0,1.1