MATLAB求解一二维热传导方程：显式、隐式与Crank-Nicolson实战

搞热传导仿真这些年，我越来越觉得热传导方程这东西就像是老朋友聚会——不管你做机械结构的热应力、电池包的温升管理，还是芯片散热设计，大家总得找个方式坐下来，认认真真聊聊温度到底怎么传。今天这篇就用MATLAB把一维和二维的导热算例从头到尾折腾一遍，重点聊聊不同计算格式的“脾气”：显式格式快但娇气，隐式格式稳但麻烦，Crank-Nicolson则是工程里的常客。我会给出可以直接跑的代码、参数选择思路，以及我踩过的一些坑，适合正在学数值传热、或者刚接触计算传热学并想脱离仿真软件自己验证离散格式的同学。

顺便说一句，这类仿真逻辑搞明白之后，后面想接光学工具箱、Simulink做联合分析，或者读grib数据算环境温度边界，都有个扎实的底子在。我不爱写那种“理论推导完就结束”的文章，咱们边写代码边聊，聊到哪儿算哪儿，但保证都能落地。

1. 热传导方程数值求解的底层逻辑

1.1 方程本身在表达什么：从傅里叶定律到能量守恒

热传导方程的起点是傅里叶定律：单位时间通过单位面积的热量，和温度梯度成正比，方向朝温度降低的方向。

[
q = -k \nabla T
]

这个式子很简单，但它背后站着一条热力学第一定律——能量守恒。把傅里叶定律代入微元体的能量平衡方程，整理出来就是热传导方程的标准形式：

[
\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q
]

如果材料物性不随温度变化，也就是常物性假设，那方程可以进一步化简成：

[
\frac{\partial T}{\partial t} = \alpha \nabla^2 T
]

其中 (\alpha = k/(\rho c_p))，叫热扩散率，单位是 (\text{m}^2/\text{s})。这个参数特别关键，它描述的是“温度变化在材料里传播的速度”。打个比方，金属的 (\alpha) 很大，你拿打火机烧铁棒一端，另一端很快就烫手；而木头的 (\alpha) 很小，同样烧一端，另一端要很久才有反应。

我见过不少刚接触仿真的人，一上来就急着敲代码，方程写出来了却不理解每一项的物理含义。这样一旦结果不对劲，连从哪里排查都不知道。所以我的习惯是先花十分钟把问题的物理场景想明白：热源在哪、初始温度是多少、边界是定温还是绝热、材料是什么。这些信息到最后都会变成代码里的参数和边界条件，一步都省不了。

1.2 离散化的核心思想：网格是把连续世界切成棋盘

解析解只存在于极少数的理想几何和边界条件下，比如无限大平板、无限长圆筒之类。工程上碰到的是一堆不规则形状、多层复合材料、非线性边界条件，只能靠数值方法。有限差分是其中最容易上手的一种思路。

有限差分的核心思想，是用离散点上的温度值来近似连续的温度场。在 (x) 方向把求解区域分成 (N_x) 个点，相邻两个点之间的距离是 (\Delta x)，时间方向上每步前进 (\Delta t)。这样原来那个偏微分方程，就变成了一个代数方程组，可以用计算机一步步推进求解。

具体来说，空间二阶导数的中心差分格式是：

[
\frac{\partial^2 T}{\partial x^2} \approx \frac{T_{i-1} - 2T_i + T_{i+1}}{\Delta x^2}
]

这个格式的截断误差是 (O(\Delta x^2))，精度不错，而且形式简单。时间方向上，根据用前向差分还是后向差分，会衍生出完全不同的计算格式，这就涉及后面要聊的“脾气”问题了。

一个建议是：不要在网格设计上拍脑袋。先把你要分析的最小温度梯度尺度想清楚。比如芯片散热，发热区域可能只有几百微米，网格明显要比发热区域小一个量级，才能捕捉到真实的温度分布。然后做网格无关性验证——把网格加密一倍，看关键位置的温度变化大不大，如果变化很小，就说明你的网格已经够用了。

1.3 数值格式的“脾气”：为什么会有显式、隐式之分

