# Euler's Method for fast moving particle trajectory

I'm trying to figure out how a magnet affects the trajectory of a particle travelling near the speed of light downwards toward the ground. The equation for the force of the magnetic field is pretty difficult to solve, so I'm using Euler's method to approximate it.

However, no matter how many timesteps I give Euler's method or what initial conditions I set, it just spits out coordinates of $$(0,0,0)$$ in the second timestep and then random strange numbers. I'm not sure if it's because the speed is so high or if there is a problem with my code. I'd like to figure out the trajectory of particles before they hit the ground and where on the ground they hit.

Here's the full code I'm using, $$x,y,z$$ are the coordinates, $$v$$ is velocity, $$B$$ is the magnetic vector field $$\mathbf{B}=\frac{\mathrm{B_r} V}{4 \pi \sqrt{x^2+y^2+z^2}^5} \left(3zx,3zy,2z^2-x^2-y^2\right)$$ with $$V$$ volume and $$B_r$$ magnetic flux, and $$F$$ is the force (which is a second order DE, given by Lorentz force law). All units are meters.

import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import odeint

c = -299792458
q = -1
m = 0.002
v_x = 0
v_y = 0
v_z = -c*0.994

x_i,y_i,z_i = 1,1,100
v_xi,v_yi,v_zi = 1,1,v_z

h = 0.05
timesteps = list(range(1,100))
x,y,z = x_i,y_i,z_i
v_xi,v_yi,z_yi = v_xi,v_yi,v_zi

height = 1

B_r = 12
Vol = 0.790524
B = (B_r*Vol)/(4*np.pi*np.sqrt(x**2 + y**2 + z**2))
B_x = B*(3*z*x)
B_y = B*(3*z*y)
B_z = B*(2*z**2 - x**2 - y**2)

def F_x(x,y,z):
return((q/(m*c))*(v_y*B_z - v_z*B_y))

def F_y(x,y,z):
return((q/(m*c))*(v_z*B_x - v_x*B_z))

def F_z(x,y,z):
return((q/(m*c))*(v_x*B_y - v_y*B_x))

s = np.zeros((len(timesteps)+1,3))
s[0,:] = [x_i,y_i,z_i]

l = np.zeros((len(timesteps)+1,3))
l[0,:] = [v_xi,v_yi,v_zi]

for i in range(1,len(timesteps)):
l[i+1,0] = l[i-1,0] + h*s[i,1]
s[i+1,0] = s[i-1,0] + h*F_x(l[i,0],l[i,1],l[i,2])
s[i+1,1] = s[i-1,1] + h*F_y(l[i,0],l[i,1],l[i,2])
s[i+1,2] = s[i-1,2] + h*F_z(l[i,0],l[i,1],l[i,2])

print(s[0:10,0:10])

• Why do you have $s^{i+1} = s^{i-1} +hF(l^i)$? At step one, $l^1=0$ afaik, since you only initialized $l^0$ and then setb$l^2 = l^0$. Explicit Euiler is $u^{k+1} = u^k +h f(t_k, u^k)$ and not $u^{k+1} = u^{k-1} + hf(t_k, u^k)$. I think your indices are off by one. But even then I think with that large $c$ and $h$ you will violate stability. Also I am not sure how python works but are $B$ and the rest really functions? They look like constants. Jul 15, 2023 at 10:02
• Yea I think I need to put the B explicitly in the F to turn them into functions. Do you know any alternatives to Euler's method which would work better here? Jul 15, 2023 at 10:50
• You have a speed of 3e+8 to cover a distance of 1e+2, this takes time 3e-7. Your h=0.05 is much too large. Jul 15, 2023 at 11:02
• Implicit Euler would handle large time steps better. There are also multistep methods. I suggest implementing explicit Euler and implicit Euler on a toy example for which you know the solution first though in order to get some practice, also wrt how you program things in python. You also didn't write your evolution equation, although it looks like you wanted to solve $d_{tt} p(t) = F(p(t))$, then substitute $d_t p(t) = v(t)$ resulting in the equation $d_t v(t) = F(p(t))$. Then $p^{k+1} = p^k + h v^k$ and $v^{k+1} = v^k + h F(p^k)$. I do not believe that this is what your code does however. Jul 15, 2023 at 11:06
• I believe you switched position and velocity and you also update only one coordinate here for some reason l[i+1,0] = l[i-1,0] + h*s[i,1], I think you should update all three. My suggestion would be for you to write out the evolution equations explicitly, then discretise along time explicitly, and only after that thinking about implementation details. Jul 15, 2023 at 11:12

I tried to modify your code as little as possible and include the comments by @lightxbulb. I changed the indices in the time stepping loop and modified $$F$$ and $$B$$ so that they are being updated in the time stepping loop.

  import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import odeint

