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