Как сохранить импульс при численном интегрировании пути заряженной частицы через магнитное поле

У меня есть скрипт Python, прикрепленный ниже, который предназначен для отслеживания траектории заряженной частицы в статическом однородном магнитном поле. Очень просто вычислить мгновенную силу, действующую на частицу (это просто произведение ее скорости на магнитное поле B), но, к сожалению, моя сверхпростая попытка численного интегрирования эффектов этой силы не сохраняет импульс. .

Я понимаю, почему это не так. Если я представлю, что при t=0 вектор скорости равен [1,0,0] (при B=[0,0,1]), то приложенная мгновенная сила будет примерно такой: F_mag = [0,a,0 ], но численное интегрирование применяет эту силу в течение более чем бесконечно малого периода, и поэтому на следующем временном шаге скорость будет [1,a,0], где a — некоторое ненулевое значение, которое явно имеет большую величину, чем [ 1,0,0]. С каждым шагом импульс будет увеличиваться на коэффициент, пропорциональный скорости. Импульс не сохраняется --- экспоненциальный взрыв (см. изображение ниже)!

Я изучил эту интересную статью , которая, похоже, касается моей проблемы. Он работает через способ вычисления следующего положения частицы с учетом ее положения на предыдущих двух временных шагах. Но я хочу определить силу , с которой магнитное поле действует на частицу (поскольку я хочу в конечном итоге включить в свою модель другие силы).

Есть ли прямой способ рассчитать силу, приложенную магнитным полем, таким образом, чтобы сохранить импульс? Я могу уменьшить размер шага, но импульс все равно будет расти экспоненциально — не идеально! Я также мог бы просто переключиться на Рунге-Кутту, а не на простое интегрирование по Эйлеру, но я подумал, что здесь может быть возможность сделать что-то более умное! Заранее спасибо за любые ваши предложения!

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