无阻尼LLG方程模拟代码输出异常，请求排查问题

阿华AIGC实验室

2026-4-21

嗨，我仔细看了你的代码，发现几个关键问题导致输出不符合预期，咱们一步步梳理清楚：

核心问题1：坐标系混淆，有效场（Heff）处理错误

你把直角坐标的[0,0,1]（沿z轴）转换成球坐标后传入函数，但这里犯了一个致命错误：球坐标的单位矢量是随磁化强度$\vec{m}$的位置变化的局部坐标系，而有效场Heff是固定在全局直角坐标系的z轴上的，不能直接转换为全局球坐标后传入函数。这导致你实际计算时用的Heff是沿$\vec{m}$的径向，而非原本的z轴，完全改变了物理场景。

正确的做法是：Heff始终用直角坐标传入，在函数内部将球坐标的$\vec{m}$转换为直角坐标后再计算叉乘。

核心问题2：球坐标下LLG方程的导数推导错误

你直接套用直角坐标叉乘的分量公式到球坐标，这是不对的——球坐标的矢量时间导数和叉乘规则和直角坐标完全不同。更简单不易出错的方式是：

把球坐标的$\vec{m}$转换为直角坐标
计算直角坐标下的叉乘$\vec{m} \times \vec{H}_{eff}$
再将直角坐标的$\dot{\vec{m}}$转换回球坐标的导数（$\dot{r}, \dot{\theta}, \dot{\phi}$）

修正后的代码

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

# 无阻尼LLG方程：球坐标输入，直角坐标Heff
def LLG_equation(t, M, gamma, Heff):
    r, theta, phi = M
    # 球坐标转直角坐标的磁化强度（单位矢量）
    mx = r * np.sin(theta) * np.cos(phi)
    my = r * np.sin(theta) * np.sin(phi)
    mz = r * np.cos(theta)
    m = np.array([mx, my, mz])
    
    # 标准无阻尼LLG方程：d m/dt = -gamma * (m × Heff)
    cross_product = np.cross(m, Heff)
    dm_dt = -gamma * cross_product  # 若要和你原方程形式一致，可去掉负号
    
    # 直角坐标导数转球坐标导数
    dr_dt = np.dot(m, dm_dt) / r  # 单位矢量时dr_dt恒为0
    # 计算dθ/dt：从mz = r cosθ求导推导而来
    if np.abs(np.sin(theta)) < 1e-10:
        dtheta_dt = 0.0  # 避免θ接近0/π时除以0
    else:
        dtheta_dt = (dr_dt * np.cos(theta) - dm_dt[2]) / (r * np.sin(theta))
    # 计算dφ/dt：从mx、my的导数联立推导而来
    numerator = mx * dm_dt[1] - my * dm_dt[0]
    denominator = r**2 * np.sin(theta)**2
    dphi_dt = numerator / denominator if np.abs(denominator) > 1e-10 else 0.0
    
    return [dr_dt, dtheta_dt, dphi_dt]

# 参数设置
gamma = 1.0  # 旋磁比
# 初始磁化强度：直角坐标归一化后转球坐标
M0_rect = np.array([1.0, 1.0, 1.0])
M0_rect /= np.linalg.norm(M0_rect)
r0 = 1.0
theta0 = np.arccos(M0_rect[2] / r0)
phi0 = np.arctan2(M0_rect[1], M0_rect[0])
M0 = np.array([r0, theta0, phi0])

Heff_rect = np.array([0.0, 0.0, 1.0])  # 沿z轴的有效场，直接用直角坐标

# 时间配置
t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# 求解方程
solution = solve_ivp(
    LLG_equation,
    t_span,
    M0,
    t_eval=t_eval,
    args=(gamma, Heff_rect),
    method='RK45',
    rtol=1e-6,
    atol=1e-6
)

# 转换回直角坐标并绘图
Mr_sol, Mtheta_sol, Mphi_sol = solution.y
mx_sol = Mr_sol * np.sin(Mtheta_sol) * np.cos(Mphi_sol)
my_sol = Mr_sol * np.sin(Mtheta_sol) * np.sin(Mphi_sol)
mz_sol = Mr_sol * np.cos(Mtheta_sol)

plt.figure(figsize=(10, 6))
plt.plot(t_eval, mx_sol, label='mx')
plt.plot(t_eval, my_sol, label='my')
plt.plot(t_eval, mz_sol, label='mz')
plt.xlabel('Time')
plt.ylabel('Magnetization Component')
plt.title('Time Evolution of Magnetization (Undamped LLG)')
plt.legend()
plt.grid(True)
plt.show()