# Python Verlet algorithm and Lorentz force trajectory

#### Arman777

Gold Member
I need to write a code to calculate the trajectory of the particle under lorentz force.Since the position depended equations are hard to solve I ll use a code, how can I appraoch this problem. I should use veloctiy-verlet algorithm or any suggestions ? You should consider that lorentz force is a velocity dependent force.

Related Programming and Computer Science News on Phys.org

#### BvU

Homework Helper
I need to write a code
So, is this homework ?

Anyway, why not concentrate a bit more on what has to be done, and a little less on which method of integration to use ? In a well-written code you can even make that a choice between various methods !

#### Arman777

Gold Member
Well I find a appraoch to solve it but I need to turn them into the code form

#### BvU

Homework Helper
Assume you have a subroutine available to do the integration. What are you going to integrate ?

Gold Member

#### anorlunda

Mentor
Gold Member
Are you sure this is not homework? What is it?

#### Mark44

Mentor
edit:I also didnt understand the v-hat or x-hat thing ?
From what I can tell, $\hat x$ and $\hat v$ are estimates of position and velocity

Are you sure this is not homework? What is it?
It might or might not be homework, but at this point @Arman777 seems to be looking for guidance on how to approach the problem. When things get to the stage of writing some code, it should be a new thread in the Homework section, under Engineering & Comp Sci.

#### Arman777

Gold Member
ut at this point @Arman777 seems to be looking for guidance on how to approach the problem.
Yes exactly.
When things get to the stage of writing some code, it should be a new thread in the Homework section, under Engineering & Comp Sci.
I agree..I mean I also think it doesnt matter for this stage..since I still dont know how to approach the problem.

From what I can tell, ^xx^\hat x and ^vv^\hat v are estimates of position and velocity
Hmm Which ones should I use then ?

in leapfrog it says inital conditions should be like in eqn(8) ( page 3 ) and also we have $x_0$
so using these 2 values as inital values I ll try to write the first 2 equations in (26) ?

#### BvU

Homework Helper
since I still dont know how to approach the problem
I, in turn, can't guess how to help you effectively. Why not set up a scenario, define some symbols and list the information considered as input.

#### Arman777

Gold Member
Okay let me try,
So I ll use leapfrog/varlet algorithm.
In leapfrog method for initial values we need $v_{-1/2}$ and $x_0$. Hence I can write
$v_{-1/2}=v_0-Δt/2a_0$
$v_0$ is given as $0$ , but $a_0$ should be also given to us right ?

and $x_0$ also given as $0$

In the Verlet/leapfrog method for damped systems we can write the system equation as
$a(x,v,t)=1/m[F(x,t)-G(v)]$

Since I am doing lorentz force it will look like
$\vec{a}(x,v,t)=q/m[\vec{E}(x,t)-\vec{G}(v)]$ where $G(v)=\vec{v}×\vec{B}$
or in component-wise
$a_x(x,v,t)=q/m[E_x(x,t)-G_x(v)]$
$a_y(x,v,t)=q/m[E_y(x,t)-G_y(v)]$
$a_z(x,v,t)=q/m[E_z(x,t)-G_z(v)]$

so for example
let us think one component x

it says that we should approximate the $a_x(x,v,t)$ to $q/m[E_x(x,t)-G_x(v_{n-1/2})]$ then we have
$v_{n+1/2}=v_{n-1/2}+Δt[E_x(x,t)-G_x(v_{n-1/2})]$

In here we should set $v_{n-1/2}=v_{-1/2}=v_0-Δt/2a_0$ I guess but I am not sure what should I do with $G_x(v_{n-1/2})$

this is the first problem

Next since I know the $v_{n+1/2}$, I can use it in the

$x_{n+1}=x_n+Δt[v_{n+1/2}]$ equation. Now I can solve the problem I guess ?

In the code process I should set a $dt=0.0001$ or something like that and multiply it by $n$ right ?
so I mean something like

Python:
dt=0.0001
for n in range(1000):
v=...n*dt...

#### BvU

Homework Helper
$v_0$ is given as 0 , but $a_0$ should be also given to us right ?
It would seem more logical to me that you are given the fields and then you have to calculate $\vec a$, so also at $t=0$.
$v_0=0$ is possible, but $> 0$ may be more usual -- depends on the application (high-energy physics, electron tube, ...).