用前向差分近似时间导数，也就是用当前时刻的温度直接算出下一时刻的温度，这叫显式格式。它的好处是简单，一步一个变量，直接代入即可，不需要解方程。

但显式格式有个著名的稳定性限制——傅里叶数必须满足：

[
Fo = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5
]

一旦超出这个门槛，计算就会发散，温度出现锯齿状振荡，最后直接变成NaN。这个限制在生活中可以这么理解：显式格式只允许信息在一个时间步内传播一个网格。如果时间步长太大，温度变化的“信号”在物理上已经跑出网格范围，数值上却来不及传递，就会乱套。

隐式格式则不同。它用下一时刻的温度来计算空间导数，每个节点的方程都牵扯到邻近节点的未知量，需要联立求解一个方程组。看起来麻烦，但换来的是无条件稳定——无论 (\Delta t) 取多大，计算都不会发散。

中间还有Crank-Nicolson格式，它把显式和隐式各取一半，时间方向上用中心差分逼近，空间上取新旧两个时刻的平均。它也是无条件稳定的，而且精度比全隐式高一阶，时间截断误差是 (O(\Delta t^2))。

我把三种格式放在一张表里，方便对比：

看到这张表，你应该就明白为什么说格式有“脾气”了。显式格式像急性子，跑得快但一受刺激就崩；隐式格式像慢性子，每步都要计算矩阵，但稳得一批；Crank-Nicolson则是那种做事稳妥但依然灵活的中间派。工程上三种都有人用，关键看你需要什么。

2. 一维导热算例：三种格式的正面交锋

2.1 算例设定与核心参数计算

用一个经典的算例把三种格式跑一遍。一根长度为 (L = 1,\text{m}) 的钢杆，初始温度全场为 (0^\circ\text{C})，左端温度突然升到 (100^\circ\text{C}) 并保持不变，右端绝热，求温度在杆内的传播过程。

钢的材料参数取典型值：导热系数 (k = 50,\text{W/(m·K)})，密度 (\rho = 7800,\text{kg/m}^3)，比热容 (c_p = 500,\text{J/(kg·K)})。于是热扩散率：

[
\alpha = \frac{k}{\rho c_p} = \frac{50}{7800 \times 500} \approx 1.28 \times 10^{-5} ,\text{m}^2/\text{s}
]

网格取 (N_x = 51)，则 (\Delta x = 0.02,\text{m})。如果显式格式取 (Fo = 0.4)，那么时间步长由稳定条件反过来算：

[
\Delta t = \frac{Fo \cdot \Delta x^2}{\alpha} = \frac{0.4 \times 0.0004}{1.28 \times 10^{-5}} \approx 12.5,\text{s}
]

我故意用这个参数，因为一维杆的导热过程在几十秒到几百秒的时间尺度上才明显。如果拿显式格式跑到两千步，也就六七个小时的实际物理时间，对应温度已经基本稳定了。

2.2 显式格式：简单直观但有个硬门槛

显式格式的代码很简洁。核心就是那个递推公式：

[
T_i^{n+1} = Fo \cdot T_{i-1}^n + (1 - 2Fo) \cdot T_i^n + Fo \cdot T_{i+1}^n
]

写成MATLAB矢量化的形式，一个循环就够了：

matlab复制% 一维显式格式求解热传导方程
L = 1.0;                  % 杆长 m
Nx = 51;                  % 网格数
dx = L/(Nx-1);            % 空间步长 m
alpha = 1.28e-5;          % 热扩散率 m^2/s
Fo = 0.4;                 % Fourier数（必须<=0.5）
dt = Fo * dx^2 / alpha;   % 时间步长 s
Nt = 2000;                % 时间步数

T = zeros(1, Nx);         % 初始温度 0°C
T(1) = 100;               % 左端定温 100°C

for n = 1:Nt
    T(2:end-1) = Fo * (T(1:end-2) + T(3:end)) + (1 - 2*Fo) * T(2:end-1);
    T(1) = 100;           % 左端边界重置
    T(end) = T(end-1);    % 右端绝热：一阶近似
end

figure;
plot(0:dx:L, T, 'LineWidth', 1.5);
xlabel('x (m)');
ylabel('Temperature (°C)');
title('一维杆稳态温度分布（显式格式）');
grid on;

