Python Verlet algorithm and Lorentz force trajectory

[Arman777]

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]

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]

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

 

[BvU]

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

 

[Arman777]

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

velocity..?
https://pdfs.semanticscholar.org/1261/d723e67a771a4cb93bed9f4a38177754735b.pdf
in here page 8 ı have my equations..my initial conditions are x=0 and v=0 so I need the last 2 equations right ?

edit:I also didnt understand the v-hat or x-hat thing ?

 

[anorlunda]

Are you sure this is not homework? What is it?

 

[Mark44]

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]

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 doesn't matter for this stage..since I still don't 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]

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.

 

[Arman777]

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

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

 

[BvU]

$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]

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 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.

 

[BvU]

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]

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

 

[BvU]

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 couldn't 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]

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

#!/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]

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]

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]

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 (?)

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,

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

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]

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]

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 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.

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.

 

[Arman777]

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]

And Euler is simpler yes but not accurate..sadly

 

[BvU]

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

 

[Arman777]

How can I make a graph ?

 

[BvU]

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

is a place to start, perhaps

 

[Arman777]

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

is a place to start, perhaps

oh thanks a lot

 

[Arman777]

I finally finished the code if anyone of you wants to see the code pm me and I ll send it

 

[BvU]

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).

 

[Arman777]

Not really but I prefer through pm