Python复现Mathematica磁化计算代码除零报错、无输出问题排查

核心错误排查 #

你的Python代码无法正常运行是由5个关键逻辑错误导致的：

• 迭代初值作用域错误：Mathematica中每个温度t1、每个磁场h11的自洽迭代都从初始值i1=7、i2=1.92开始，你将i1/i2定义在所有循环外部，上一轮迭代收敛的结果会带入下一轮计算，极易迭代到0值触发除零，或进入震荡死循环导致长时间无输出。
• 布里渊函数公式错误：Mathematica中计算x1/x2、m1/m2时用的是固定角动量量子数7、1.92作为系数，你错误替换成了迭代变量i1/i2，完全偏离了物理公式。
• 收敛判断逻辑错误：Mathematica直接判断新计算的磁化强度与上一轮迭代值的差值小于阈值，你额外给m1/m2套了绝对值，收敛条件不匹配会导致无法正常跳出循环。
• 缺少死循环保护：while True没有设置最大迭代次数，一旦迭代震荡会永久卡住。
• 数据存储与后续逻辑缺失：存储数据时没有保存温度t1的取值，也没有实现Mathematica中的插值求导、绘图逻辑，无法复现最终结果。

修正后完整可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

# 物理常量与参数定义，和Mathematica完全一致
u = 9.27 * 10.0**(-21)
k = 1.38 * 10.0**(-16)
j1 = 3.5
j2 = 1.2
nrr = 9740
nrm = -131180
nmm = 1072770
na = 6 * 10.0**(23)
x = 0.5

# 自定义coth函数，加极小值截断避免除零报错
def coth(x):
    x = np.where(np.abs(x) < 1e-12, 1e-12, x)
    return 1 / np.tanh(x)

# 遍历Mathematica指定的两个磁场值
for h11 in [10000, 20000]:
    mgData = []
    mg1Data = []
    mg2Data = []
    # 遍历温度点1~600K
    for t1 in range(1, 601):
        # 每个温度点必须重置迭代初值
        i1 = 7.0
        i2 = 1.92
        converged = False
        # 最大迭代10000次，防止死循环
        for iter_step in range(10000):
            h1 = h11 + (1 - x) * nrr * i1 + 2 * nrm * i2
            h2 = h11 + (1 - x) * nrm * i1 + 2 * nmm * i2
            # 注意：系数为固定值7、1.92，不是迭代变量i1/i2
            x1 = (7 * u * h1) / (k * t1)
            x2 = (1.92 * u * h2) / (k * t1)
            x3 = 2 * j1
            x4 = 2 * j2
            q1 = (x3 + 1) / x3
            q2 = (x4 + 1) / x4
            # 布里渊函数计算磁化强度
            m1 = 7 * (q1 * coth(q1 * x1) - (1/x3) * coth(x1/x3))
            m2 = 1.92 * (q2 * coth(q2 * x2) - (1/x4) * coth(x2/x4))
            m = np.abs((1 - x) * m1 + 2 * m2)
            # 收敛判断和Mathematica逻辑完全一致
            if np.abs(m1 - i1) < 1e-7 and np.abs(m2 - i2) < 1e-7:
                mgData.append([t1, m])
                mg1Data.append([t1, m1])
                mg2Data.append([t1, m2])
                converged = True
                break
            # 更新迭代值
            i1, i2 = m1, m2
        if not converged:
            print(f"警告：H={h11}Oe, T={t1}K迭代未收敛")
    # 转numpy数组方便处理
    mgData = np.array(mgData)
    mg1Data = np.array(mg1Data)
    mg2Data = np.array(mg2Data)
    # 插值求dM/dT，和Mathematica Interpolation逻辑匹配
    m_interp = interp1d(mgData[:,0], mgData[:,1], kind='cubic', fill_value="extrapolate")
    t_seq = np.arange(1, 601, 1)
    dm_dt = np.gradient(m_interp(t_seq), t_seq)
    # 绘制曲线，和Mathematica ListPlot结果对应
    plt.plot(mgData[:,0], mgData[:,1], label=f'总磁化强度 H={h11}Oe')
    plt.plot(mg1Data[:,0], mg1Data[:,1], '--', label=f'次晶格磁化强度M1 H={h11}Oe')
    plt.plot(mg2Data[:,0], mg2Data[:,1], ':', label=f'次晶格磁化强度M2 H={h11}Oe')

# 绘图格式设置
plt.xlabel('温度 K')
plt.ylabel(r'M $\mu_B$/f.u')
plt.legend()
plt.gcf().set_frameon(True)
plt.show()

补充说明 #

• 代码中使用numpy向量化计算替代mpmath的高精度计算，运行速度提升数十倍，计算精度和Mathematica默认双精度计算完全一致。
• 磁熵计算部分需要0到目标磁场之间足够多间隔的磁场点M-T数据，才能完成dM/dT对磁场的数值积分，如果你补充全磁场点数据，直接用scipy.integrate.simpson做数值积分即可复现Mathematica中ds的计算结果。

内容的提问来源于stack exchange，提问作者RaaaX