注意边界条件的写法。左端是定温边界，所以每一时间步都要把第一个节点的温度强行赋为100。右端绝热意味着温度梯度为零，最直接的做法是让最后一个节点的温度和倒数第二个节点一样，这就是 (T(\text{end}) = T(\text{end}-1))。这是一阶近似，误差 (O(\Delta x))，在网格足够细时影响不大。但如果你追求更高精度，可以用虚拟节点法，让 (T(\text{end}+1) = T(\text{end}-1))，再代入二阶导数为零的条件，这样精度能提到二阶。

显式格式跑起来确实快，但它的硬门槛在于稳定性。假设你把网格加密一倍，(\Delta x) 变成0.01，那在保持 (Fo) 不变的情况下 (\Delta t) 会缩小到原来的四分之一。也就是说，网格翻倍后你要跑四倍的时间步数，每个时间步的网格点数量也翻倍，总计算量变成八倍。这种线性增长在一维还能接受，到了二维、三维就非常痛苦。这也是为什么工程上很少用显式格式做大规模导热分析。

2.3 隐式格式与追赶法：稳定性的代价是矩阵求解

隐式格式的时间导数是后向差分：

[
\frac{T_i^{n+1} - T_i^n}{\Delta t} = \alpha \frac{T_{i-1}^{n+1} - 2T_i^{n+1} + T_{i+1}^{n+1}}{\Delta x^2}
]

整理一下，得到每个内部节点的方程：

[
-Fo \cdot T_{i-1}^{n+1} + (1 + 2Fo) \cdot T_i^{n+1} - Fo \cdot T_{i+1}^{n+1} = T_i^n
]

所有节点放在一起，就是一个三对角方程组。大名鼎鼎的Thomas追赶法就是专门解这类方程组的，计算量只有 (O(N))，效率极高。MATLAB里不需要自己写托马斯算法，直接用 \ 运算符或者 gallery('tridiag') 生成矩阵即可：

matlab复制% 一维隐式格式求解热传导方程
L = 1.0;
Nx = 51;
dx = L/(Nx-1);
alpha = 1.28e-5;
Fo = 0.8;                 % 隐式可以取更大的Fo
dt = Fo * dx^2 / alpha;
Nt = 1000;

T = zeros(1, Nx);
T(1) = 100;

% 内部节点对应的三对角矩阵（Nx-2个内部节点）
A = gallery('tridiag', Nx-2, -Fo, 1+2*Fo, -Fo);

for n = 1:Nt
    rhs = T(2:end-1)';
    % 边界修正：把已知边界值移到右侧
    rhs(1) = rhs(1) + Fo * T(1);
    rhs(end) = rhs(end) + Fo * T(end);
    T(2:end-1) = A \ rhs;
    % 边界重置
    T(1) = 100;
    T(end) = T(end-1);
end

figure;
plot(0:dx:L, T, 'LineWidth', 1.5);
xlabel('x (m)');
ylabel('Temperature (°C)');
title('一维杆稳态温度分布（隐式格式）');
grid on;

这里我把 (Fo) 取到了0.8，是显式格式承受不了的值，但隐式格式跑得毫无压力。工程上用它有个额外的好处：时间步长可以选得更大，能快速跳过前面那段温度变化剧烈的阶段，直接看长时间行为的趋势。

不过隐式格式也有自己的坑。它的时间精度只有一阶，如果时间步长取得太大，瞬态过程会被严重抹平。比如我想看温度波前从左端传到右端的细节，用大 (\Delta t) 的隐式格式跑出来，可能波前已经糊成一团了。所以隐式格式适合看稳态结果、或者做长时间尺度分析，不适合盯瞬态细节。

2.4 Crank-Nicolson：工程实用平衡点

Crank-Nicolson格式的聪明之处在于，它对时间导数用中心差分，相当于在 (n+1/2) 时刻离散方程。空间二阶导数取新旧时刻的平均值，这样整体在时间方向上的精度是二阶。

离散方程可以写成隐式的形式：

