python - 用 sympy 简化 ODE，然后用 Scipy 进行数值求解

Question

我想用数值求解 E/M 场中相对论粒子的运动方程：

其中, , 和,
。

我已经开始使用 Sympy 进行简化来求解. 使用

""" Define symbols"""
t,m,q, c = smp.symbols('t m q c', Real =True)
x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz, gamma = smp.symbols('x y z vx vy vz  E_x E_y E_z B_x B_y B_z gamma ',cls=smp.Function, real=True)

""" Define functions"""

x    = x(t);                y   = y(t);               z    = z(t);                 # Position   
vx   = vx(t);              vy   = vy(t);              vz   = vz(t);                # Velocity  
Ex   = Ex(x, y, z, t);     Ey   = Ey(x, y, z, t);     Ez   = Ez(x, y, z, t);       # Components of E field 
Bx   = Bx(x, y, z, t);     By   = By(x, y, z, t);     Bz   = Bz(x, y, z, t);       # Components of B field
x_d  = smp.diff(x,t);      y_d  = smp.diff(y,t);      z_d  = smp.diff(z,t);        # dx_{i}/dt ,    x_{i} = x, y, z
vx   = x_d;                vy   = y_d;                vz   = z_d;                  # v_{i} ,           i  = x, y, z 
vx_d = smp.diff(x,t,2);    vy_d = smp.diff(y,t,2);    vz_d = smp.diff(z,t,2);      # d^{2}x_{i}/dt^2 , i  = x, y, z 

""" Define Matrices"""
r      = smp.Matrix([x, y, z])
E      = smp.Matrix([Ex, Ey, Ez])
B      = smp.Matrix([Bx, By, Bz])
v      = smp.Matrix([vx, vy, vz])

""" Define gamma """
v_norm = v.norm()
gamma  = gamma(v_norm)
gamma  = 1/smp.sqrt(1 - (v.norm()/c)**2)

""" Define momentum """
p    = m*gamma*v

""" 2nd order ODE of motion """
de1 = smp.diff(p,t)-q*(E +v.cross(B)/c)

""" Solve for dv_{i}/dt, i = x, y, z """ 
sols   = smp.solve([de1[0],de1[1], de1[2]], (vx_d, vy_d, vz_d), Simplify =False, rational=False)

"""Create numpy functions that we can use with numerical methods"""
list_p = (c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz)


dvx_dt_f  = smp.lambdify(list_p,sols[vx_d], modules=smp)
dx_dt_f   = smp.lambdify(vx,vx, modules=smp)

dvy_dt_f  = smp.lambdify(list_p,sols[vy_d], modules=smp)
dy_dt_f   = smp.lambdify(vy, vy, modules=smp)

dvz_dt_f  = smp.lambdify(list_p,sols[vz_d], modules=smp)
dz_dt_f   = smp.lambdify(vz, vz, modules=smp)

然后我想使用 lamdified 表达式并使用 Scipy 的 ivp 以以下形式集成它们。

def dSdt(S,t):
    vx, x, vy, y, vz, z = S
    return[dvx_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dx_dt_f(vx),
           dvy_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dy_dt_f(vy),
           dvz_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dz_dt_f(vz)
          ]

对于初始条件，v0 = (x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz)

然而，我没有根据数值积分所需的 Vx、Vy、Vz 来获得 Sympy 解，而是根据时间导数 $\frac{dx}{dt} 来获得它。

我可以做些什么来改变 $v_x、v_y、v_z$ 的导数？

标签： pythonscipysympynumerical-integrationnumerical-analysis