Sympy简化相对论粒子运动ODE后Scipy数值求解的符号替换问题

问题根源 #

你当前的代码中将vx直接赋值为smp.diff(x,t)，这就导致后续所有用到vx的位置都会被自动替换为x的一阶导数符号，因此最终求解得到的加速度表达式中自然不会保留vx、vy、vz作为独立变量，而是直接保留导数符号。

修正步骤 #

1. 调整符号定义，解绑速度和位置导数的绑定 #

不要将位置、速度、场分量定义为时间的函数，直接定义为普通实值符号，仅保留加速度为待求解的符号：

import sympy as smp
import numpy as np

# 定义常数符号
t, m, q, c = smp.symbols('t m q c', real=True)
# 定义独立输入变量：位置、速度、电磁场分量
x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz = smp.symbols('x y z vx vy vz E_x E_y E_z B_x B_y B_z', real=True)
# 定义待求解的加速度分量符号
ax, ay, az = smp.symbols('a_x a_y a_z', real=True)

2. 用独立速度符号定义动量和运动方程 #

直接用vx,vy,vz计算洛伦兹因子和动量，求时间导数时手动应用链式法则替换：

v = smp.Matrix([vx, vy, vz])
a = smp.Matrix([ax, ay, az])

# 定义洛伦兹因子
gamma = 1 / smp.sqrt(1 - (v.norm()/c)**2)
# 动量对时间求导，应用链式法则：dp/dt = m*gamma*a + m*gamma**3*(v·a)/c² * v
dp_dt = m * gamma * a + m * gamma**3 * (v.dot(a))/c**2 * v
# 洛伦兹力
F = q * (smp.Matrix([Ex, Ey, Ez]) + v.cross(smp.Matrix([Bx, By, Bz])))

# 运动方程：dp_dt = F
eqs = [dp_dt[i] - F[i] for i in range(3)]

3. 求解加速度分量并生成数值函数 #

求解得到的ax,ay,az表达式仅包含你定义的普通符号，没有导数项，直接生成lambdify函数即可：

# 求解加速度分量
sols = smp.solve(eqs, (ax, ay, az), simplify=False, rational=False)

# 生成数值函数
input_params = (c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz)
dvx_dt_f = smp.lambdify(input_params, sols[ax], modules='numpy')
dvy_dt_f = smp.lambdify(input_params, sols[ay], modules='numpy')
dvz_dt_f = smp.lambdify(input_params, sols[az], modules='numpy')

4. 调整数值积分函数 #

注意如果电磁场是随位置/时间变化的，需要在dSdt中先根据当前的位置、时间计算出瞬时场值，再传入加速度函数：

from scipy.integrate import odeint

# 示例：假设是恒定均匀场，实际使用时替换为你的场函数
E_field = (1e5, 0, 0)
B_field = (0, 0, 1)
# 物理常数
c_val = 3e8
m_val = 9.1e-31
q_val = 1.6e-19

def dSdt(S, t):
    vx, x, vy, y, vz, z = S
    # 计算当前时刻的电磁场值，这里用恒定场示例
    Ex_val, Ey_val, Ez_val = E_field
    Bx_val, By_val, Bz_val = B_field
    # 计算加速度
    ax = dvx_dt_f(c_val, m_val, q_val, x, y, z, vx, vy, vz, Ex_val, Ey_val, Ez_val, Bx_val, By_val, Bz_val)
    ay = dvy_dt_f(c_val, m_val, q_val, x, y, z, vx, vy, vz, Ex_val, Ey_val, Ez_val, Bx_val, By_val, Bz_val)
    az = dvz_dt_f(c_val, m_val, q_val, x, y, z, vx, vy, vz, Ex_val, Ey_val, Ez_val, Bx_val, By_val, Bz_val)
    return [ax, vx, ay, vy, az, vz]

# 示例初值：[vx0, x0, vy0, y0, vz0, z0]
S0 = [0, 0, 1e6, 0, 0, 0]
t_arr = np.linspace(0, 1e-8, 1000)
sol = odeint(dSdt, S0, t_arr)

内容的提问来源于stack exchange，提问作者jokerp