热输运解析解精度不足：Python计算优化方案问询

问题背景 #

我在用Python计算热输运问题的Barends解析解（公式4）以验证数值模型，通过scipy.integrate.quad计算积分后，与单调递减的前缀项相乘得到ΔT值。但结果出现不光滑的问题，且时间t越大，该现象越明显，怀疑是数值精度导致的，想知道如何提升乘法运算及整体计算的精度。

附可复现代码：

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

# 示例参数
alpha = 1.0
k = 0.5
x = 0.1
t_list = np.logspace(0, 3, 100)  # 覆盖1到1000的时间范围

def integrand(s, x, t, alpha, k):
    # Barends公式4对应的被积函数（示例形式）
    return np.exp(-s**2/(4*alpha*t)) / (s**2 + k**2)

def delta_T(t, x, alpha, k):
    # 单调递减的前缀项
    prefix = (1/(2*np.pi)) * np.exp(-k*x)
    # 计算积分
    integral, _ = quad(integrand, 0, np.inf, args=(x, t, alpha, k))
    return prefix * integral

# 计算所有时间点的ΔT
delta_T_results = [delta_T(t, x, alpha, k) for t in t_list]

# 绘图展示结果
plt.loglog(t_list, delta_T_results)
plt.xlabel('Time t')
plt.ylabel('ΔT')
plt.title('Barends解析解计算结果')
plt.grid(True)
plt.show()

解决方案 #

1. 提升quad积分的精度 #

结果不光滑的核心原因大概率是积分计算的误差，而非乘法本身。quad默认的精度参数（epsabs=1e-8、epsrel=1e-8）可能不足以满足大t场景下的计算需求。可以调小这两个参数，强制积分器使用更高精度：

integral, _ = quad(integrand, 0, np.inf, args=(x, t, alpha, k), epsabs=1e-12, epsrel=1e-12)

更小的精度阈值会让quad迭代更多次，捕捉到被积函数的细微变化，减少积分误差。

2. 使用更高精度的浮点数类型 #

Python默认的float是64位双精度，对于大数乘小数的场景（大t时积分结果很大，前缀项很小，相乘时有效数字损失），可以切换到128位扩展精度浮点数：

import numpy as np

# 改用float128定义所有参数和变量
alpha = np.float128(1.0)
k = np.float128(0.5)
x = np.float128(0.1)
t_list = np.logspace(0, 3, 100, dtype=np.float128)

def integrand(s, x, t, alpha, k):
    return np.exp(-s**2/(4*alpha*t)) / (s**2 + k**2)

def delta_T(t, x, alpha, k):
    prefix = (1/(2*np.pi)) * np.exp(-k*x)
    integral, _ = quad(integrand, 0, np.inf, args=(x, t, alpha, k), epsabs=np.float128(1e-12), epsrel=np.float128(1e-12))
    return prefix * integral

# 计算后转换回float64用于绘图
delta_T_results = np.array([delta_T(t, x, alpha, k) for t in t_list], dtype=np.float64)
t_plot = t_list.astype(np.float64)

plt.loglog(t_plot, delta_T_results)
# ... 绘图代码不变

注意：float128的计算速度会比float64慢，且部分库对其支持有限，需根据实际需求权衡。

3. 改写表达式避免数值抵消 #

大t时，积分结果是大数，前缀项是小数，两者相乘会出现有效数字抵消的问题。可以将前缀项融入被积函数，直接计算积分，避免单独的大数乘小数操作：

def delta_T(t, x, alpha, k):
    def integrand_with_prefix(s, x, t, alpha, k):
        # 将前缀项与被积函数合并
        prefix = (1/(2*np.pi)) * np.exp(-k*x)
        return prefix * np.exp(-s**2/(4*alpha*t)) / (s**2 + k**2)
    integral, _ = quad(integrand_with_prefix, 0, np.inf, args=(x, t, alpha, k), epsabs=1e-12, epsrel=1e-12)
    return integral

这种方式让积分器直接计算最终的小数值，减少中间步骤的精度损失。

4. 拆分积分区间优化计算 #

当t很大时，被积函数的主要贡献集中在s很小的区间（因为指数项exp(-s²/(4αt))会快速衰减）。可以手动拆分积分区间，对核心区间用更高精度计算，其余区间做近似：

def delta_T(t, x, alpha, k):
    prefix = (1/(2*np.pi)) * np.exp(-k*x)
    # 拆分区间：0到s_cutoff，s_cutoff到无穷大
    s_cutoff = 10 * np.sqrt(alpha*t)  # 指数项衰减到可忽略的区间边界
    integral1, _ = quad(integrand, 0, s_cutoff, args=(x, t, alpha, k), epsabs=1e-12, epsrel=1e-12)
    # 对s_cutoff到无穷大的区间，用渐近近似计算（示例）
    integral2 = np.sqrt(np.pi * alpha * t) * np.exp(-k**2 * alpha * t) / (2 * k)
    return prefix * (integral1 + integral2)

这种方法利用被积函数的特性减少计算量，同时提升核心区间的计算精度。

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