Verlet algorithm and Lorentz force trajectory

In summary: Verlet algorithm?First, in order to approximate the ##a_x(x,v,t)## to ##q/m[E_x(x,t)-G_x(v_{n-1/2})]##, you would use the following equation:v = v_0 - Δt/2a_0Next, to solve for x_n+1, you would use the following equation:x_n+1 = x_n + Δt[v_n+1/2]
  • #1
Arman777
Insights Author
Gold Member
2,168
192
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.
 
Technology news on Phys.org
  • #2
Arman777 said:
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 !
 
  • #3
Well I find a appraoch to solve it but I need to turn them into the code form
 
  • #4
Assume you have a subroutine available to do the integration. What are you going to integrate ?
 
  • #6
Are you sure this is not homework? What is it?
 
  • #7
Arman777 said:
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

anorlunda said:
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.
 
  • #8
Mark44 said:
ut at this point @Arman777 seems to be looking for guidance on how to approach the problem.
Yes exactly.
Mark44 said:
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 doesn't matter for this stage..since I still don't know how to approach the problem.

Mark44 said:
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) ?
 
  • #9
Arman777 said:
since I still don't 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.
 
  • #10
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...
 
  • #11
Arman777 said:
##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, ...).

Arman777 said:
##\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)##

Arman777 said:
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.
 
  • #12
BvU said:
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)
BvU said:
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.

BvU said:
Your plan is one way. Another is to start with an Euler step and then change to your favorite method.
I don't 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.
 
  • #13
Arman777 said:
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.

Arman777 said:
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 ?
 
  • #14
BvU said:
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 couldn't finish it and gives me error. I should probably ask it my teacher
 
  • #15
Arman777 said:
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 .

Arman777 said:
I can share my code but its not even wrong since I couldn't finish it
Might help. At least we can comment a bit
Arman777 said:
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:
 
  • #16
BvU said:
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..
BvU said:
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()
 
  • #17
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?).
 
  • #18
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
 
  • #19
BvU said:
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 doesn't much matter I guess (?)

BvU said:
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##)

BvU said:
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 ?

BvU said:
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

BvU said:
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(')
 
  • #20
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.

Arman777 said:
it doesn't 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.

Arman777 said:
Not really ? Then how can I increase the time ?
My bad. I meant:
VHx=VHx + dt * ( E[0] + y[0] ) * ( q / m )
Remember:
BvU said:
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.

Arman777 said:
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

  • upload_2018-4-30_17-40-46.png
    upload_2018-4-30_17-40-46.png
    3.4 KB · Views: 549
  • #21
BvU said:
Aren't we having fun ?
ahh no not reall :/
BvU said:
For n = 0 you can simply write
VH = Vi
y = cross(VH,B)before entering the loop.
hmm I ll try that
BvU said:
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
 
  • #22
And Euler is simpler yes but not accurate..sadly
 
  • #23
I know, but that we can deal with once it's running as we intend...
 
  • #24
How can I make a graph ?
 
  • #25
I used excel -- my python is still in statu nascendi :rolleyes:

is a place to start, perhaps
 
Last edited:
  • #26
BvU said:
I used excel -- my python is still in statu nascendi :rolleyes:

is a place to start, perhaps

oh thanks a lot :smile:
 
  • #27
I finally finished the code if anyone of you wants to see the code pm me and I ll send it
 
  • #28
Is it too long to post here ? Still interested to 'hear how you are doing (and @isukatphysics69 might find some goodies in this thread too).
 
  • #29
Not really but I prefer through pm
 
  • Like
Likes BvU

1. What is the Verlet algorithm?

The Verlet algorithm is a numerical method used to solve equations of motion in classical mechanics. It is commonly used in molecular dynamics simulations to calculate the positions and velocities of particles over time.

2. How does the Verlet algorithm work?

The Verlet algorithm uses a time-stepping approach to calculate the positions and velocities of particles. It involves using the current position and velocity of a particle to predict its position at the next time step, and then using this predicted position to update the velocity. This process is repeated for each time step, resulting in a trajectory of the particle's motion over time.

3. What is the Lorentz force in classical mechanics?

The Lorentz force is the force experienced by a charged particle moving in an electromagnetic field. It is given by the equation F = q(E + v x B), where q is the charge of the particle, E is the electric field, v is the velocity of the particle, and B is the magnetic field.

4. How is the Lorentz force incorporated into the Verlet algorithm?

The Lorentz force is incorporated into the Verlet algorithm by adding an additional force term to the equations of motion. This force term is calculated using the current position and velocity of the particle, as well as the electric and magnetic fields at that point in space.

5. What are some applications of the Verlet algorithm and Lorentz force trajectory?

The Verlet algorithm and Lorentz force trajectory have many applications in physics and engineering. They are commonly used in molecular dynamics simulations to study the behavior of particles in a variety of systems, such as gases, liquids, and solids. They are also used in simulations of particle accelerators, plasma physics, and astrophysics.

Similar threads

  • Programming and Computer Science
Replies
2
Views
2K
  • Programming and Computer Science
Replies
2
Views
7K
  • Programming and Computer Science
Replies
19
Views
2K
  • Programming and Computer Science
Replies
2
Views
856
  • Introductory Physics Homework Help
Replies
3
Views
757
  • Programming and Computer Science
Replies
3
Views
1K
  • Introductory Physics Homework Help
Replies
3
Views
365
  • Differential Geometry
Replies
1
Views
1K
  • Programming and Computer Science
Replies
4
Views
1K
  • Introductory Physics Homework Help
Replies
7
Views
931
Back
Top