I am tring to solve the equation of motion of charged particle in planetary magnetic field to see the path of the particle using Forward Euler's and RK5 method in python (as an excercise in learning Numerical methods) I encounter two problems:

1. The 'for loop' in the RK4 method does not update the new values. It give the values of the first iteration for all iteration.
2. With the change of the sing of 'β = charge/mass' the path of particle which is expected does not change. It seems the path is unaffected by the nature(sign) of the particle. What does this mean physically or mathematically? The codes are adapted from : python two coupled second order ODEs Runge Kutta 4th order and Applying Forward Euler Method to a Three-Box Model System of ODEs I would be immensely grateful if anyone explain to me what is wrong in the code. thank you. The Code are as under:

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

CodePudding user response：

You are not iterating more than once inside the for loop in rk4 because it returns after the first iteration.

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 
# This is the problem line, the return was tabbed in, to be inside the for block, so the block executed once and returned.
return t,p

For physics questions please try a different forum.