【微磁学：扒一扒mumax3的内核】LLG方程的多种求解方法

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

最近我很想把现微磁学模拟器的内容拆开看一遍，就像小时候拆收音机一样。目前感兴趣的微磁学模拟器为mumax3，现下免费且GPU加速，微磁学科研普遍使用的软件。会开一些文章记录细节，欢迎同道中人一起交流！！
第一篇给LLG方程。
本篇内容将讲解LLG方程的求解，给出不同的求解方式，并且用一个例子可视化LLG方程求解的过程。

一、LLG方程和微磁学模拟之间的联系

LLG方程是描述磁矩动态行为的基本方程，本质上就是一个偏微分方程，描述磁化矢量的时空变化中的表达方式。
微磁学模拟是一种计算方法，用于模拟磁性材料的微观磁性行为。在微磁学模拟中，材料被划分为许多小单元（通常是有限差分网格），每个单元内的磁矩被假定为均匀的。微磁学模拟通过求解在每个网格点上的LLG方程来模拟磁性材料的磁化过程。成熟的微磁学模拟器比如mumax3还会给出数据处理的结果。
一言以蔽之，微磁学模拟本质上就是解LLG方程，且在很多个格子内求解多个LLG方程。
LLG方程说：“微磁学模拟器是解千千万万个我构成的。”

二、LLG方程的形式

$$\frac{\mathrm{d}\mathbf{m}}{\mathrm{d}t}=-\gamma (\mathbf{m}\times {\mathbf{H}}_{\mathbf{eff}})+\alpha \mathbf{m}\times (\mathbf{m}\times {\mathbf{H}}_{\mathbf{eff}})$$

这个表达式含有两个主要项。

1. 进动项。描述了磁矩向量 $\mathbf{m}$在有效场${\mathbf{H}}_{\mathbf{eff}}$中的进动。$\gamma$是旋磁比，代表磁矩对磁场的响应强度。 因为进动项的存在，磁矩将在有效磁场的作用下进行旋转。
2. 阻尼项。磁矩在进动时，还有一个与上面提到的项的叉积成比例的阻尼作用力。作用为减小磁矩的进动速率，最后使得磁矩趋于稳定在有效磁场的方向。

三、微磁学模拟中的LLG方程求解部分

对于LLG方程的求解，微磁学模拟软件往往会给出不同的求解算法以求得更快更好的解，比如在mumax3设计的文章《The design and verification of Mumax3》中有提到：

在mumax3的api中我们可以找到如何设置不同的求解器：

四、代码部分

因为mumax3的源代码部分为GO语言，而且不同模块之间的联合看起来很复杂，这部分我得慢慢剥。为了理解LLG方程的不同算法的求解，我这里先写了一些代码去理解不同算法的求解方式，可作为参考。求解的算法对应的是上面引用的Mumax3中的solver，可作为参考。
关于求解偏微分方程的基本思路，可以看我这篇：[有限差分法的求解思路](https://blog.youkuaiyun.com/CRUSH8496052/article/details/138095763?

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

def euler_step(M, H_eff, alpha, gamma, dt):
    """执行单个Euler步骤。"""
    dM_dt = -gamma * np.cross(M, H_eff) + alpha * np.cross(M, np.cross(M, H_eff))
    return M + dt * dM_dt

def rk4_step(M, H_eff, alpha, gamma, dt):
    """执行单个RK4步骤。"""
    k1 = -gamma * np.cross(M, H_eff) + alpha * np.cross(M, np.cross(M, H_eff))
    k2 = -gamma * np.cross(M + 0.5 * dt * k1, H_eff) + alpha * np.cross(M + 0.5 * dt * k1, np.cross(M + 0.5 * dt * k1, H_eff))
    k3 = -gamma * np.cross(M + 0.5 * dt * k2, H_eff) + alpha * np.cross(M + 0.5 * dt * k2, np.cross(M + 0.5 * dt * k2, H_eff))
    k4 = -gamma * np.cross(M + dt * k3, H_eff) + alpha * np.cross(M + dt * k3, np.cross(M + dt * k3, H_eff))
    return M + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)

def heun_step(M, H_eff, alpha, gamma, dt):
    """执行单个Heun步骤。"""
    # 初始斜率
    k1 = -gamma * np.cross(M, H_eff) + alpha * np.cross(M, np.cross(M, H_eff))
    # 使用Euler方法预测的终点
    M_euler = M + dt * k1
    # 校正斜率
    k2 = -gamma * np.cross(M_euler, H_eff) + alpha * np.cross(M_euler, np.cross(M_euler, H_eff))
    # 校正后的平均斜率
    average_k = 0.5 * (k1 + k2)
    return M + dt * average_k

