# Is my differential equation solving code wrong?

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,

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(
LLG_equation,
t_span,
M0,
t_eval=t_eval,
args=(gamma, Heff),
method='RK45',
rtol=1e-6,
atol=1e-6
)

# 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.xlabel('Time')
plt.ylabel('Component Values')
plt.title('Time Evolution of Magnetization Components')
plt.legend()
plt.grid(True)
plt.show()

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.

• 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$. Oct 10, 2023 at 9:29
• I am really embarrassed to ask, but I don't see that I have done that. I used global basis vectors didn't I?
– User
Oct 10, 2023 at 10:39
• $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$. Oct 10, 2023 at 10:51
• 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.
– User
Oct 10, 2023 at 11:01
• 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. Oct 10, 2023 at 14:30

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.
• 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? Oct 10, 2023 at 9:36