$\vec a(\vec x, \vec v,t)=q/m[\vec E(\vec x,t)−\vec G(v)] For the Lorentz force I'm more used to$F_L = q(\vec E + \vec v \times \vec B)$not sure what should I do with G_x(v_{n-1/2}) Your plan is one way. Another is to start with an Euler step and then change to your favorite method. #### Arman777 Gold Member It would seem more logical to me that you are given the fields and then you have to calculate →aa→\vec a, so also at t=0t=0t=0. v0=0v0=0v_0=0 is possible, but >0>0> 0 may be more usual -- depends on the application (high-energy physics, electron tube, ...). Well its how the question asked yes E and B should be given but it asks about the trajectory so in other words x(t) For the Lorentz force I'm more used to FL=q(→E+→v×→B)FL=q(E→+v→×B→)F_L = q(\vec E + \vec v \times \vec B) okay so ? I just used the notation that the link has. Your plan is one way. Another is to start with an Euler step and then change to your favorite method. I dont think that would work. Since I lorentz force is a velocity depended force and it would be a waste of time I guess to study that and then again come to the this system again. #### BvU Science Advisor Homework Helper Since the Lorentz force is a velocity dependent force It is velocity dependent, whether you go forward or backward. Integration method shouldn't make a difference. Once you have done one (half size) Euler step (backward -- as you propose -- or forward) you have the values of the variables you need to proceed the integrations following one of the other schemes. E and B should be given but it asks about the trajectory so in other words x(t) So we understand that you have$\vec E(\vec r, t),\ \vec B(\vec r, t),\ \vec r_0,\ \vec v_0\ $given. (or perhaps the fields do not depend on time ?). It's a 3D problem, correct ? So how many integrations do you have to keep track of ? re your python code: ' v=...n*dt... ' where would you need n*dt ? #### Arman777 Gold Member It is velocity dependent, whether you go forward or backward. Integration method shouldn't make a difference. Once you have done one (half size) Euler step (backward -- as you propose -- or forward) you have the values of the variables you need to proceed the integrations following one of the other schemes. So we understand that you have$\vec E(\vec r, t),\ \vec B(\vec r, t),\ \vec r_0,\ \vec v_0\ $given. (or perhaps the fields do not depend on time ?). It's a 3D problem, correct ? So how many integrations do you have to keep track of ? re your python code: ' v=...n*dt... ' where would you need n*dt ? Well I am trying but its bit of hard..I need to recall the function and its values everytime so its hard..I can share my code but its not even wrong since I couldnt finish it and gives me error. I should probably ask it my teacher #### BvU Science Advisor Homework Helper Well I am trying but its bit of hard I sympathize . I suppose no one claimed it would be a breeze. However: if it were super easy you wouldn't learn much Nevertheless, one shouldn't make things more complicated than necessary -- which is why I tried to formulate a complete problem statement iin #13. Are$\vec B$and$\vec E$in fact constant in space ? In time? I notice your other thread where a lot of groundwork seems to be laid already (or is it totally unrelated ?) and you had some excellent help form @Orodruin and @Ray Vickson . I can share my code but its not even wrong since I couldnt finish it Might help. At least we can comment a bit I should probably ask it my teacher Always a good idea. I hope his comments are along the lines of the comments you had here #### Arman777 Gold Member I notice your other thread where a lot of groundwork seems to be laid already (or is it totally unrelated ?) and you had some excellent help form @Orodruin and @Ray Vickson . It was related but then I learned that we cannot do like that.. Are →BB→\vec B and →EE→\vec E in fact constant in space ? In time? yes they are constant okay I ll share my code Python: #!/usr/bin/python print("Lorentz Force Trajectory Calculator") q=-1.6e-19 dt=0.0001 m=9.10938e-31 file=open("Dreams.txt","w") Ex,Ey,Ez=(input("Enter the X,Y,Z components of the Electric field respectively in V/m: ")).split(",") Bx,By,Bz=(input("Enter the X,Y,Z components of the Magnetic field in Tesla: ")).split(",") E0=[Ex,Ey,By] B0=[Bx,By,Bz] E=[] B=[] for i in range(3): E.append(float(E0[i])) B.append(float(B0[i])) Vi=[-dt/2,-dt/2,-dt/2] Ri=[0,0,0] P=[] for n in range(1000): def cross(V,B): Px=(V[1]*B[2]-V[2]*B[1]) Py=(V[2]*B[0]-V[0]*B[2]) Pz=(V[0]*B[1]-V[1]*B[0]) P=[Px,Py,Pz] return P y=cross(Vi,B) VHx=Vi[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m) VHy=Vi[1]+(n-1/2)*dt*(E[1]-y[1])*(q/m) VHz=Vi[2]+(n-1/2)*dt*(E[2]-y[2])*(q/m) VH=[VHx,VHy,VHz] Rx=Ri[0]+n*dt*VH[0] Ry=Ri[1]+n*dt*VH[1] Rz=Ri[2]+n*dt*VH[2] file.write("%.3f %4.3f %4.3f\n" %(Rx,Ry,Rz)) y=cross(VH,B) file.close() #### BvU Science Advisor Homework Helper Well, I got what I asked for ... (still struggling to get pycharm debugging to work decently on my desktop ) But commenting is easy: Why is the def cross inside the loop ? Why is Vi = [-dt/2, -dt/2, -dt/2] ? With 'VHx=Vi[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m)' you take a huge step when n is large I suppose you mean Vi[1]=Vi[0]+*dt*(E[0]-y[0])*(q/m) ? Idem R. You want to integrate, i.e. R(t+dt) = R(t) + dt * dR/dt ! And I'd start with q=m=1, E=[1,0,0], B = [0,1,0], Vi = [0,0,1] and see what comes out (a straight line?). #### BvU Science Advisor Homework Helper Oh, and I see where your minus suign comes from. But now you use it twice for B y=cross(Vi,B) has the same arguments with every call ! y=cross(VH,B) is overwritten bfre y is used #### Arman777 Gold Member Why is the def cross inside the loop? It was out at first then I put it inside... I thought after the y=cross(VH,B) line it will define the "cross " again but it doesnt much matter I guess (?) Why is Vi = [-dt/2, -dt/2, -dt/2] ? Our initial velocity is$<0,0,0>$and in the formula it says$v_{-1/2}=v_0-Δt/2a_0$(I assumed$a_0=0$) With 'VHx=Vi[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m)' you take a huge step when n is large I suppose you mean Vi[1]=Vi[0]+*dt*(E[0]-y[0])*(q/m) ? Not really ? Then how can I increase the time ? And I'd start with q=m=1, E=[1,0,0], B = [0,1,0], Vi = [0,0,1] and see what comes out (a straight line?). I have no idea what would it looks like Oh, and I see where your minus sign comes from. But now you use it twice for B y=cross(Vi, B) has the same arguments with every call! y=cross(VH,B) is overwritten before y is used That's the real problem kind of. Actually, I think the only problem in my code is how to change the variables...for example at the end of the first loop we have VH vector right, so in the second loop this vector should be "crossed" by B. Here what I am trying to do So at first, it is like, Python: for n=0 y=cross(Vi,B) VHx=Vi[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m) VHy=Vi[1]+(n-1/2)*dt*(E[1]-y[1])*(q/m) VHz=Vi[2]+(n-1/2)*dt*(E[2]-y[2])*(q/m) VH=[VHx,VHy,VHz] then it should be like Python: for n=1 y=cross(VH,B) VH'x=VH[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m) VH'y=VH1]+(n-1/2)*dt*(E[1]-y[1])*(q/m) VH'z=VH[2]+(n-1/2)*dt*(E[2]-y[2])*(q/m) VH'=[VH'x,VH'y,VH'z] Python: for n=2 y=cross(VH',B) VH''x=VH'[0]+(n-1/2)*dt*(E[0]-y[0])*(q/m) VH''y=VH'1]+(n-1/2)*dt*(E[1]-y[1])*(q/m) VH''z=VH'[2]+(n-1/2)*dt*(E[2]-y[2])*(q/m) VH''=[VHx'',VH''y,VH''z] Notice the primes(') #### BvU Science Advisor Homework Helper Aren't we having fun ? My debugger works and so does the program (but I'm too dumb to use anything more sophisticated than Euler) with q=m=1, dt=0.007, E=[0,0,0], B=[0, 1, 0], Vi=[0, 0.1, 1] I get a 'helix' (picture is the x-z plane; y is boring) Note: q=1 and m =1 and the other simple numbers make it a lot easier to check results. I increased dt to get a$> 2 \pi ## time.

