企业官网建设流程全解析

简介：这是一份基于 MATLAB 的一维热传导偏微分方程求解源码，围绕傅里叶热传导方程展开，适合正在学习数学物理方程、数值计算或传热学基础的高校学生与科研人员使用。压缩包内共 4 个 M 文件，分别对应主程序、方程定义、初始条件与边界条件设置，整体仅 1KB，代码精简、逻辑清晰，稍加修改即可适配其他热传导初边值问题。已有 249 人学习下载。读者通过研读代码，可以掌握有限差分法对空间二阶导数和时间导数的离散化思路，理解不同边界条件与初始条件在 PDE 数值求解中的具体写法，并建立一维热传导问题的完整求解框架；对希望快速入门 MATLAB 求解偏微分方程的初学者而言，这份小巧的源码也是一份可直接运行和二次开发的学习样例。

1. 时间步一放大就发散，问题不出在代码上

一维热传导方程是所有偏微分方程（PDE）数值解里最像“Hello World”的问题，但它一点不比流体力学里的方程好伺候。只改一个时间步长 dt，显式格式的解就能从平滑曲线变成锯齿噪声，几分钟内温度冲到 1e30，这种崩溃不是浮点误差，也不该靠减小 dt 去“硬解”。真正的原因是离散格式的稳定性条件被突破了，而大多数现成求解器不会主动告诉你这一点。

本篇文章要做的，就是把“一维热传导 + 偏微分方程数值解”这条线完整打通：从连续方程出发，用有限差分把它离散成可迭代的代数方程，给出显式和隐式（Crank-Nicolson）两套可运行的 Python 实现，然后把边界条件、源项和稳定性判定这些工程问题讲清楚。

这个领域的核心难点不是推导，而是三个问题：如何选离散格式、如何设 dx 和 dt、如何确认结果可信。下面五章会逐个击破。适合刚上手 PDE 数值计算的工程师，也适合那些调用过现成求解器、但遇到发散现象后需要搞清内部机制的人。

2. PDE 离散化：从连续方程到可迭代的差分格式

2.1 一维热传导方程的定解问题与量纲

一维热传导方程的标准形式是：

ut = α · uxx

其中 u(x,t) 表示 t 时刻位置 x 处的温度，α 是热扩散系数，单位是 m²/s。它的物理含义是：某一点温度的变化率，取决于该点温度分布的二阶导数。温度分布是下凹的，温度就上升；是上凸的，温度就下降——这是抛物型 PDE 最直观的特征，也是“扩散”一词的由来。铝的 α 大约在 9.7e-5，水约 1.4e-7，选择不同的 α，时间尺度完全不同，量纲检查是第一道防线：α·dt/dx² 必须是无量纲的。

要完整求解这个 PDE，还需要三类输入。初始条件 u(x,0) 写明初始温度分布；边界条件说明端点处是固定温度（狄利克雷）还是固定热流（诺伊曼）；几何参数确定杆长 L 和网格划分数 nx。三者组合不同，工程上的含义完全不同——绝热边界和恒温边界算出来的是两种物理过程，不能混用。

2.2 显式欧拉迭代：最直觉的离散路径

求解 PDE 数值解的常用做法是有限差分法，用差分商替代偏导数。对时间用向前差分，对空间用二阶中心差分，得到显式格式：

u[i]^(n+1) = u[i]^n + r · (u[i-1]^n − 2u[i]^n + u[i+1]^n)

其中 r = α·dt/dx²。这个式子之所以叫“显式”，是因为新时间层的值只依赖旧时间层，直接赋值就能推进。以下是完整的可运行代码：

import numpy as np def solve_explicit(u0, alpha, dx, dt, nt): r = alpha * dt / dx**2 if r > 0.5: raise ValueError(f"r = {r:.3f} 超出显式格式稳定性上限 0.5") u = u0.copy() nx = len(u) for n in range(nt): u_new = u.copy() u_new[1:-1] = u[1:-1] + r * (u[:-2] - 2*u[1:-1] + u[2:]) # 狄利克雷边界，两端恒温为 0 u_new[0] = 0.0 u_new[-1] = 0.0 u = u_new return u

这段代码的核心是向量化差分计算。u[:-2]、u[1:-1]、u[2:] 构成三组错位的数组切片，分别对应左邻居、自身和右邻居，避免使用 Python for 循环遍历每个网格点。内部点更新完后再强制覆盖边界值，确保边界条件在每个时间步都严格成立。

r 是显式格式唯一的控制参数。代码开头就做了检查，因为 r>0.5 时误差会逐层放大，这不是减小 dt 能解决的，必须换格式。物理直觉是：一个时间步内热量从一个网格传到另一个网格的距离不能超过半个网格间距。

2.3 稳定性红线：为什么 r>0.5 数值解会发散

显式格式的稳定性受 CFL 条件约束，对一维热传导方程具体表现为 r ≤ 0.5。这个界限可以直接从误差传播矩阵的特征值推出来：放大因子 G = 1 − 4r·sin²(k·dx/2)，当 r>0.5 时，高频分量的放大因子绝对值会超过 1，误差指数增长。观察高频分量的时候，数值解出现的就是高频锯齿振荡：

