# Using ODE to plot particle-motion with scipy.integrate.solve_ivp

## My Problem:

A positively charged particle (mass = 2 * 10-27 kg) is moving along the x-axis. It is travelling in a homogenous magnetic field such that the field axis in z-direction. The energy of the particle is 2 MeV and B = 4 T. Use a ODE solver to plot the motion of the particle for 1 microseconds.

## My attempt of solving the problem

#### Note that I have used question marks where I am unsure.

import numpy as np
from scipy.integrate import solve_ivp

initialZ = [?, ?, ?, ?, ?, ?] # = [positionX, positionY, positionZ, velocityX, velocityY, velocityZ]
t0 = 0
tf = 1*(10**-6) # 1 microsecond = 1*10^6 seconds
times = (t0, tf)

def ivf(t, Z):
x, y, z = Z[0], Z[1], Z[2]
u, v, w = Z[3], Z[4], Z[5]
return np.array([u, v, w, ?, ?, ?])

s = solve_ivp(ivf, times, initialZ)


## My question

#### what should the question marks (?) in the code be replaced with?

I have tried to solve the ODE as an initial value problem. I tried to determine inital velocity by equating the lorentz force and centripetal force. I am very new to differential equations and it is therefore difficult to know If I am doing things in the correct way. (note that the first three values of my initialZ vector represents positions x, y and z, and the last three values represent velocity in the x, y and z direction). I am grateful for any help or guidance.
• You can use zeros for the initial coordinates. The initial velocity is stated to be in the x direction, so the other components are zero. Your function ivf() returns the time derivatives of vector Z, so here you'd use your 2nd Newton's law. The only caveat here is to check if at this energy your particle needs to be treated relativistically. Jul 2 at 7:25
• You first need to write your system of equations as a system of first-order differential equations. That will give you 6 equations. Those are the 6 initial conditions and the 6 variables inside ivf. Also, maybe converting your problem to nondimensional form might prove useful. Jul 2 at 16:34
• Cross-post from stackoverflow.com/questions/68217126/…, a little more on-topic here. Jul 3 at 6:32

Your particle is a rounded proton (mass m = 2e-27 kg instead of 1.672e-27 kg). The equation of motion is $$\dot x=v,~~~ m\dot v = q\,v\times B,$$ where $$B=(0,0,B_z)$$ with $$B_z=4T=4N/(m\,A)$$ and $$q=1e=1.602·10^{-19} C$$, $$C=A\,s$$

This then gives for the acceleration

m=2e-27
e_charge = 1.6e-19
q=+1*e_charge
Bz = 4

ax = q/m*vy*Bz; ay = -q/m*vx*Bz; az = 0


For the initial conditions, without the field the particle would move along the $$x$$ axis with constant speed $$v_x$$, $$v_y=v_z=0$$. As the system is translation invariant, set the position to the origin at time $$t=0$$. The speed is then given via the kinetic energy formula $$\frac12mv_x^2=E=2MeV\implies v_x=\sqrt{\frac{4MeV}{m}}$$ or

vx = 2e+3*sqrt(e_charge/m)