it doesnt much matter I guess (?)
I suppose a thousand definitions take a thousand times longer than just one ?

Then: if a0 = 0 then dt/2 times a0 can not be nonzero.

Not really ? Then how can I increase the time ?
VHx=VHx + dt * ( E[0] + y[0] ) * ( q / m )
Remember:
You want to integrate, i.e. R(t+dt) = R(t) + dt * dR/dt !
For n = 0 you can simply write
VH = Vi
y = cross(VH,B)​
before entering the loop.

Notice the primes(')
I notice. But I don't think Python likes them. You would need a two-dimensional array. That's wasteful overkill.

My point is you can simply overwrite variables when you don't need the old value any more. That's why I like to start with Euler. For a two-step thing you indeed need to keep some VHx_old idem y z. But after the next step you can re-use these VHx_old etc
Again, the choice of integrator hinders you here when trying to get the right algorithm.

#### Attachments

• 4.4 KB Views: 183

#### Arman777

Gold Member
Aren't we having fun ?
ahh no not reall :/
For n = 0 you can simply write
VH = Vi
y = cross(VH,B)before entering the loop.
hmm I ll try that
I notice. But I don't think Python likes them. You would need a two-dimensional array. That's wasteful overkill.

My point is you can simply overwrite variables when you don't need the old value any more. That's why I like to start with Euler. For a two-step thing you indeed need to keep some VHx_old idem y z. But after the next step you can re-use these VHx_old etc
Again, the choice of integrator hinders you here when trying to get the right algorithm.
I ll look more

#### Arman777

Gold Member
And Euler is simpler yes but not accurate..sadly

#### BvU

Homework Helper
I know, but that we can deal with once it's running as we intend...

#### Arman777

Gold Member
How can I make a graph ?

#### BvU

Homework Helper
I used excel -- my python is still in statu nascendi

is a place to start, perhaps

Last edited:

"Verlet algorithm and Lorentz force trajectory"

### Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving