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

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

得到这个：

为什么轨迹的形状不同？

Answer 1

是的，你是对的。初始条件似乎是问题所在。您在 dy4 和 dy5 中除以零，这导致了 RuntimeWarning: divide by zero encountered.

将起始条件替换为：

y0 = [10, 20., 20., 1., 1.e3, 0]

给出以下输出：

Answer 2

中心力在半径方向确实有F=-c/r^2、c=e**2/(4*pi*eps_0*me)的大小。要获得矢量值力，需要将其与远离中心的方向矢量相乘。这给出了 F=-c*x/r^3 的向量力，其中 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 s 之后它将因此有 1e8 m 的距离，这段代码的结果证实了这一点。

关于添加的双曲线图swing-by，请注意开普勒定律的圆锥截面参数化为

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

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

How to set initial condition to get the trajectory as a solution of the equation of the motion?

python

numpy

scipy

ode

odeint