Python求解洛伦兹加速度运动方程：RK4与β符号问题排查

核心问题 #

• RK4方法的迭代未正常执行，所有时间步输出的都是第一次迭代的结果
• 修改β（电荷/质量比）的符号后，粒子轨迹未发生预期变化，似乎不受电荷符号影响

代码错误解析 #

1. RK4函数循环提前返回 #

rk4函数的for循环中，return t,p语句被放在循环内部，导致循环仅执行一次就直接返回结果，后续迭代步骤完全未运行，这是所有输出都是第一次迭代结果的直接原因。

2. β被硬编码导致符号修改无效 #

LzForce函数内部硬定义了β = +9.36e10，无论外部如何修改β的值，函数都会使用内部固定值，因此改变电荷符号的操作根本没有传递到运动方程中，自然不会影响粒子轨迹。

3. 额外问题：初始条件被覆盖 #

ForwardEuler和rk4函数都在内部重新定义了p0，覆盖了用户传入的初始条件参数，导致传入的初始值失效，不符合函数设计的预期。

物理层面说明 #

正常情况下，带电粒子在磁场中受洛伦兹力作用，电荷符号改变会导致洛伦兹力方向反转，粒子轨迹应呈现镜像变化。当前轨迹未变化完全是代码错误导致的，并非物理规律异常。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos

scales = np.array([1e7, 0.1, 1, 1e-5, 10, 1e-5])

def LzForce(t,p, beta):
    # 解包缩放后的状态量
    r,x,θ,y,ϕ,z = p*scales
    
    # 物理常数
    R = 60268e3  # 行星半径（米）
    g_20 = 1583e-9
    Ω = 9.74e-3    # 自转角速度（度/秒）
    B_θ = (R/r)**4*g_20*cos(θ)*sin(θ)
    B_r = 2*(R/r)**4*g_20*(0.5*(3*cos(θ)**2-1))

    # 定义运动微分方程
    drdt = x
    dxdt = r*(y**2 +(z+Ω)**2*sin(θ)**2 - beta*z*sin(θ)*B_θ)
    dθdt = y
    dydt = (-2*x*y + r*(z+Ω)**2*sin(θ)*cos(θ) + beta*r*z*sin(θ)*B_r)/r 
    dϕdt = z
    dzdt = (-2*x*(z+Ω)*sin(θ) - 2*r*y*(z+Ω)*cos(θ) + beta*(x*B_θ - r*y*B_r))/(r*sin(θ))

    return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])/scales

def ForwardEuler(fun,t0,p0,tf,dt, beta):
    t = np.arange(t0,tf+dt,dt)
    p = np.zeros([len(t), len(p0)])
    p[0] = p0

    for i in range(len(t)-1):
        p[i+1,:] = p[i,:] + fun(t[i],p[i,:], beta) * dt

    return t, p
   
def rk4(fun,t0,p0,tf,dt, beta):
    t = np.arange(t0,tf+dt,dt)
    p = np.zeros([len(t), len(p0)])
    p[0] = p0
     
    for i in range(len(t)-1):
        k1 = dt * fun(t[i], p[i], beta)    
        k2 = dt * fun(t[i] + 0.5*dt, p[i] + 0.5 * k1, beta)
        k3 = dt * fun(t[i] + 0.5*dt, p[i] + 0.5 * k2, beta)
        k4 = dt * fun(t[i] + dt, p[i] + k3, beta)
        p[i+1] = p[i] + (k1 + 2*(k2 + k3) + k4)/6
    # 将return移到循环外部
    return t,p
    
# 参数设置
dt = 0.5
tf = 1000
p0 = np.array([6.6e+07, 0.0, 88.0, 0.0, 0.0, 22e-3])
t0 = 0
beta_pos = +9.36e10  # 正电荷
beta_neg = -9.36e10  # 负电荷

# 用Forward Euler求解正负电荷情况
t,p_Euler_pos = ForwardEuler(LzForce,t0,p0,tf,dt, beta_pos)
t,p_Euler_neg = ForwardEuler(LzForce,t0,p0,tf,dt, beta_neg)

# 用RK4求解正负电荷情况
t ,p_RK4_pos = rk4(LzForce,t0, p0 ,tf,dt, beta_pos)
t ,p_RK4_neg = rk4(LzForce,t0, p0 ,tf,dt, beta_neg)

# 绘制Forward Euler结果（正负电荷对比）
fig,ax=plt.subplots(2,3,figsize=(12,6))
plt.suptitle("Forward Euler 正负电荷轨迹对比", y=1.02)
for idx, (a,s_pos,s_neg) in enumerate(zip(ax.flatten(), p_Euler_pos.T, p_Euler_neg.T)):
    labels = ["r","x","θ","y","ϕ","z"]
    a.plot(t,s_pos, label=f"β={beta_pos}")
    a.plot(t,s_neg, label=f"β={beta_neg}", linestyle='--')
    a.set_xlabel('时间（秒）')
    a.set_ylabel(labels[idx])
    a.grid()
    a.legend()
plt.tight_layout(); plt.show()

# 绘制RK4结果（正负电荷对比）
fig,ax=plt.subplots(2,3,figsize=(12,6))
plt.suptitle("RK4 正负电荷轨迹对比", y=1.02)
for idx, (a,s_pos,s_neg) in enumerate(zip(ax.flatten(), p_RK4_pos.T, p_RK4_neg.T)):
    labels = ["r","x","θ","y","ϕ","z"]
    a.plot(t,s_pos, label=f"β={beta_pos}")
    a.plot(t,s_neg, label=f"β={beta_neg}", linestyle='--')
    a.set_xlabel('时间（秒）')
    a.set_ylabel(labels[idx])
    a.grid()
    a.legend()
plt.tight_layout(); plt.show()

修正后的效果 #

• RK4函数的循环会完整执行所有时间步，生成正确的迭代轨迹
• 修改β符号后，正负电荷的轨迹会呈现预期的差异，符合洛伦兹力的物理规律

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