The following exercise in setting up a numerical integration of
the trajectory of a charged particle in a static magnetic field
has probably been duplicated thousands of times. First the
Lorentz equation is written down:
F=mddtv=eE+v×B
F
m
t
v
e
E
v
B
(1)
and then simplified (we will assume planar motion in a uniform
magnetic field as the first step to get our feet wet) to the
xx-
yy
plane with
E=0
E
0
and
B=
B
0
ez
B
B
0
e
z
ddt
v
x
=
v
y
em
B
0
t
v
x
v
y
e
m
B
0
(2)
and
ddt
v
y
=-
v
x
em
B
0
t
v
y
v
x
e
m
B
0
(3)
Now comes the time to convert this into a numerical scheme,
where the trajectory is sampled at time steps
Δt
Δ
t
.
The numerical methods textbooks make the next step
obvious. First you replace the acceleration with
ddt
v
x
⇒
v
x
'
-
v
x
Δt
t
v
x
v
x
'
v
x
Δ
t
(4)
where
v
x
'
v
x
'
is the new velocity after one time step. After rearrangement, we
get the equations
v
x
'
=
v
x
+
v
y
em
B
0
Δt
v
x
'
v
x
v
y
e
m
B
0
Δ
t
(5)
and
v
y
'
=
v
y
-
v
x
p
v
y
'
v
y
v
x
p
(6)
where we have replaced the dimensionless constant
e
B
0
Δtm
e
B
0
Δ
t
m
with the shorter symbol:
pp.
These are very simple equations to iterate to recursively
generate the trajectory in velocity space by simple repetition.
To get the trajectory in configuration space, we can
simultaneously iterate the equations
x
'
=x+
v
x
Δt
x
'
x
v
x
Δ
t
(7)
and
y
'
=y+
v
y
Δt
y
'
y
v
y
Δ
t
(8)
As thousands of would-be programmers have found, this set of
equations does not quite give the desired behavior. The orbit
is not the closed circle that one gets analytically. Nor is the
energy conserved as required by these differential equations.
So thousands of thinkers have reasoned that these failures are a
natural result of trying to do this problem numerically. The
quantity pp is the cyclotron
frequency times the time step, so it is approximately the
angular step taken by the particle in trying to go around in a
circle. Clearly if this step is substantial, the circle is only
crudely outlined.
This reasoning is reinforced by the standard solution to
numerical inaccuracies: take smaller steps. If we reduce
Δt
Δ
t
(and with it
p
p),
these two errors are in fact reduced and a nicer circle is
transcribed. And there are somewhat more sophisticated ways to
doing the numerical integrations (half-step techniques and
Runge-Kutta strategies) that help but do not cure these
problems. But of course, if the problem is intrinsic to the
numerics, no cure is possible.
It's not an accident that almost every numerical treatment of
some complex problem ends up taking hours to days. These are
simply the longest run times that human operators are willing to
tolerate, and eventually that, and not accuracy, is the limiting
factor. And taking tiny steps to get even that accuracy is
surely a factor in long run times.
We will now proceed to solve this problem exactly!
The real problem with the above algorithm that it is
incorrect. First, let us deal with the coordinates to begin with
instead of setting up a two-step integration. Then the
numerical replacement for the accelerations (Equation 4) are
ddt
v
x
⇒
x
2
-2
x
1
+
x
0
Δt2
t
v
x
x
2
2
x
1
x
0
Δ
t
2
(9)
where we have changed notation and given the coordinates
numerical indices to indicate successive values (0,1,2). Thus
at any step, we have the values
x
0
x
0
and
x
1
x
1
and will use them to calculate the new value,
x
2
x
2
.
This is just the differences of two successive velocities,
divided by
Δ
t
Δt,
and is the standard replacement. To run the program we simply
replace
x
0
x
0
with
x
1
x
1
and
x
1
x
1
with
x
2
x
2
,
and repeat the calculation as many times as needed.
But now we have to calculate what this acceleration is equal to,
"
v
y
p
v
y
p
".
Well the constant pp is no problem,
but what about the yy-veclocity? It
must be some combination of the corresponding coordinates
y
0
y
0
,
y
1
y
1
,
and
y
2
y
2
.
So let's look at the acceleration term. This is most accurately
calculated at what xx and yy position? By symmetry
it has to be
x
1
x
1
and
y
1
y
1
.
So the symmetric question is: at position 1, what is the most accurate numerical velocity? Well it can't be
y
1
-
y
0
y
1
y
0
or
y
2
-
y
1
y
2
y
1
because those are the velocities to one or the other side of
y
1
y
1
.
The answer has to be
y
2
-
y
0
2
y
2
y
0
2
from symmetry (in particular, time-reversal symmetry where we imagine the direction of motion is reversed).
Thus we have the Lorentz equations (here we have absorbed the
fraction 1/212 into
pp)
x
2
-2
x
1
+
x
0
=
y
2
-
y
0
p
x
2
2
x
1
x
0
y
2
y
0
p
(10)
and
y
2
-2
y
1
+
y
0
=-
x
2
-
x
0
p
y
2
2
y
1
y
0
x
2
x
0
p
(11)
Now we can directly solve for the coordinates. First we have to
chose some initial conditions. For a circular orbit, there is
always a point at which the trajectory is going in, say, the
+x
x
direction, so we can take, without loss of generality,
x
0
=-1
x
0
-1
and
x
1
=+1
x
1
+1
,
and we could assume that the
yy-coordinate was zero at this
point, then
y
0
=0
y
0
0
,
and
y
1
=0
y
1
0
.
But now we notice that the new positions
x
2
x
2
and
y
2
y
2
are mixed together, so we have to do a (very) little algebra to
sort them out,
x
2
-
y
2
p=
c
1
=2
x
1
-
x
0
-
y
0
p
x
2
y
2
p
c
1
2
x
1
x
0
y
0
p
(12)
and
x
2
p+
y
2
=
c
2
=2
y
1
-
y
0
+
x
0
p
x
2
p
y
2
c
2
2
y
1
y
0
x
0
p
(13)
from which we find
x
2
=
c
1
+
c
2
p1+p2
x
2
c
1
c
2
p
1
p
2
(14)
and
y
2
=-
c
1
p-
c
2
1+p2
y
2
c
1
p
c
2
1
p
2
(15)
To show how well the routine works, Figure 1 shows the trajectory in the equatorial plane of a dipole magnetic field. The particle is started at an abitrary distance (here "10") from the dipole center (at the origin) and the gradient drift causes the cyclotron motion (the small circles) to encircle the dipole and continue the drift motion, which clearly overlaps the original, even though the steps are taken so large that the trajectory appears to consist of sharp bends (roughly one radian per step).
One of the lessons one finds repeatedly in numerical schemes is
that very tiny changes make a huge improvements in how well the
scheme works. One way to search for such changes is to insist
that the scheme reproduce the exact solution for some simple
special case.
What was the problem with the "obvious" formulation given at the
beginning? It was the implicit assumption that the
xx and
yy-velocities multiplying
B
0
B
0
were the previous velocities. That is equivalent to using
y
1
-
y
0
y
1
y
0
as the estimator. This had the "simplification" that the two new
velocities could be calculated separately. We didn't have to
solve the simultaneous equations. But that also meant that we
didn't get the right answer!
To fix this, we could also have simply written the Lorentz
equations as
v
x
'
=
v
x
+
v
y
'
+
v
y
p
v
x
'
v
x
v
y
'
v
y
p
(16)
and
v
y
'
=
v
y
-
v
x
'
+
v
x
p
v
y
'
v
y
v
x
'
v
x
p
(17)
where again we absorb the factor of
1/212 (to average the
before and after velocities) into
pp. Let us rearrange this as if we
were solving for the "after" velocities:
v
x
'
-
v
y
'
p=
v
x
+
v
y
p
v
x
'
v
y
'
p
v
x
v
y
p
(18)
and
v
y
'
+
v
x
'
p=
v
y
-
v
x
p
v
y
'
v
x
'
p
v
y
v
x
p
(19)
but now square each side and add. We then get
v
x
'
2+
v
y
'
21+p2=
v
x
2+
v
y
21+p2
v
x
'
2
v
y
'
2
1
p
2
v
x
2
v
y
2
1
p
2
(20)
and we can see that the particle energy is exactly conserved
even if the value of
pp is changed
with each step, as in the case of gradient drift. If we leave
off the
v
'
v
'
terms on the right side of the Lorentz equations, the right side
of the energy equation would lack the
p2
p
2
in the
1+p2
1
p
2
,
so violation of energy conservation is guaranteed! By reducing
the step size,
p2
p
2
declines as the square, and the violation becomes less evident.
The reason Runge-Kutta helps is that the routine effectively
attempts to predict the value of
y
2
y
2
by using the trends in previous values, but there is no need for
such extrapolating when we can solve for the exact value.
If we add in the ignored motion in the
zz-direction, we simply have
ddt
v
z
=0
t
v
z
0
(21)
or
v
z
'
-
v
z
Δt=0
v
z
'
v
z
Δ
t
0
(22)
or
v
z
'
=
v
z
v
z
'
v
z
(23)
which is the remaining conservation equation. Thus in a
uniform magnetic field the particle circles the direction of
the field while drifting freely along that direction (here the
zz-direction). But at any
instant, that is what the particle is doing regardless of how
complicated the magnetic field is: the field points in some
direction and the particle (instantaneously) circles about
that direction while drifting freely parallel to it. So the
general case of motion in 3D is no different. The only
difference is that the coordinates are not necessarily aligned
with
zz pointing along the
field. In principle, one could rotate the coordinates to give
this alignment, calculate the new coordinates as if they were
the 2D ones used above, and rotate them back to the original
orientation. But we do that automatically when we write out
the full set of Lorentz forces:
v
x
'
=
v
x
+
v
y
'
+
v
y
p
z
-
v
z
'
+
v
z
p
y
v
x
'
v
x
v
y
'
v
y
p
z
v
z
'
v
z
p
y
(24)
and so on for the other two components. It is not too easy to factor the equations in an obvious energy-conserving relation like
Equation 20. Instead, we can just rotate the coordinates so that
zz points along
BB and then we have the same "2D" set of equations that we started with. To prove the general case, we notice that the 3D equations become, written in matrix form:
1-pypypz1-px-pypx1v
x'v
y'v
z'=1py-py-pz1pxpy-px1v
xv
yv
z
1
py
py
pz
1
px
py
px
1
v
x'
v
y'
v
z'
1
py
py
pz
1
px
py
px
1
v
x
v
y
v
z
(25)
and we can define an antisymmetric matrix
AA such that
A=0-pypypz0-px-pypx0
A
0
py
py
pz
0
px
py
px
0
in which case the
above equation can be written more compactly as
I+AV'=I-AV
I
A
V'
I
A
V
where
II is the
3×33
3 unit matrix and
VV and
V'V'
are
3×1 3
1 vectors.
The solution for
V'V'
is simply
V'=I+A-1I-AV=BV
V'
I
A
I
A
V
B
V
In other words, BB is the inverse of I+A
I
A
then matrix multiplied by I-A
I
A
.
The resultant program, as originally written, is reproduced below. Here the matrix product to the left of Equation 25 is multiplied out and the velocities have been written as coordinate differences, so
vx=x1-x02dt
vx
x1
x0
2
dt
,
vx'
=x2-x12dt
vx'
x2
x1
2
dt
, etc.
Trajectory in a dipole magnetic field
% DIPOLEORBIT
% Calculates orbits in a dipole field starting
% at unit radius moving with scale velocity V0
% less than unity and pitch angle alpha to the
% equatorial plane
% fcm 19 Nov 1999
% initializing
clear X, Y, Z;
N = 1270; % number of time steps
v0 = 2.0; % initial velocity
alpha = 0.0; % pitch angle
t = 0;
dt = 0.1; % time step
x0 = 1; % intial position
y0 = 0;
z0 = 0;
a = v0*dt*cos(alpha);
b = v0*dt*sin(alpha);
x1 = 1 + a*sin(dt);
y1 = a*(1 - cos(dt) ); % exact solution
z1 = b;
for n = 1:N
r2 = x1*x1 + y1*y1 + z1*z1; % calculate dipole magnetic fields
r = sqrt(r2);
r5 = r*r2*r2;
A = dt*3*z1*x1/(2*r5); % Bx term
B = dt*3*z1*y1/(2*r5);
C = dt*(3*z1*z1 - r2)/(2*r5);
Dx = 2*x1 - x0 - C*y0 + B*z0;
Dy = 2*y1 - y0 - A*z0 + C*x0;
Dz = 2*z1 - z0 - B*x0 + A*y0;
M = [ 1 -C B; C 1 -A; -B A 1 ];
MI = inv(M);
D = [ Dx; Dy; Dz ];
P = MI*D;
x2 = P(1);
y2 = P(2);
z2 = P(3);
X(n) = x0;
Y(n) = y0;
Z(n) = z0;
x0 = x1;
y0 = y1;
z0 = z1;
x1 = x2;
y1 = y2;
z1 = z2;
end
plot(X,Y)
Here with the pitch angle taken to be zero, one gets a
gradient-drift orbit similar to
the figure, but with a faster
drift (10 loops around).
Since
detI+A=detI-A=1+px2+py2+pz2
I
A
I
A
1
px
2
py
2
pz
2
, it follows that
detB=1
B
1, since
detA-1=1detA
A
1
A and
detCD=detCdetD
C
D
C
D
. The remaining thing to prove is that the transpose of B
B (BT
B) is the inverse of BB, namely
BBT=I
B
B
I
. This will prove that
BB is a rotation matrix and therefore simply rotates
VV into
V'
V' without a change in amplitude. Because
AA is antisymmetric,
AT=-A
A
A
and for this special case we can write (using
CDT=DTCT
C
D
D
C
)
BT=I+AI-A-1
B
I
A
I
A
where we have simply reversed the signs of
AA.
Formally,
I-A-1=I+A+AA2+AAA6+...
I
A
I
A
A
A
2
A
A
A
6
...
, and taking the transpose simply changes the sign of AA in every term of the series.
The other special consideration is that
II and
AA commute with one another, which allows us to rearrange BT
B to read
BT=I-A-1I+A
B
I
A
I
A
and the inverses now cancel their corresponding terms to give
BBT=I
B
B
I
. The critical point is that we could write the two matrices on each side of Equation 25 in terms of
II and just one other matrix, AA. Otherwise the argument will not work.
