Wie man Impuls erhält, wenn man den Weg eines geladenen Teilchens durch ein Magnetfeld numerisch integriert

Ich habe das unten angehängte Python-Skript, das die Flugbahn eines geladenen Teilchens in einem statischen, gleichmäßigen Magnetfeld verfolgen soll. Es ist sehr einfach, die augenblickliche Kraft auf das Teilchen zu berechnen (es ist nur das Kreuzprodukt seiner Geschwindigkeit mit dem Magnetfeld B), aber leider erhält mein supereinfacher Versuch, die Auswirkungen dieser Kraft numerisch zu integrieren, keinen Impuls .

Ich verstehe, warum es nicht geht. Wenn ich mir vorstelle, dass bei t=0 der Geschwindigkeitsvektor [1,0,0] ist (mit B=[0,0,1]), dann wird die augenblicklich aufgebrachte Kraft so etwas wie F_mag = [0,a,0] sein ], aber die numerische Integration wendet diese Kraft für mehr als eine infinitesimale Periode an, und so wird die Geschwindigkeit beim nächsten Zeitschritt [1,a,0] sein, wobei a ein Wert ungleich Null ist, der eindeutig eine größere Größe hat als [ 1,0,0]. Jeder Schrittimpuls erhöht sich um einen Faktor, der proportional zur Geschwindigkeit ist. Der Impuls wird nicht konserviert --- exponentielle Explosion (siehe Bild unten)!

Ich habe diesen interessanten Artikel gelesen , der sich anscheinend mit meinem Problem befasst. Es funktioniert durch eine Möglichkeit, die nächste Position des Partikels zu berechnen, wenn seine Position in den vorherigen zwei Zeitschritten gegeben ist. Aber ich möchte die Kraft ausarbeiten , die das Magnetfeld auf Partikel ausübt (da ich schließlich andere Kräfte in mein Modell einbeziehen möchte).

Gibt es eine einfache Möglichkeit, die vom Magnetfeld ausgeübte Kraft so zu berechnen, dass der Impuls erhalten bleibt? Ich kann die Schrittgröße kleiner machen, aber die Dynamik wächst immer noch exponentiell – nicht ideal! Ich könnte auch einfach zu Runge-Kutta statt einfacher Euler-Integration wechseln, aber ich dachte, dass es hier vielleicht die Möglichkeit gibt, etwas Klügeres zu tun! Vielen Dank im Voraus für alle Vorschläge, die Sie haben könnten!

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