基于Python求解热方程并复现指定温度波动曲线

我正在求解热方程，目标是复现指定的温度波动曲线图。目前已有有限差分法的Python实现代码，同时通过热方程解析解代码绘制出了目标曲线，现在需要调整有限差分法代码，使其求解结果与解析解一致。

原有限差分法代码 #

import numpy as np
import matplotlib.pyplot as plt

L = 0.04     # 土壤深度
Te = 28    # 外部温度
Tn = 36.6      # x=L处的温度
D = 0.75    # 扩散系数 (m²/s)
tau = 232920  # 实验总时长
dt = 600   
dx = 0.001  
Tm=38   # 温度波动幅值
T=24*3600 # 周期
w=(2*np.pi)/T

Nt = int(tau/dt)+1  
Nx = int(L/dx)+1    

x = [j*dt for j in range(Nx)]       
T = np.zeros((Nt,Nx))               
for j in range(Nx):
    T[0,j] = Te    # 初始条件：t=0时所有位置温度为Te
T[0,0] = Te         # x=0处初始温度
for i in range(1,Nt): 
    T[i,0] = Te+Tm*np.cos(w*i)     # x=0处的边界条件
    T[i,-1] = Tn    # x=L处的边界条件


for i in range(Nt-1):   
    for j in range(1,Nx-1):    
        T[i+1,j]=T[i,j]+D*dt/(dx*dx)*(T[i,j+1]+T[i,j-1]-2*T[i,j])

t=np.arange(0, 232920, 10*60)
plt.plot(t, T[i,:], label='temperature fluctuation')
plt.xlabel('time (s)')
plt.ylabel('temperature (c)')
plt.legend(loc='lower right')
plt.show()

解析解代码 #

import numpy as np
import matplotlib.pyplot as plt
from math import *


t0=28
tm=38
x= 4*10e-2  # 对应有限差分法中的L=0.04m
T=24*3600
w=(2*pi)/T
D=0.75
delta=sqrt(2*D/w)
t=np.arange(0, 232920, 10*60)
temp = t0 + tm*np.exp(-x/delta)*np.cos(w*t-(x/delta))

plt.plot(t, temp, label='temperature fluctuation')
plt.xlabel('time (s)')
plt.ylabel('temperature (c)')
plt.legend(loc='lower right')
plt.show()

问题分析与修正后的有限差分法代码 #

原有限差分法代码存在多个问题导致结果无法匹配解析解：

• 变量名冲突：周期变量T与温度数组T重名，导致后续计算错误
• 空间坐标计算错误：x = [j*dt for j in range(Nx)]应该改为基于空间步长dx计算位置
• 边界条件时间项错误：w*i应为w*i*dt，因为i是时间步索引，实际时间是i*dt
• 数值稳定性问题：显式格式的CFL条件D*dt/(dx²) ≤ 0.5不满足，原参数下该值为450000，远大于阈值，导致数值发散，需改用隐式格式
• 绘图数据错误：原代码绘制的是最后一个时间步的空间温度分布，而解析解是固定位置（x=0.04m）的时间序列，应取T[:, -1]的时间序列

以下是修正后的代码，采用Crank-Nicolson隐式格式保证稳定性：

import numpy as np
import matplotlib.pyplot as plt

L = 0.04     # 土壤深度
Te = 28      # 外部温度
Tn = 36.6    # x=L处的温度
D = 0.75     # 扩散系数 (m²/s)
tau = 232920 # 实验总时长
dt = 600     # 时间步长
dx = 0.001   # 空间步长
Tm = 38      # 温度波动幅值
period = 24*3600 # 周期
w = (2 * np.pi) / period

Nt = int(tau / dt) + 1  
Nx = int(L / dx) + 1    

# 生成正确的空间坐标数组
x = np.linspace(0, L, Nx)       
# 温度数组初始化
T = np.zeros((Nt, Nx))               
# 初始条件：t=0时所有位置温度为Te
T[0, :] = Te    

# 构建Crank-Nicolson格式的系数矩阵
r = D * dt / (2 * dx**2)
A = np.diag(1 + 2*r * np.ones(Nx-2)) + np.diag(-r * np.ones(Nx-3), k=1) + np.diag(-r * np.ones(Nx-3), k=-1)

for i in range(Nt-1):
    # 边界条件
    T[i+1, 0] = Te + Tm * np.cos(w * i * dt)  # x=0处的时变边界
    T[i+1, -1] = Tn                           # x=L处的固定边界
    
    # 构建右端向量
    b = T[i, 1:-1] * (1 - 2*r) + r * (T[i, 2:] + T[i, :-2])
    # 处理边界条件对右端向量的影响
    b[0] += r * T[i+1, 0]
    b[-1] += r * T[i+1, -1]
    
    # 求解线性方程组
    T[i+1, 1:-1] = np.linalg.solve(A, b)

# 生成与解析解对齐的时间序列
t = np.arange(0, tau, 10*60)
# 提取x=L处的温度时间序列（对应解析解的位置）
fd_temp = T[::10, -1]  # 每10个时间步取一个点，和解析解采样间隔一致

# 绘制对比图
plt.plot(t, fd_temp, label='有限差分法结果')
# 绘制解析解
delta = np.sqrt(2*D/w)
analytical_temp = Te + Tm*np.exp(-L/delta)*np.cos(w*t - (L/delta))
plt.plot(t, analytical_temp, label='解析解', linestyle='--')

plt.xlabel('时间 (s)')
plt.ylabel('温度 (℃)')
plt.legend(loc='lower right')
plt.show()

修改说明 #

1. 变量名修正：将周期变量改为period，避免与温度数组冲突
2. 空间坐标修正：用np.linspace生成正确的空间位置数组
3. 边界条件时间项修正：使用w*i*dt计算实际时间对应的余弦项
4. 稳定性处理：采用Crank-Nicolson隐式格式，该格式无条件稳定，不受CFL条件限制
5. 绘图数据修正：提取x=L处的温度时间序列，与解析解的位置对应，并保持相同的时间采样间隔
6. 添加对比绘图：同时绘制有限差分法结果和解析解，方便验证一致性

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