使用FFT求解热方程的Python实现异常问题排查

核心问题 #

• 热方程符号错误：热传导方程的正确形式是 $\frac{\partial u}{\partial t} = c \frac{\partial^2 u}{\partial x^2}$（其中c为热扩散系数），你代码里写成了 du_dt = -c**2 * dd_u，负号导致方程变为逆热方程——这是一个完全不稳定的方程，数值解会指数级爆炸，直接引发"步长小于数值间距"的错误。
• 初始条件注释矛盾：你注释说明要设置x=2到x=10的脉冲，但定义域L=4，x的最大值仅为4。代码实际是设置x=1到x=3的脉冲，注释与逻辑不匹配（虽非数值爆炸主因，但会导致初始条件不符合预期）。

修正后的代码 #

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

# Constants
c = 1  # 热扩散系数
L = 4    # 定义域长度
N = 100  # 离散点数
dx = L / N
x = np.linspace(0, L, N, endpoint=False)

kappa = 2 * np.pi * fft.fftfreq(N, d=dx)

# 初始条件：x=1到x=3之间设置为1的脉冲
u0 = np.zeros_like(x)
start_idx = int(1 / dx)  # 对应x=1的索引
end_idx = int(3 / dx)    # 对应x=3的索引
u0[start_idx:end_idx] = 1  

def RHSheatspatial(t, u, kappa, c):
    uhat = fft.fft(u)
    dd_uhat = -np.power(kappa, 2) * uhat  # 傅里叶空间的二阶导数
    dd_u = fft.ifft(dd_uhat)
    du_dt = c * np.real(dd_u)  # 修正符号：符合热方程物理意义
    return du_dt

t = np.linspace(0, 10, 1000) 

# 求解IVP
sol = solve_ivp(RHSheatspatial, (0, 10), u0, method='RK45', args=(kappa, c), t_eval=t, rtol=1.0e-6, atol=1.0e-8)
print(sol)

# 绘图
plt.figure(figsize=(12, 5))

# 初始条件
plt.subplot(1, 2, 1)
plt.plot(x, u0, label='初始条件', color='blue')
plt.title('初始条件')
plt.xlabel('x')
plt.ylabel('T')
plt.grid(True)
plt.legend()

# t=10时的解
plt.subplot(1, 2, 2)
plt.plot(x, sol.y[:, -1], label='t=10时的解', color='red')
plt.title('t=10时的温度分布')
plt.xlabel('x')
plt.ylabel('T')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

修正说明 #

1. 符号修正：将du_dt = -c**2 * dd_u改为du_dt = c * np.real(dd_u)，回归热方程的物理本质——热量从高温区域向低温区域扩散，解会逐渐平滑衰减，不会出现数值爆炸。
2. 注释匹配：调整初始条件注释，使其与代码实际逻辑一致，明确脉冲范围为x=1到x=3。
3. 数值稳定性优化：直接对傅里叶逆变换结果取实部，避免冗余操作，保证计算精度。

运行结果 #

修正后，solve_ivp会正常完成积分，解呈现出脉冲逐渐扩散、高度降低的合理热传导行为，不会再出现数值爆炸和步长错误。

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