我特灵解决带电粒子在行星磁场中的运动方程，以查看粒子的路径使用前向欧拉和RK 5方法在python（作为练习学习数值方法）我遇到两个问题：

1. RK 4方法中的“for循环”不更新新值。它为所有迭代给予第一次迭代的值。
   1.随着“β =电荷/质量”符号的改变，粒子的预期路径并没有改变。似乎该路径不受粒子性质（符号）的影响。这在物理上或数学上意味着什么？代码改编自：python two coupled second order ODEs Runge Kutta 4th order和Applying Forward Euler Method to a Three-Box Model System of ODEs我将不胜感激，如果任何人向我解释什么是错误的代码。谢谢。代码如下：

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

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

def LzForce(t,p):
  # assigning each ODE to a vector element
    r,x,θ,y,ϕ,z = p*scales
    
    # constants
    R = 60268e3  # metre
    g_20 = 1583e-9
    Ω = 9.74e-3    # degree/second
    B_θ = (R/r)**4*g_20*cos(θ)*sin(θ)
    B_r = 2*(R/r)**4*g_20*(0.5*(3*cos(θ)**2-1))
    β = +9.36e10

    # defining the ODEs
    drdt = x
    dxdt = r*(y**2 +(z+Ω)**2*sin(θ)**2-β*z*sin(θ)*B_θ)
    dθdt = y
    dydt = (-2*x*y +r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r 
    dϕdt = z
    dzdt = (-2*x*(z+Ω)*sin(θ)-2*r*y*(z+Ω)*cos(θ)+β*(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):
    r0 = 6.6e+07
    x0 = 0. 
    θ0 = 88.
    y0 = 0.
    ϕ0 = 0.
    z0 = 22e-3
    p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])
    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,:]) * dt

    return t, p
   
def rk4(fun,t0,p0,tf,dt):
    # initial conditions
    r0 = 6.6e+07
    x0 = 0. 
    θ0 = 88.
    y0 = 0.
    ϕ0 = 0.
    z0 = 22e-3
    p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])
    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])    
        k2 = dt * fun(t[i] + 0.5*dt, p[i] + 0.5 * k1)
        k3 = dt * fun(t[i] + 0.5*dt, p[i] + 0.5 * k2)
        k4 = dt * fun(t[i] + dt, p[i] + k3)
        p[i+1] = p[i] + (k1 + 2*(k2 + k3) + k4)/6
        return t,p
    
dt = 0.5
tf = 1000
p0 = [6.6e+07,0.0,88.0,0.0,0.0,22e-3]
t0 = 0

 #Solution with Forward Euler
t,p_Euler = ForwardEuler(LzForce,t0,p0,tf,dt)

#Solution with RK4    
t ,p_RK4 = rk4(LzForce,t0, p0 ,tf,dt)

print(t,p_Euler)
print(t,p_RK4)

# Plot Solutions
r,x,θ,y,ϕ,z = p_Euler.T
fig,ax=plt.subplots(2,3,figsize=(8,4))
plt.xlabel('time in sec')
plt.ylabel('parameters')
for a,s in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,s); a.grid()
plt.title("Forward Euler", loc='left')
plt.tight_layout(); plt.show()

r,x,θ,y,ϕ,z = p_RK4.T
fig,ax=plt.subplots(2,3,figsize=(8,4))
plt.xlabel('time in sec')
plt.ylabel('parameters')
for a,q in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,q); a.grid()
plt.title("RK4", loc='left')
plt.tight_layout(); plt.show()

[RK4 solution plot][1]
[Euler's solution methods][2]

''''RK4 does not give iterated values.
The path is unaffected by the change of sign which is expected as it is under Lorentz force''''

  [1]: https://i.stack.imgur.com/bZdIw.png
  [2]: https://i.stack.imgur.com/tuNDp.png