4
$\begingroup$

I have the python script attached below, which is intended to track the trajectory of a charged particle in a static, uniform magnetic field. It is very simple to calculate the instantaneous force upon the particle (it is just the cross-product of its velocity with the magnetic field, B), but unfortunately my super-simple attempt at numerically integrating the effects of this force does not conserve momentum.

I understand why it doesn't. If I imagine that at t=0, the velocity vector is [1,0,0] (with B=[0,0,1]), then the instantaneous force applied will be something like F_mag = [0,a,0], but numerical integration applies this force for more than an infinitessimal period and so at the next timestep, the velocity will be [1,a,0], where a is some non-zero value, which clearly has a greater magnitude than [1,0,0]. Every step momentum will increase by a factor proportional to velocity. Momentum is not conserved --- exponential explosion (see image below)!

I have studied this interesting article that appears to be addressing my issue. It works through a way to calculate the next position of the particle given its position in the previous two time steps. But I want to work out the force that the magnetic field is applying to particle (as I want to eventually include other forces in my model).

Is there a straight-forward way to calculate the force applied by the magnetic field in such a way as so conserve momentum? I can make the step size smaller, but still momentum grows exponentially -- not ideal! I could also just switch over to Runge-Kutta rather than simple Euler integration, but I thought there might be the opportunity here to do something more clever! Thanks in advance for any suggestions you might have!

from pylab import * 

DT = 0.1 
N_ITS = int(100.0/DT)

pos    = array([0.0,0.0,0.0],'f')
vel    = array([1.0,0.0,0.0],'f') 
mass   = 1.0 
charge = 1.0

# uniform magnetic field
B = array([0.0,0.0,1.0])

# track momentum for plotting
mom_h = []
# track position for plotting
pos_h = []

for _ in xrange(N_ITS) :
    ## total momentum
    mom = linalg.norm(vel)
    mom_h.append(mom)
    pos_h.append(array(pos))

    ### fixed magnetic field
    F_mag = array(-cross(vel,B))

    ### apply force
    vel += DT*F_mag/mass

    ### update position
    pos += vel*DT

figure(figsize=(6,10))
subplot2grid((2,1),(0,0)) 
title('position') 
pos_h = array(pos_h)    
plot(pos_h[:,0],pos_h[:,1])

subplot2grid((2,1),(1,0))
title('momentum') 
plot(mom_h)

show()

enter image description here

$\endgroup$
  • $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ – Qmechanic Apr 18 '16 at 6:40
  • $\begingroup$ @Qmechanic Given the role of momentum conservation in the question, it might be on topic here. I'm looking at this meta question and I think this is very close to the boundary between on topic and off topic, but it doesn't quite hit any of the criteria we identified as being off topic. (I've gone back and forth on this decision over the years.) $\endgroup$ – David Z Apr 18 '16 at 7:03
  • $\begingroup$ I deleted some comments that were answering the question. $\endgroup$ – David Z Apr 18 '16 at 7:04
  • $\begingroup$ @DavidZ Um.. why? This is quite annoying. I saw some useful information posted here in response to my question, but now I can't access that information any more? Who does it benefit to delete those comments? $\endgroup$ – weemattisnot Apr 18 '16 at 21:28
  • $\begingroup$ @weemattisnot It helps the site stay clean. Those people should have posted answers, not comments, and they are still welcome to do so. $\endgroup$ – David Z Apr 19 '16 at 9:21
2
$\begingroup$

With forward Euler, you're simply out of luck. Because of the way the scheme is constructed, you always make an error in momentum in the same direction, which then compounds, leading to exponential behavior. You need a symplectic integrator, the simplest being leapfrog, or any other Verlet integrator. These still won't conserve momentum on a timestep to timestep basis, but will on average over an "orbit", i.e. a turn around the helix for your particle. From there, you can set the timestep to keep the momentum error below some tolerance, if you wish.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.