如果发现温度分布出现 ±1 交替的锯齿状，而网格画出来不光滑，大概率是 r 越过了 0.5。可以在这个时间步的循环里加一个检查：

if np.max(np.abs(u)) > 1e8: print(f"第 {n} 步数值解发散，当前最大值 {np.max(np.abs(u)):.3e}") break

显式格式的代价在这里很清晰：想精确一点就把 dx 减半，但 r 不变的约束下 dt 必须缩到原来的 1/4，总计算量变成原来的 8 倍。实际项目里，这就是我一般会转向隐式格式的直接原因——不是为了炫耀数学技巧，而是为了少等那几个小时。

3. Crank-Nicolson 隐式迭代：PDE 解算的核心格式

3.1 隐式格式为什么能摆脱时间步限制

显式格式的稳定性限制来自把空间差分项放在了旧时间层。如果把 uxx 的一部分或者全部放到 n+1 层，就得到了隐式格式。完全隐式格式（全取 n+1 层）的放大因子是 G = 1 / (1 + 4r·sin²(k·dx/2))，对任意 r 都小于 1，这就是无条件稳定的来由。

但完全隐式只有一阶时间精度，收敛慢。工程上更常用的折中方案是 Crank-Nicolson：把空间二阶导数在 n 层和 n+1 层取平均，时间二阶精度，空间二阶精度，同时继承无条件稳定性。离散形式是：

u[i]^(n+1) − (r/2) · (u[i-1]^(n+1) − 2u[i]^(n+1) + u[i+1]^(n+1)) = u[i]^n + (r/2) · (u[i-1]^n − 2u[i]^n + u[i+1]^n)

这个等式左边涉及三个相邻点的未知值，不能逐个点推进，必须联立求解一个三对角线性方程组。

3.2 用稀疏矩阵构造三对角线性系统

把离散方程整理成矩阵形式 Au^(n+1) = b。矩阵 A 的主对角是 1+r，上下次对角是 −r/2。这个结构在问题规模增大时非常重要：如果 nx=10000，A 有 1e8 个元素，但非零元素只有约 3e4 个，用稀疏矩阵存储并求解，内存和时间开销都只有稠密方法的几个百分点。这是大规模 PDE 求解的底线思维：随时用稀疏结构，不要因为规模小就忽视这个习惯。

构造稀疏矩阵的推荐方式是使用 scipy.sparse。先定义三对角结构，再处理边界行。这里用 lil_matrix 修改行，因为它在按行赋值时的开销远低于 csr 或 csc 格式：

from scipy.sparse import diags, lil_matrix def build_cn_matrix(nx, r): A = lil_matrix((nx, nx)) main_diag = 1 + r off_diag = -r / 2 A.setdiag(main_diag) A.setdiag(off_diag, k=-1) A.setdiag(off_diag, k=1) # 恒温边界：两端点的方程变为 u[0]=0, u[-1]=0 A[0, :] = 0 A[0, 0] = 1 A[-1, :] = 0 A[-1, -1] = 1 return A.tocsr()

参数说明：setdiag 的三个参数分别设置主对角线和两条次对角线，k=-1 表示下对角线，k=1 表示上对角线。边界行的处理方式是逻辑核心：把矩阵第一行和最后一行替换为 u[0]=0、u[-1]=0，这等价于在离散方程组中直接嵌入狄利克雷边界条件，比每步迭代后覆盖边界值要优雅得多，尤其是在处理非零边界值时只需要改这一行。

3.3 完整迭代循环与参数表

右侧向量 b 的计算要把旧时间层 u 的扩散项搬过去，同时处理边界值。狄利克雷边界值为 0 时，b[0] 和 b[-1] 直接赋 0。完整求解过程如下：

def solve_cn(u0, alpha, dx, dt, nt): nx = len(u0) r = alpha * dt / dx**2 A = build_cn_matrix(nx, r) u = u0.copy() for n in range(nt): b = np.zeros_like(u) b[1:-1] = (1 - r) * u[1:-1] + (r / 2) * (u[:-2] + u[2:]) b[0] = 0.0 b[-1] = 0.0 u = spsolve(A, b) return u

对 scipy 提供的稀疏直接求解器 spsolve 来说，三对角矩阵的求解复杂度是 O(n)，即使跑一万步、上千个网格点，整个过程也在秒级完成。spsolve 用的是稀疏 LU 分解，如果 A 在迭代过程中保持不变（r 不变时成立），还可以分解一次后反复使用，代码里可以复用分解结果来加速：

from scipy.sparse.linalg import splu lu = splu(A) u = lu.solve(b)

这项优化对长时间模拟尤其有效。下表对比了两种格式与参数的关系，作为选型的快速依据：

Crank-Nicolson 虽然也需要解方程，但换来的时间步自由度足够弥补这个开销。实际模拟中可以用比显式格式大几十倍的 dt，总计算量反而下降一个量级。

4. 边界条件与工程参数：让热传导模型贴近真实工况

