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

2 Answers 2


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


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.