不使用Numpy的Python电子轨迹计算异常求助

核心错误点 #

• 洛伦兹力计算用错速度变量：循环中计算磁场力时，使用了初始值为[0,0,0]的V，而非当前时刻的速度v[-1]。这导致磁场力始终基于零速度计算，完全偏离实际物理规律，直接引发速度指数增长和轨迹异常。
• 位置更新索引错误：R[j] = r[-1] + V[j]*dt + 0.5*a[j]*dt**2中，r[-1]是完整的位置列表，应该取对应维度的分量r[-1][j]，否则会出现列表与数值相加的错误，导致Z轴轨迹数据异常。
• 变量初始化逻辑混乱：循环内的V和R应该基于当前时刻的速度v[-1]和位置r[-1]计算，而非从零开始赋值后逐步覆盖。

修正后的代码 #

def cross(a, b):#function for cross product
    c = [a[1]*b[2] - a[2]*b[1],
         a[2]*b[0] - a[0]*b[2],
         a[0]*b[1] - a[1]*b[0]]
    return c

# 补充关键物理参数（可根据需求调整）
q = -1.6e-19  # 电子电荷量
m = 9.1e-31   # 电子质量
dt = 1e-5     # 时间步长

E = [1,1,0]
B = [0,0,10]
r = [[0,0,0]]
v = [[0,0,0]]
force = [0,0,0]#list for instantaneous values of force

t = 0
while(t < 1):
    # 获取当前时刻的速度和位置
    current_v = v[-1]
    current_r = r[-1]
    
    a = [0,0,0]
    V = [0,0,0]
    R = [0,0,0]
    
    # 先计算洛伦兹力对应的加速度
    v_cross_B = cross(current_v, B)
    for j in range(3):
        force[j] = q * (v_cross_B[j] + E[j])
        a[j] = force[j] / m
    
    # 用当前加速度更新速度和位置（欧拉法）
    for j in range(3):
        V[j] = current_v[j] + a[j] * dt
        R[j] = current_r[j] + current_v[j] * dt + 0.5 * a[j] * dt**2
    
    v.append(V)
    r.append(R)
    t += dt

额外说明 #

• 补充了q、m、dt等必要物理参数，未定义这些参数会直接导致代码运行失败。
• 将速度和位置的计算拆分：先基于当前速度计算加速度，再用加速度更新速度和位置，逻辑更清晰，避免了循环内变量覆盖导致的错误。
• 如果需要更高精度的轨迹计算，可改用蛙跳法（Leapfrog），这种方法对保守力系统的精度优于欧拉法，更适合模拟洛伦兹力下的运动。

内容的提问来源于stack exchange，提问作者Cengizhan Koyutürk