Numerical code to solve LLG is not preserving norm

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


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.