一维热扩散解析解与数值解匹配问题排查及实现指导

问题排查与修正方案 #

1. 数值解边界条件错误 #

你标注右壁温度为100°C，但代码中错误将左壁（T[:,0]）设为100，且右壁（数组索引N-1）未设置边界条件，始终保持初始0值。结合你提供的解析解模型，正确边界条件应为：

• 右壁（x=L，对应N-1索引）固定100°C：T[:, N-1] = 100
• 左壁（x=0，对应0索引）为绝热边界（温度导数为0），有限差分中用T[k+1, 0] = T[k, 1]（等价于T_{-1}=T_1，模拟绝热效果）

2. 解析解公式实现错误 #

代码中正弦项的计算逻辑完全偏离公式：
原公式为sin(nπx/(2L))，你写成了(n * np.pi * x) / 2 * L，这等价于(nπxL)/2，应修正为(n * np.pi * x) / (2 * L)。

3. 数值解初始条件缺失 #

代码中初始化T为全0，但未显式声明初始时刻温度分布，建议补充T[0, :] = 0明确初始温度为0°C，与解析解的初始条件匹配。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

N = 20  # Number of points
L = 1   # Length of the rod (m)
total_time_steps = 300  # Number of time steps
a = 0.001  # Diffusion coefficient
dx = L / N
dt = 1  # Time step size

r = (a * dt) / dx**2
print(f'Convergence factor is: {r:.4f} < 0.5' if r < 0.5 else f'Convergence factor {r:.4f} >=0.5, numerical solution may diverge')

# Initialize arrays
T_num = np.zeros((total_time_steps, N))
T_analytical = np.zeros((total_time_steps, N))
x = np.linspace(0, L, N)

# Set boundary and initial conditions
T_num[:, N-1] = 100  # Right wall fixed at 100°C
T_num[0, :] = 0      # Initial temperature of the rod is 0°C

for k in range(total_time_steps - 1):
    # Update inner points
    for i in range(1, N-1):
        T_num[k+1, i] = T_num[k, i] * (1 - 2*r) + r*(T_num[k, i+1] + T_num[k, i-1])
    # Update left wall (adiabatic boundary: dT/dx=0 → T0 = T1)
    T_num[k+1, 0] = T_num[k, 1]
    
    # Calculate analytical solution for current time (actual time = k*dt)
    current_t = k * dt
    T_analytical[k] = 100  # Base term T1
    # Sum over odd n only (n=1,3,5...)
    for n in range(1, 200, 2):
        exponent = -a * current_t * ((n * np.pi) / (2*L))**2
        sin_term = np.sin((n * np.pi * x) / (2*L))
        T_analytical[k] += (0 - 100) * (4 / np.pi) * (1/n) * np.exp(exponent) * sin_term

# Handle the last time step for analytical solution
current_t = (total_time_steps-1)*dt
T_analytical[-1] = 100
for n in range(1, 200, 2):
    exponent = -a * current_t * ((n * np.pi) / (2*L))**2
    sin_term = np.sin((n * np.pi * x) / (2*L))
    T_analytical[-1] += (0 - 100) * (4 / np.pi) * (1/n) * np.exp(exponent) * sin_term

# Create animation
fig, ax = plt.subplots()
ax.set_xlim(0, L)
ax.set_ylim(0, 110)  # Slightly higher to avoid clipping
ax.set_xlabel('Position (m)')
ax.set_ylabel('Temperature (°C)')
line_num, = ax.plot([], [], lw=2, color='red', label='Numerical')
line_anal, = ax.plot([], [], lw=2, color='blue', label='Analytical')
ax.legend()

def update(frame):
    line_num.set_data(x, T_num[frame])
    line_anal.set_data(x, T_analytical[frame])
    ax.set_title(f'Time: {frame*dt} s')
    return line_num, line_anal

ani = animation.FuncAnimation(fig, update, frames=total_time_steps, interval=100, blit=True)
ani.save('heat_diffusion_comparison.gif', writer='pillow', fps=10)
plt.show()

额外说明 #

• 你的参数计算出的收敛因子r=0.4，满足显式有限差分法r<0.5的稳定条件，无需调整步长。
• 解析解的级数求和取到n=199（奇数）已足够收敛，无需更大的n值。
• 动画y轴设为0-110，避免温度接近100°C时被边界截断。

内容的提问来源于stack exchange，提问作者Felipe da Costa Kraus