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

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

文章目录

前言
一、LLG方程和微磁学模拟之间的联系
二、LLG方程的形式
三、微磁学模拟中的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_{eff} })+ \alpha\mathbf{m}\times(\mathbf{m}\times\mathbf{H_{eff}})$

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

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

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

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

在这里插入图片描述
在mumax3的api中我们可以找到如何设置不同的求解器：

四、代码部分

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

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Ddef 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_dtdef 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_kdef 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()