I am trying to simulate LLG equation without damping. The equation is

$$\frac{d\vec{m}}{dt} = \vec{m}\times\vec{H}$$

I am solving in spherical coordinates as LLG equation is known to have problems in rectangular coordinates.

My output should look like sine wave for all three components of the $\vec{m}$. Instead I am getting an output like this,

enter image description here

My python code is

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D

# LLG equation without damping (alpha = 0)
def LLG_equation(t, M, gamma, Heff):
    Mr, Mtheta, Mphi = M
    Hr, Htheta, Hphi = Heff

    # dMr_dt = -gamma * (Mtheta * Hphi - Mphi * Htheta)
    dMr_dt = 0 #Because m is always a unit length vector
    dMtheta_dt = -gamma * (Mphi * Hr - Mr * Hphi)
    dMphi_dt = -gamma * (Mr * Htheta - Mtheta * Hr)

    return [dMr_dt, dMtheta_dt, dMphi_dt]

# Parameters
gamma = 1.0  # Gyromagnetic ratio

M0_rect = np.array([1.0, 1.0, 1.0])
r = np.linalg.norm(M0_rect)
M0_rect /= r
r = np.linalg.norm(M0_rect)
theta = np.arccos(M0_rect[2] / r)
phi = np.arctan2(M0_rect[1], M0_rect[0])
M0 = np.array([r, theta, phi])

Heff = np.array([0.0, 0.0, 1.0])
r = np.linalg.norm(Heff)
theta = np.arccos(Heff[2] / r)
phi = np.arctan2(Heff[1], Heff[0])
Heff = np.array([r, theta, phi])