[
-\frac{Fo}{2} \cdot T_{i-1}^{n+1} + (1 + Fo) \cdot T_i^{n+1} - \frac{Fo}{2} \cdot T_{i+1}^{n+1} = \frac{Fo}{2} \cdot T_{i-1}^n + (1 - Fo) \cdot T_i^n + \frac{Fo}{2} \cdot T_{i+1}^n
]

等式两边都有系数矩阵，所以代码里需要两个三对角矩阵：

matlab复制% 一维Crank-Nicolson格式求解热传导方程
L = 1.0;
Nx = 51;
dx = L/(Nx-1);
alpha = 1.28e-5;
Fo = 0.4;
dt = Fo * dx^2 / alpha;
Nt = 1000;

T = zeros(1, Nx);
T(1) = 100;

half = Fo/2;
A = gallery('tridiag', Nx-2, -half, 1+Fo, -half);
B = gallery('tridiag', Nx-2, half, 1-Fo, half);

for n = 1:Nt
    rhs = B * T(2:end-1)';
    % 边界修正
    rhs(1) = rhs(1) + half * T(1);
    rhs(end) = rhs(end) + half * T(end);
    T(2:end-1) = A \ rhs;
    
    T(1) = 100;
    T(end) = T(end-1);
end

figure;
plot(0:dx:L, T, 'LineWidth', 1.5);
xlabel('x (m)');
ylabel('Temperature (°C)');
title('一维杆稳态温度分布（Crank-Nicolson格式）');
grid on;

这个代码有个容易踩坑的地方：边界修正里用的是 (half)，也就是 (Fo/2)，不是 (Fo)。我见过不少人在从全隐式改到CN格式的时候，忘了把边界项从 (Fo) 换成 (Fo/2)，导致前几个时间步看着还行，后面温度曲线就开始扭曲。这种隐蔽错误最磨人，排查了半天才发现是边界修正写错了。

CN格式在工程上的定位很明确：当你在乎瞬态过程的精度，又不希望被显式格式的稳定性条件卡死，就用它。虽然每步要多做一次矩阵乘法，但换来二阶精度，还是值得的。

三种格式跑出来的稳态结果其实一模一样，都是从那根杆的左端到右端，温度从100度逐渐降低到某个值。真正的差别在瞬态阶段——显式格式在满足稳定条件时能给出清晰的波前传播过程；隐式格式在这个时间尺度上会显得比较钝；CN格式则兼顾了清晰和稳定。

3. 二维导热算例：从一维走到二维的复杂度跃升

3.1 二维方程的离散与直接求解的困境

二维热传导方程比一维复杂不止一个量级：

[
\frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right)
]

离散化之后，每个内部节点的方程会牵扯到周围四个邻居节点。如果用全隐式格式，形成的矩阵不再是三对角，而是一个五对角矩阵。网格规模是 (N_x \times N_y)，矩阵的规模就是 ((N_x N_y)^2)。假设网格是 (100 \times 100)，矩阵就有1亿个元素，虽然大部分是零，但直接用稠密矩阵求解，内存直接爆掉。

换成显式格式呢？二维显式的稳定性条件更苛刻：

[
Fo_x + Fo_y \le 0.5
]

也就是说两个方向的傅里叶数加起来要小于0.5。一维显式要求单个 (Fo \le 0.5)，二维直接把限制减半，(\Delta t) 的大小进一步缩水。所以做二维瞬态导热，全显式和全隐式都不是好选择——显式太脆，全隐式太重。

这个困境其实在工程里很常见。等你把网格加到 (200 \times 200) 甚至更高，全隐式每步要解一个巨大的稀疏矩阵系统，虽然MATLAB的稀疏矩阵求解效率不错，但架不住成千上万的迭代步，算一次要等到天荒地老。

3.2 ADI格式：把二维问题拆成两个一维问题

工程上破解二维困境的经典解法是ADI，也就是交替方向隐式格式。核心思想很巧妙：把一个大矩阵求解问题，拆成两个方向上一维的三对角求解问题。

具体做法是每个时间步分成两个半步：

第一个半步，在 (x) 方向用隐式格式，在 (y) 方向用显式格式。这时你需要解的是 (N_y-2) 个三对角方程组，每个方程组的规模是 (N_x-2)。

