Python Verlet algorithm and Lorentz force trajectory

  • Thread starter Arman777
  • Start date
1,506
110
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.
 

BvU

Science Advisor
Homework Helper
12,021
2,645
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 !
 
1,506
110
Well I find a appraoch to solve it but I need to turn them into the code form
 

BvU

Science Advisor
Homework Helper
12,021
2,645
Assume you have a subroutine available to do the integration. What are you going to integrate ?
 
1,506
110

anorlunda

Mentor
Insights Author
Gold Member
6,791
3,748
Are you sure this is not homework? What is it?
 
32,519
4,237
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.
 
1,506
110
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

Science Advisor
Homework Helper
12,021
2,645
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.
 
1,506
110
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

Science Advisor
Homework Helper
12,021
2,645
##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.
 
1,506
110
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
12,021
2,645
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 ?
 
1,506
110
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
12,021
2,645
Well I am trying but its bit of hard
I sympathize :wink: . I suppose no one claimed it would be a breeze.

However: if it were super easy you wouldn't learn much :cool:

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 :rolleyes:
 
1,506
110
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
12,021
2,645
Well, I got what I asked for ... o0)
(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
12,021
2,645
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
 
1,506
110
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
12,021
2,645
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'

upload_2018-4-30_17-40-46.png


(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 ?
My bad. I meant:
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

1,506
110
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
 
1,506
110
And Euler is simpler yes but not accurate..sadly
 

BvU

Science Advisor
Homework Helper
12,021
2,645
I know, but that we can deal with once it's running as we intend...
 

BvU

Science Advisor
Homework Helper
12,021
2,645
I used excel -- my python is still in statu nascendi :rolleyes:

is a place to start, perhaps
 
Last edited:

Want to reply to this thread?

"Verlet algorithm and Lorentz force trajectory" You must log in or register to reply here.

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
Top