Cómo conservar el impulso al integrar numéricamente la trayectoria de una partícula cargada a través de un campo magnético

Tengo el script de python adjunto a continuación, que está destinado a rastrear la trayectoria de una partícula cargada en un campo magnético estático y uniforme. Es muy simple calcular la fuerza instantánea sobre la partícula (es solo el producto cruzado de su velocidad con el campo magnético, B), pero desafortunadamente mi intento súper simple de integrar numéricamente los efectos de esta fuerza no conserva el impulso. .

Entiendo por qué no. Si imagino que en t=0, el vector de velocidad es [1,0,0] (con B=[0,0,1]), entonces la fuerza instantánea aplicada será algo así como F_mag = [0,a,0 ], pero la integración numérica aplica esta fuerza durante más de un período infinitesimal y, por lo tanto, en el siguiente paso de tiempo, la velocidad será [1,a,0], donde a es un valor distinto de cero, que claramente tiene una magnitud mayor que [ 1,0,0]. Cada impulso de paso aumentará en un factor proporcional a la velocidad. El impulso no se conserva --- ¡explosión exponencial (vea la imagen a continuación)!

He estudiado este interesante artículo que parece estar abordando mi problema. Funciona a través de una forma de calcular la siguiente posición de la partícula dada su posición en los dos pasos de tiempo anteriores. Pero quiero calcular la fuerza que el campo magnético está aplicando a la partícula (ya que eventualmente quiero incluir otras fuerzas en mi modelo).

¿Existe una forma sencilla de calcular la fuerza aplicada por el campo magnético de tal manera que se conserve el impulso? Puedo hacer que el tamaño del paso sea más pequeño, pero aún así el impulso crece exponencialmente, ¡no es lo ideal! También podría simplemente cambiar a Runge-Kutta en lugar de la simple integración de Euler, ¡pero pensé que podría haber una oportunidad aquí para hacer algo más inteligente! ¡Gracias de antemano por cualquier sugerencia que pueda tener!

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