c = -299792458
q = -1
m = 0.002
v_x = 0
v_y = 0
v_z = -c*0.994

x_i,y_i,z_i = 1,1,100
v_xi,v_yi,v_zi = 1,1,v_z

h = 1e-8
timesteps = list(range(1,int(4e5)))
x,y,z = x_i,y_i,z_i
v_xi,v_yi,z_yi = v_xi,v_yi,v_zi

height = 1

B_r = 12
Vol = 0.790524
B = lambda x,y,z : (B_r*Vol)/(4*np.pi*np.sqrt(x**2 + y**2 + z**2))
B_x = lambda x,y,z : B(x,y,z)*(3*z*x)
B_y = lambda x,y,z : B(x,y,z)*(3*z*y)
B_z = lambda x,y,z : B(x,y,z)*(2*z**2 - x**2 - y**2)

def F_x(x,y,z,v_x,v_y,v_z):
return((q/(m*c))*(v_y*B_z(x,y,z) - v_z*B_y(x,y,z)))

def F_y(x,y,z,v_x,v_y,v_z):
return((q/(m*c))*(v_z*B_x(x,y,z) - v_x*B_z(x,y,z)))

def F_z(x,y,z,v_x,v_y,v_z):
return((q/(m*c))*(v_x*B_y(x,y,z) - v_y*B_x(x,y,z)))

p = np.zeros((len(timesteps)+1,3))
p[0,:] = [x_i,y_i,z_i]

v = np.zeros((len(timesteps)+1,3))
v[0,:] = [v_xi,v_yi,v_zi]

for i in range(len(timesteps)):
p[i+1,:] = p[i,:] + h*v[i,:]
v[i+1,0] = v[i,0] + h*F_x(*p[i,:],*v[i,:])
v[i+1,1] = v[i,1] + h*F_y(*p[i,:],*v[i,:])
v[i+1,2] = v[i,2] + h*F_z(*p[i,:],*v[i,:])

plt.title("Position")
plt.plot(p[:,0], p[:,1])
plt.savefig("position.png")


In the code I used Python's lambda functions as a way to define functions in one line. Furthermore,

*p[i,:],*v[i,:]


uses the unpacking operator to abbreviate

p[i,0],p[i,1],p[i,2],v[i,0],v[i,1],v[i,2]


Now I get the solution trajectory for x- and y-positions:

On page 53 in Princeton Companion Applied Math Book there is a nice comparison of the Forward Euler, Backward Euler, Trapezoidal rule and Leapfrog method for the ODE $$\partial_t u = v$$ and $$\partial_t v = -u$$. For the initial condition $$u(0) = 1$$ and $$v(0) = 0$$, the solution is $$u(t) = \cos(t)$$ and $$v(t) = -\sin(t)$$. If you plot $$(u(t), v(t))$$ in the xy-plane, the analytical solution describes a circle of radius 1. With the forward Euler method (as in your code), the solution spirals outward away from the circle. With the backward Euler method, the solution gets damped and spirals inward towards the origin. For the Trapezoidal rule and the Leapfrog method the solution roughly follows the circle. Therefore, I would recommend either using the forward Euler method with a strict time step restriction or even better use a higher order or even energy preserving time stepping method. Personally, I would go for the trapezoidal rule / Crank-Nicolson scheme which for the ODE $$\partial_t u = F(u)$$ reads $$\frac{u^{k+1} - u^{k}}{h} = \frac{1}{2}(F(u^{k+1})+F(u^{k}))$$ or equivalently $$u^{k+1} - \frac{1}{2}F(u^{k+1}) = u^{k} + \frac{1}{2}F(u^{k})$$. But this is just my personal preference.

• Isn't position and velocity still switched? Or did I misunderstand the evolution equation they are trying to solve? My understanding is that you really ought to have $s$ being velocity and $l$ position in the Euler update, and the initial conditions are then also swapped. Jul 15, 2023 at 12:35
• Yes, you are right I think I forgot to change this. I think that the evolution equations are $\partial_t p = v$ and $\partial_t v = \frac{q}{mc}v \times B(p)$. Jul 15, 2023 at 13:09
• The new code should then be the following in my opinion s[i+1,:] = s[i,:] + h*l[i,:] and l[i+1,0] = l[i,0] + h*F_x(*s[i,:],*l[i,:]) (resp. y and z). Am I correct? Jul 15, 2023 at 13:10
• Yes, I believe so, at least that's what I get from the evolution equation I wrote. I would also rename those to $p$ and $v$ instead of $s$ and $l$, even if just to be consistent with the naming of the initial condition variables $v_x, v_y, v_z$. Jul 15, 2023 at 13:38
• Thank you @lightxbulb. I fixed this now in the code and in my answer above. Jul 15, 2023 at 13:50