4.1 狄利克雷与诺伊曼边界在代码里的差别

实际工程中，恒温边界很少见，更常见的是绝热边界、热流边界和换热边界，这三者都属于诺伊曼类型，约束的是空间导数 ux，而不是温度值本身。

绝热边界最简单的理解是：端点与外界没有热交换，温度梯度为零，即 ux(0,t)=0。离散化后用一阶差分近似，得到 u[0]=u[1]。这个处理在显式格式里只需在每步迭代后加一句 u[0] = u[1]。在 Crank-Nicolson 格式中，需要修改矩阵的第一行：

A[0, 0] = 1 A[0, 1] = -1 # 替代原来的 u[0] = 0 写法 b[0] = 0

这行的意思是 u[0] − u[1] = 0，对应的正是差分近似。表面上看只是一个赋值，但矩阵行的含义变了：狄利克雷边界是硬性固定温度，诺伊曼边界是让方程自己去满足端点的热流约束。CN 里如果还把边界行写成 u[0] = 0，实际模拟的就是端点一直保持 0 度，与绝热假设完全矛盾，而且结果不会立刻报错，只在出图时发现温度分布不对，这种坑最难排查。

4.2 热流边界与源项的叠加方法

有恒定热流 q 注入时，边界条件变为 −k·ux = q，其中 k 是导热系数。离散后得到 u[0] − u[1] = −q·dx/k。在 CN 矩阵里实现为：

A[0, 0] = 1 A[0, 1] = -1 b[0] = q * dx / k

如果体力本身有热源或散热，比如电热丝加热或化学反应放热，还要在方程右边加源项 f(x,t)。在离散方程中，源项要同时影响时间层 n 和 n+1，最简单且精确的做法是取平均值：把 b 的定义式修改为在右端加上 0.5·Δt·f。

4.3 初始条件构造与 dt/dx 的设定依据

初始条件决定了模拟的物理场景。如果是模拟“一根均匀温度的杆，中间某点突然受热”，可以用高斯函数描述这个局部热脉冲：

x = np.linspace(0, 1.0, 201) x0, sigma = 0.5, 0.05 u0 = np.exp(-(x - x0)**2 / (2 * sigma**2))

这个初始分布展示了一步到位的高斯脉冲平滑过程，也是后面第五章验证解析解的标准测试场景。网格设定时先确定 dx，再根据稳定性和精度要求调整 dt。显式格式优先检查 r ≤ 0.5；CN 格式理论上无限制，但 dt 过大时时间离散误差会恶化解的质量，一般建议 r 控制在 5 到 10 以内，经验规则是保证 dt 不超过最小物理特征时间（如脉冲扩散跨越一个网格的时间）的几分之一。跑完一次后，把 dt 减半再跑一次，如果最终分布几乎不变化，就说明当前时间步已经合适，这是标准的网格无关性检查。

5. 验证数值解的正确性：解析解、能量守恒与收敛阶

比起解方程，更难的是确认方程解得对。解析解是最好的基准工具。一维无限域上高斯脉冲的解析解是：

u(x,t) = σ/√(σ² + 2αt) · exp(−(x−x0)²/(2(σ² + 2αt)))

这个解给出任意时刻的温度分布，并且刻画了脉冲“变矮变宽”的过程。数值解做对比检验时，可以设置一个足够大的杆长，让边界效应在模拟时间内不波及计算域，把数值解与解析解放到同一图上对比，并计算两者之间的 L2 误差：

u_analytic = sigma / np.sqrt(sigma**2 + 2*alpha*T) * \ np.exp(-(x - x0)**2 / (2 * (sigma**2 + 2*alpha*T))) err = np.sqrt(np.mean((u_cn - u_analytic)**2))

除了对标解析解，还有更通用的诊断方法：能量守恒。在一维热传导方程中，如果边界是绝热的，总热量 Σu·dx 应该保持不变，否则就是边界条件出错。这个检查只需要在迭代循环外对最终结果做一次 sum，是排查边界处理的利器。

验证收敛阶是进阶做法。把网格数从 41 依次翻倍到 81、161、321，记录误差值，如果每次网格加密一半时误差约缩小到 1/4，就说明空间离散达到了二阶精度，这就是 Crank-Nicolson 应有的表现：

for nx in [41, 81, 161, 321]: dx = L / (nx - 1) dt = 0.5 * dx**2 / alpha # 保持 r 恒定 # 运行求解器，计算误差

注意保持 r 恒定这个细节：如果同时改变 dx 和 dt，收敛阶会被时间误差污染。网格加细时 dt 也应按比例缩小，才能单独验证空间方向的收敛性。对显式格式做同样的测试时，误差会约缩小到 1/2，这是它时间方向只有一阶精度的信号。

当你面对一个新来的 PDE 求解代码，最有效的验收流程就是这三连：先看能量是否守恒，再对标一个可解析解的特殊工况，最后用网格减半法确认收敛阶。这样就不会把一个写错的格式和写错的参数留着当隐患了。

本文还有配套的精品资源，点击获取