第二个半步，改在 (y) 方向隐式，(x) 方向显式。同样解 (N_x-2) 个三对角方程组，每个规模是 (N_y-2)。

两个半步合起来，时间精度可以达到二阶，而稳定性也是无条件稳定。这种“分而治之”的思路，让二维问题在计算量和稳定性之间找到了极好的平衡点。

核心代码框架大致是这样：

matlab复制% 二维ADI（P-R格式）核心循环
% 假设T是(Nx, Ny)矩阵，边界已按问题设定初始化
Fox = alpha * dt / (2 * dx^2);
Foy = alpha * dt / (2 * dy^2);

Ax = gallery('tridiag', Nx-2, -Fox, 1+2*Fox, -Fox);
Ay = gallery('tridiag', Ny-2, -Foy, 1+2*Foy, -Foy);

for n = 1:Nt
    % 第一个半步：x方向隐式，y方向显式
    for j = 2:Ny-1
        rhs = T(2:Nx-1, j) + Foy * (T(2:Nx-1, j-1) - 2*T(2:Nx-1, j) + T(2:Nx-1, j+1));
        % 加入x方向的边界修正，例如定温边界的贡献
        rhs(1) = rhs(1) + Fox * T(1, j);
        rhs(end) = rhs(end) + Fox * T(Nx, j);
        T(2:Nx-1, j) = Ax \ rhs;
    end
    
    % 第二个半步：y方向隐式，x方向显式
    for i = 2:Nx-1
        rhs = T(i, 2:Ny-1) + Fox * (T(i-1, 2:Ny-1) - 2*T(i, 2:Ny-1) + T(i+1, 2:Ny-1));
        % 加入y方向的边界修正
        rhs(1) = rhs(1) + Foy * T(i, 1);
        rhs(end) = rhs(end) + Foy * T(i, Ny);
        T(i, 2:Ny-1) = Ay \ rhs;
    end
end

这个代码框架里有一个很容易出错的地方：两个半步中的显式项系数用的是 (Fox) 和 (Foy)，不是 (2Fox)。因为每个半步的时间长度是 (\Delta t / 2)，代入定义后，显式项和隐式项里的系数正好是 (\alpha \Delta t / (2 \Delta x^2))。如果搞混了，结果会差上一倍。我最初写ADI的时候，在这里卡了差不多两天，最后是拿一个已知解析解的例子对比才揪出来。

关于边界条件的处理，ADI有个微妙的地方：每个半步结束时，边界上的温度怎么取？常见做法是边界值保持不变，或者根据物理条件显式更新。但如果边界本身也在变化，两个半步之间需要小心处理，否则会引入额外误差。工程上我一般先跑一个简单的等温边界算例，验证码跑出来的结果和解析解一致，再做复杂边界。

3.3 稳态与瞬态的边界处理技巧

如果你只需要二维稳态温度场，其实用不着ADI这种瞬态推进方法，直接解稳态方程更省事。稳态方程是拉普拉斯方程：

[
\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = 0
]

常用迭代法求解：Jacobi、Gauss-Seidel、SOR。这三种方法的收敛速度差别很大。Jacobi迭代收敛最慢，Gauss-Seidel会快一些，SOR加一个松弛因子 (\omega)，通常取1.2到1.8，能显著加速收敛。但注意，最优松弛因子依赖网格和几何形状，没有通用值，需要试算。

我有个习惯：先跑几百步SOR迭代看残差下降速度，如果残差曲线在振荡，说明 (\omega) 取大了；如果下降太慢，说明 (\omega) 偏小。这样调两三次就能找到合适的值。判断收敛的准则一般是检查相邻两次迭代的温度差最大值，小于某个阈值比如 (10^{-6}) 就算收敛。

瞬态计算里的边界处理则要区分清楚物理类型。第一类边界也叫Dirichlet边界，给定壁面温度值；第二类边界是给定热流密度，绝热边界就是热流为零的特殊情况；第三类边界是流体与壁面之间的对流换热，常见于散热器设计，边界方程里会多出一个对流换热系数 (h) 和流体温度 (T_\infty)。这些在代码里对应不同的赋值方式，稍不注意就会出错。

4. 工程仿真中的参数处理与边界条件选择

