一维热传导PDE求解：显式欧拉法正常、隐式欧拉法初边值设置异常咨询

当前实现的核心异常根源 #

• 离散推导完全错位：从控制方程移项组装矩阵、右端项的时候，把隐式格式的系数关系搞反了，导致矩阵本身描述的方程和原偏微分方程不匹配，求解出来的结果自然没有物理意义。
• 边界条件未真正嵌入求解系统：写矩阵的时候，x=0的定温节点、域末端的零梯度节点都没有替换对应的方程，只是把内部节点的系数抄给了边界行，导致边界约束根本没生效。
• 右端项组装逻辑混乱：把已经确定的边界节点值也当成未知量纳入了pS的计算，没有将边界已知值直接赋值给方程右端，进一步放大了求解误差。

隐式欧拉法的正确实现逻辑 #

首先明确标准隐式离散形式：针对控制方程∂T/∂t = alpha * (∂²T/∂x²)，隐式欧拉格式将时间导数取n时刻差值、空间二阶导取n+1时刻的网格值，得到：

(T^{n+1}_i - T^n_i)/dt = alpha * (T^{n+1}_{i+1} - 2T^{n+1}_i + T^{n+1}_{i-1})/dx²
定义网格傅里叶数r = alpha*dt/dx²，将所有n+1时刻的项移到方程左端，n时刻项移到右端，整理得到内部节点的通用方程：
r*T^{n+1}_{i-1} - (1 + 2r)*T^{n+1}_i + r*T^{n+1}_{i+1} = -T^n_i

初边值的正确嵌入方式 #

边界条件必须直接写入三对角矩阵的对应行，不要单独留空再移到右端，这是抛物型方程数值求解最稳妥的方案：

1. x=0位置的定温边界（T=100℃）
   • 对应矩阵第1行的方程替换为1*T^{n+1}_1 = 100，即矩阵元素MAT(1,1) = 1、MAT(1,2) = 0，对应右端项pS(1) = 100。每一个时间步求解时，这一行会强制第一个节点温度固定为100，满足定温边界要求。
2. 域末端的零梯度边界（T(x=L) = T(x=L-1)）
   • 对应矩阵最后一行的方程替换为-1*T^{n+1}_{x_count-1} + 1*T^{n+1}_{x_count} = 0，即矩阵元素MAT(end, end-1) = -1、MAT(end, end) = 1，对应右端项pS(end) = 0。这一行直接强制最后一个节点和相邻节点温度相等，满足零梯度边界的物理要求。

核心疑问解答 #

1. 如何实现时间充分长后全域温度收敛到100℃的预期效果？
   • 只要将x=0的定温边界正确嵌入矩阵求解系统，末端采用零梯度（绝热）边界，系统的稳态解就必然收敛到全域100℃。因为域末端没有热损失通路，定温边界持续向域内传热，最终整个计算域的温度会和边界温度持平，该收敛特性和dx、dt的取值无关——隐式欧拉格式本身是无条件稳定的，不会出现结果对步长敏感的随机波动。
2. 是否应当将I.C.与B.C.直接写入三对角矩阵的对应位置？
   • 应当这么做。之前尝试该方案失败的核心原因，不是方法本身有问题，而是写矩阵时没有把内部节点的热传导方程替换成边界条件，只是在原有内部节点的系数上叠加了边界赋值，导致方程本身仍然是错误的。直接把边界方程写入矩阵行的方案，比“把边界项移到右端、矩阵行留空”的方案更不容易出现符号错误，数值稳定性更强。

代码修正的关键调整点 #

• 三对角矩阵MAT只需要在时间循环外生成一次即可，因为隐式格式的矩阵系数只和dx、dt、热扩散系数有关，不随时间、温度值变化，不需要每个时间步重复赋值。
• 修正矩阵系数的推导错误：原来的主对角线系数B符号和量级完全错位，替换为基于r的正确系数，即内部节点主对角线为-(1 + 2*r)，上下次对角线为r。
• 右端项pS仅对内部节点用n时刻的温度值计算，边界节点的右端项固定为已知的边界值，不需要用n时刻温度计算。
• 去掉循环中重复给T2边界节点提前赋值的逻辑，因为求解线性方程组后，第一节点会自动根据矩阵行的约束得到100℃，最后一个节点会自动根据零梯度约束和相邻节点温度相等，提前赋值反而会引入干扰。

内容的提问来源于stack exchange，提问作者Josh E.