Numerical code to solve LLG is not preserving norm

Question

I am new to this thread. I am trying to do a simple exercise on solving the LLG equation. The equation reads:

$\frac{d\vec{m}}{dt} = -\gamma(\vec{m} \times\vec{H})$.

Given a normalized input state ($m_x$, $m_y$, $m_z$) = (1,1,1) and $\vec{H} = 2\vec{z}$, the expected output should preserve $m_z$ while $m_x$ and $m_y$ should oscillate sinusoidally. The problem is that this code is not preserving the norm. If I force it to preserve norm after each iteration, then the $m_z$ drops.

Any help is much appreciated.

Attaching the code for convenience:

# solving LLG without damping
%pylab inline
gamma = 2.87e10;
tnodes = 100; tmax = 200e-12;
tstep = float(tmax)/tnodes;
trange = linspace(0, tmax, tnodes);
m = zeros((3, 100));
def update_m(tstep, minit, Heff):
    mfinal = minit - tstep * gamma * cross(minit, Heff);
    #mfinal = normalize_m(mfinal)
    return mfinal; 
def normalize_m(m):
    return m/norm(m);
m_0 = (1,1,1);
Heff = (0, 0, 2);
m[:,0] = normalize_m(m_0);

for i in arange(1, len(trange)):
    m[:, i] = update_m(tstep, m[:, i-1], Heff);
    
mx = m[0, :]; my = m[1, :]; mz = m[2, :];
plot(trange/1e-12, mx);
plot(trange/1e-12, my);
plot(trange/1e-12, mz);

Answer 1

Your forward Euler method will not preserve the norm. You need to use a different integrator (such as symplectic integrators).

Answer 2

Just to elaborate a bit on the previous answer, you have a linear ODE $$\dot m = Am$$ where $A$ is an anti-symmetric matrix, i.e. $A^* = -A$. That means that all the eigenvalues of $A$ are purely imaginary. The true solution $m(t) = e^{tA}m(0)$ consists in applying the matrix $e^{tA}$, which you can prove is a unitary matrix and in particular all its eigenvalues are on the unit circle in the complex plane.

If you use a 1st-order explicit scheme with timestep $\delta t$, i.e. $$m(t + \delta t) = m(t) + \delta t\cdot A\cdot m(t) = (I + \delta t\cdot A)m(t),$$ then you're repeatedly applying the matrix $I + \delta t\cdot A$. This is not unitary matrix and so you shouldn't expect it to preserve norms.

Now consider instead the midpoint method $$\frac{m(t + \delta t) - m(t)}{\delta t} = \frac{A\cdot m(t + \delta t) + A\cdot m(t)}{2}$$ which we can rewrite as $$m(t + \delta t) = \left(I - \frac{\delta t}{2}A\right)^{-1}\left(I + \frac{\delta t}{2}A\right)m(t).$$ Using the anti-symmetry of $A$, you can show that this is a unitary matrix. It's equivalent to applying the 1-1 Padé approximation of the exponential function to matrices.

The moral of the story is that the order of accuracy of a numerical method only tells you part of the story. You also need to choose methods that approximate $e^{tA}$ in the part of the complex plane where the spectrum of $A$ lies. If instead $A$ were symmetric and negative-definite, I'd be telling you to use something like the implicit backward scheme instead of the midpoint scheme for exactly this reason.