4.1 热扩散率、网格与时间步长的配合

工程仿真里，参数之间的配合比单点参数重要得多。热扩散率、网格尺寸和时间步长这三者，被傅里叶数紧紧绑在一起。

[
Fo = \frac{\alpha \Delta t}{\Delta x^2}
]

这张关系表经常被我贴在工位上：

从这个表可以看得很明白：为什么金属导热用显式格式这么痛苦。铝的 (\alpha) 那么大，(\Delta x) 只要取到1毫米，显式格式允许的最大时间步长就只有几毫秒。如果你要模拟一个几秒钟的过程，几千步只是起步。这就是为什么金属热分析的瞬态过程，工程上普遍用隐式或者CN格式。

做具体算例时还有个经验：别一开始就用细网格。先用粗网格把整体趋势跑通，验证代码逻辑和边界条件没有错误，然后逐步加密网格做网格无关性验证。这样能省下大量调试时间。

4.2 三类边界条件的物理含义与代码实现

第一类边界条件，定温边界，代码里最简单，就是每步赋值：

matlab复制T(1,:) = 100;    % 左边界保持100°C
T(end,:) = 50;   % 右边界保持50°C

第二类边界条件，定热流边界。在左边界有一恒定热流 (q) 进入时，离散形式是：

[
k \frac{T_2 - T_1}{\Delta x} = q
]

反解出 (T_1 = T_2 - q \Delta x / k)。绝热边界就是 (q = 0) 的特例，直接让 (T_1 = T_2)。

第三类边界条件，对流边界。壁面与流体之间的换热满足：

[
-k \frac{\partial T}{\partial x}\bigg|{\text{wall}} = h (T\text{wall} - T_\infty)
]

离散化后，壁面温度被流体温度和导热、对流两者之间的相对强弱决定。这个公式在散热分析里特别常见，比如电子设备外壳自然冷却、汽车散热器表面换热，都是第三类边界。

第三类边界的一个常见陷阱是 (h) 的量纲和对流系数的取值范围。自然对流 (h) 大约在 (5 \sim 25 , \text{W/(m}^2\cdot\text{K)})，强制对流可以到几十甚至上百。如果数值取得不靠谱，计算结果会完全偏离现实。我的建议是：宁可先查文献、查工程手册里的经验值，也别靠感觉填数。

4.3 后处理与结果检查：怎么判断算对了

很多人跑完代码看一个最大温度就完事，这是远远不够的。后处理不只是画图，更是你验证计算结果合理性的关键手段。

第一件事是看温度场的形状。一维算例，画 (T-x) 曲线，观察有没有局部尖峰或锯齿。二维算例，用 surf 或者 pcolor 画热力图，配合 contourf 画等温线，能直观地捕捉到温度分布是否平滑、热点位置是否合理。