def dormand_prince_step(M, H_eff, alpha, gamma, dt):
    """执行单个Dormand-Prince步骤。"""
    k1 = -gamma * np.cross(M, H_eff) + alpha * np.cross(M, np.cross(M, H_eff))
    k2 = -gamma * np.cross(M + 0.2 * dt * k1, H_eff) + alpha * np.cross(M + 0.2 * dt * k1, np.cross(M + 0.2 * dt * k1, H_eff))
    # 需要更多阶段以实现完整
    # 简化以展示结构
    return M + dt * (0.2 * k1 + 0.8 * k2)

def bogacki_shampine_step(M, H_eff, alpha, gamma, dt):
    """执行单个Bogacki-Shampine步骤。"""
    k1 = -gamma * np.cross(M, H_eff) + alpha * np.cross(M, np.cross(M, H_eff))
    k2 = -gamma * np.cross(M + 0.5 * dt * k1, H_eff) + alpha * np.cross(M + 0.5 * dt * k1, np.cross(M + 0.5 * dt * k1, H_eff))
    k3 = -gamma * np.cross(M + 0.75 * dt * k2, H_eff) + alpha * np.cross(M + 0.75 * dt * k2, np.cross(M + 0.75 * dt * k2, H_eff))
    return M + (dt * (2 * k1 + 3 * k2 + 4 * k3) / 9)


def solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, method):
    """
    使用指定的方法数值解决LLG方程。

    参数:
        M0: numpy数组，初始磁化向量
        H_eff: numpy数组，有效磁场向量
        alpha: 浮点数，自旋耗散参数
        gamma: 浮点数，旋磁比
        dt: 浮点数，时间步长
        T:浮点数，总时间
        method: 函数，要使用的数值方法（euler_step 或 rk4_step）

    返回:
        numpy数组，每个时间步骤的磁化向量
    """
    num_steps = int(T / dt)  # 时间步数
    M = np.zeros((num_steps, len(M0)))  # 存储每个时间步的磁化向量
    M[0] = M0  # 设置初始磁化

    for i in range(1, num_steps):
        M[i] = method(M[i-1], H_eff, alpha, gamma, dt)

    return M

# 示例使用
M0 = np.array([1, 0, 0])  # 初始磁化向量
H_eff = np.array([0, 0, 1])  # 有效磁场向量
alpha = 0.05  # 自旋耗散参数
gamma = 1.0  # 旋磁比
dt = 0.0005  # 时间步长
T = 5.0  # 总时间

euler_solution = solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, euler_step)
rk4_solution = solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, rk4_step)
heun_solution = solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, heun_step)
dormand_prince_solution = solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, dormand_prince_step)
bogacki_shampine_solution = solve_llg_equation(M0, H_eff, alpha, gamma, dt, T, bogacki_shampine_step)

# 绘制3D轨迹
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 绘制不同的解法轨迹
ax.plot(euler_solution[:, 0], euler_solution[:, 1], euler_solution[:, 2], label='RK1', color='blue')
ax.plot(heun_solution[:, 0], heun_solution[:, 1], heun_solution[:, 2], label='RK12', color='green')
ax.plot(bogacki_shampine_solution[:, 0], bogacki_shampine_solution[:, 1], bogacki_shampine_solution[:, 2], label='RK32', linestyle='--', color='red')
ax.plot(rk4_solution[:, 0], rk4_solution[:, 1], rk4_solution[:, 2], label='RK4', linestyle='-.', color='black')
ax.plot(dormand_prince_solution[:, 0], dormand_prince_solution[:, 1], dormand_prince_solution[:, 2], label='RK45', linestyle=':', color='purple')

# 设置标签、图例和网格
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()
ax.grid(True)

# 设置视角
ax.view_init(elev=30, azim=120)

# 显示图像
plt.show()
plt.show()

运行结果：

可以看到不同的算法求得的结果都是正确的。

五、总结

今天的总结内容主要是关于微磁学仿真器的核心部分LLG方程的解法，自己尝试了不同的求解算法，接下来需要从本身的设计构造去看。尝试扒一扒官网的源代码。