# Time values
t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the LLG equation numerically using scipy's solve_ivp
solution = solve_ivp(
    args=(gamma, Heff),

# Extract the components
Mr_solution, Mtheta_solution, Mphi_solution = solution.y

# Calculate the x, y, z components of magnetization
mx_solution = Mr_solution * np.sin(Mtheta_solution) * np.cos(Mphi_solution)
my_solution = Mr_solution * np.sin(Mtheta_solution) * np.sin(Mphi_solution)
mz_solution = Mr_solution * np.cos(Mtheta_solution)

# Create a 2D plot for each component
plt.figure(figsize=(10, 6))

plt.plot(t_eval, mx_solution, label='mx')
plt.plot(t_eval, my_solution, label='my')
plt.plot(t_eval, mz_solution, label='mz')

plt.ylabel('Component Values')
plt.title('Time Evolution of Magnetization Components')

Is my code flawed? Please help. EDIT : Vector transformations $$\begin{align*} \mathbf{M} &= M_r \mathbf{e}_r + M_\theta \mathbf{e}_\theta + M_\phi \mathbf{e}_\phi \\ \mathbf{H} &= H_r \mathbf{e}_r + H_\theta \mathbf{e}_\theta + H_\phi \mathbf{e}_\phi \\ \mathbf{M} \times \mathbf{H} &= \left( M_\theta H_\phi - M_\phi H_\theta \right) \mathbf{e}_r + \left( M_\phi H_r - M_r H_\phi \right) \mathbf{e}_\theta + \left( M_r H_\theta - M_\theta H_r \right) \mathbf{e}_\phi \end{align*}$$

Here since $\vec{m}$ has magnitude 1, and therefore has zero derivative, $M_\theta H_\phi - M_\phi H_\theta$ is zero.

  • 1
    $\begingroup$ You should now notice that you expressed $H$ in the local basis at the point $M$ and not in the global polar coordinates. It is also easy to verify that $M_\theta=M_\phi=0$ and $M_r=r$. $\endgroup$ Oct 10, 2023 at 9:29
  • $\begingroup$ I am really embarrassed to ask, but I don't see that I have done that. I used global basis vectors didn't I? $\endgroup$
    – User
    Oct 10, 2023 at 10:39
  • $\begingroup$ $e_r,e_\theta,e_\phi$ are aligned (normal and two times tangent) to the tangent plane to the sphere of constant radius $r=r_M$ at $M$ and the slope of the polar coordinate grid there. As this basis depends on $M$, I call it local. The coordinates $r,\theta,\phi$ belong to the global, fixed spherical or polar coordinate system, based on the global and fixed Cartesian basis $e_1,e_2,e_3$. $\endgroup$ Oct 10, 2023 at 10:51
  • $\begingroup$ Oh, yes, YES!! It all makes sense. Thank you so much. I see that my code is wrong, because $\vec{v} = x \bold{i} + y \bold{j} + z \bold{k}$ is just $\vec{\rho}$ and other components are zero. Just to make my dumbass brain more clear, my code is wrong where it converted the vectors isn't it? I just converted the coordinate points and not the vector. I used to hate this section when I studied electromagnetic theory. $\endgroup$
    – User
    Oct 10, 2023 at 11:01
  • 1
    $\begingroup$ Now this seems to be plain wrong. The relation between Cartesian and spherical coordinate systems is not linear and also not pseudo-linear (that is, of the form $A(\vec x)\vec x$). At least not in a natural fashion, and if you force it, you make the calculations more complicated. $\endgroup$ Oct 10, 2023 at 14:30

1 Answer 1


The cross product with $H$ acts like a complex unit in the plane that has $H$ as normal, especially if you select $H$ to be a unit vector. Components parallel to $H$ remain unchanged. What you should thus see is a rotation in the $xy$ plane of unit speed.

I have doubt that in the equations in polar coordinates you get the naked angles in the product. The original system is invariant when the angles are shifted by $2\pi$, this is not the case for your system.

The current point is $m=r\cdot e(\theta,\phi)$. The local rectangular frame consists of $e_r=e(\theta,\phi)$, $e_\theta=e(\theta+\frac\pi2,\phi)$ and $e_\phi=e(\theta,\phi+\frac\pi2)$. Thus $$ \dot m=e_r\dot r+r e_\theta\dot\theta+r e_\phi\dot\phi=R(\theta,\phi)\pmatrix{\dot r\\r\dot\theta\\r\dot\phi},~~~~R=R(\theta,\phi)=(e_r,e_\theta,e_\phi). $$ $R$ is the rotation matrix taking the local basis vectors as columns. It thus transforms local coordinates into global Cartesian coordinates. Its inverse $R^T$ then inversely transforms global coordinates into local coordinates. Using the transformation laws of the cross product transforms the ODE system to $$ \pmatrix{\dot r\\r\dot\theta\\r\dot\phi}=(R^Tm)\times (R^TH)=(re_1)\times (R^TH) =r\pmatrix{0\\-(R^TH)_3\\(R^TH)_2} $$ $e_1,e_2,e_3$ being the Cartesian basis vectors. As $$ (R^TH)_2=\langle R^TH,e_2\rangle =\langle H,Re_2\rangle = \langle H,e_\theta\rangle $$ and similarly $(R^TH)_3=\langle H,e_\phi\rangle$, a simple computation algorithm follows, but the polar coordinates of $H$ are not very helpful here.

  • 1
    $\begingroup$ I thank you for your input. But my Linear Algebra skills are less than mediocre. Is there anyway you could explain this in a bit more? $\endgroup$
    – User
    Oct 10, 2023 at 9:22
  • 2
    $\begingroup$ From what point on? I used that $\cos(x+k\frac\pi2)$, $k=0,1,2,...$ gives the sequence of derivatives $\cos x,\,-\sin x,\,-\cos x,\,\sin x,\,...$ and similar for the sine to express the derivatives of $e(\theta,\phi)$. Everything else are variants of transforming between the local and the global basis, what is intuitiv varies for everyone. Your addition is fully compatible with what I wrote. Note that $H_\phi=⟨H,e_ϕ⟩=(R^TH)_3$ etc. So where do you want to start me explaining? $\endgroup$ Oct 10, 2023 at 9:36

Your Answer

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

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