matlab复制% 二维温度场可视化示例
figure;
surf(X, Y, T');
xlabel('x (m)');
ylabel('y (m)');
zlabel('Temperature (°C)');
shading interp;
colorbar;
title('二维平板温度分布');

第二件事是检查能量守恒。对于一个内部无热源、边界绝热的问题，初始时刻的总热能应该等于任意时刻的总热能。离散后可以做一个全局的能量积分，算完所有时间步后对比始末总能量，偏差应该在1%以内。如果偏差很大，多半是边界处理或者矩阵求解出了问题。

第三件事是拿解析解验证。一维半无限大物体在表面温度阶跃变化时，温度分布有著名的误差函数解；二维矩形平板在特定边界条件下也有傅里叶级数解。你的数值解可以和这些解对比，误差在几个百分点以内，说明代码是可靠的。这个验证步骤，我认为是数值仿真中最重要的一步，没有之一。做过这个验证，后面换个复杂几何心里才有底。

第四件事是动态观察。用MATLAB的 animatedline 或者逐帧更新图像，把瞬态过程播放出来。你会发现温度波前怎么传播、哪些区域先热起来、哪些区域始终冷却，这些信息对工程优化非常有价值。唯一要注意的是，MATLAB图形更新比较费时，大网格时可以先算完存数据，再统一后处理。

5. 常见问题与调试实录

5.1 发散、振荡、守恒性问题的排查

每个跑过数值仿真的人，都经历过对着屏幕上的NaN发呆的时刻。我做热传导仿真这些年，踩过的坑集中在这几类，整理成一个速查表：

关于矩阵奇异的问题多说两句。隐式格式的三对角矩阵通常是对角占优的，但你如果边界条件处理不当，比如在内部节点里错误地包含了边界未知量，矩阵结构就被破坏了。MATLAB在求解奇异矩阵时会给出警告，很多人忽略了这个警告，直接拿到一堆NaN。我现在的习惯是：解矩阵之前先用 condest(A) 或者 rcond(A) 检查一下矩阵的条件数，心里有底之后再做下一步。

关于单位不一致，我在刚工作那会儿犯过低级错误：把热扩散率写成 (1.28 \times 10^{-5} , \text{m}^2/\text{s})，结果时间步长取成了秒，但网格尺寸用毫米。一换算，结果差了六个量级，温度场莫名其妙。后来我给自己定了个规矩：所有代码统一用SI单位制，材料参数用一个 params.m 文件集中管理，写清每个变量的单位，跑完计算再根据需要转换成工程单位。这样做以后，这种低级错误基本绝迹了。

5.2 从算例到工程应用的扩展建议

这套一维二维基板热传导的代码思路，往工程方向扩展的空间其实很大。

首先是材料物性随温度变化的情况。高温环境的导热问题，比如热成型、焊接热循环，材料的导热系数、比热容都不能再当作常数。这时候在每个时间步根据当前温度更新物性参数，方程变成非线性，求解时需要迭代。隐式格式和CN格式每步要多几轮迭代，计算量增加，但思路没变。

其次是带内热源的场景。电池放电产热、芯片焦耳热、相变材料吸热等，都可以在离散方程里加一个源项 (Q)。源项如果本身依赖温度，比如辐射换热，那就把它线性化处理，放一半在已知项、一半在未知项，格式仍然稳定。这个技巧在工程热仿真里特别常用。

第三是几何复杂度的提升。坐标轴对齐的矩形网格只能处理简单几何，碰到带圆角、斜面的结构，就得用非均匀网格或者有限元、有限体积方法。不过，有限差分算例的价值在于帮你建立数值离散的感觉——什么是网格无关性、什么是稳定性条件、什么是边界条件离散。这些底层素养，放到任何数值工具里都管用。

另外提一句MATLAB生态。这套热传导求解逻辑，之后可以和MATLAB里的光学工具箱配合，做光热耦合分析；也可以跟Simulink搭个联合模型，把温度场的影响传给控制系统的反馈。把基础算例跑明白，就好比给自己的工具箱配齐了扳手和螺丝刀，后面想组装什么复杂结构都顺手。

在实际操作中，我一般会把这套格式封装成几个函数——solver_1D_explicit.m、solver_1D_implicit.m、solver_1D_cn.m、solver_2D_adi.m，每个函数接收物性参数、网格参数、边界条件，返回温度场的时间序列。这样以后换一种材料、换一组边界，只需要改参数传入，不折腾核心代码。文件命名规则也尽量统一，省得过了三个月回头看，自己都分不清哪个函数是哪一版。

初学者最容易犯的错误，是直接在一个脚本里从参数定义写到后处理，几十行代码揉成一团。出问题时根本分不清是物理设定错了还是代码逻辑错了。我比较推荐的做法是把参数区、求解区、后处理区用注释分成清晰区块，每个区块单独检查。参数区打印一遍，确认数值合理；求解区先跑一个已知解析解的算例验证；后处理区再画图。这套工作流看起来麻烦，但在调试时省下来的时间，远远超过多写几行注释的代价。

最后再分享一个我个人的小习惯：做仿真之前先把物理问题用文字写下来，包括初始条件、边界条件、材料参数、期望观测的时间尺度、想得到的关键输出指标。写到一半如果发现某些条件自己定不出来，说明你对问题的理解还有漏洞，这时候回头补物理远比你跑出一个似是而非的结果要高效。

这些算例本身不难，但把格式的脾气摸清楚、把参数的坑绕过去，需要一点一点积累。希望这篇折腾记录能让你少走点弯路。