C++有限差分实现热传导与泊松方程：边界、稳定性与收敛阶

简介：偏微分方程数值解差分方法C++程序代码面向高校数值分析、科学计算课程的学习者及自学者，聚焦椭圆型、抛物型和双曲型方程的离散求解，把网格划分、Dirichlet与Neumann边界处理、差分算子构造、隐式迭代和显式推进等环节落实到可编译的代码中。压缩包共6个文件，均为cpp源程序，整体约6KB，体量轻便，适合逐个编译调试、对照课堂推导修改参数，并观察雅可比迭代、高斯-塞德尔迭代以及CFL条件对计算结果的影响。已有2036人学习，说明它在大学课程实验与自学场景中具备一定参考价值。读者可借助这些程序理解截断误差、舍入误差与数值稳定性的实际表现，把格式推导转化为运行结果，并用于课程作业、实验报告或进一步扩展数值算例。

1. 从课程实验的一张差分离散表说起：这套 C++ 代码能跑什么

如果你正在上 c++ 大学课程里的偏微分方程数值解实验，拿到一套差分方法 C++ 程序代码，第一反应通常不是“真香”，而是“这堆数组到底对应哪张差分表”。这套资源解决的正是这个断点：把一维/二维热传导、泊松方程这类经典模型，用有限差分写成可编译、可改参数、可做收敛阶验证的 C++ 工程。它适合课程实验、期末大作业和自学数值分析的人，也适合已经写过矩阵求解、但没把网格、边界、稳定性串起来的开发者。代码不是黑匣子，重点是把差分离散、系数矩阵、时间推进和误差验证四件事接上。

2. 差分格式落成 C++：网格、系数矩阵与时间推进的选型

差分方法的核心不是背公式，而是把连续算子换成网格上的代数关系。选型要围绕三个问题：网格怎么排、未知量怎么存、时间推进用显式还是隐式。课程实验常见的模型是一维热传导和二维泊松方程，前者看时间推进稳定性，后者看线性系统求解。下面按这三个问题拆，先把最小闭环写出来，再谈参数和边界。

2.1 先定网格和离散格式：一维/二维热传导与泊松方程的最小闭环

一维热传导写成 (u_t=\alpha u_{xx})，显式 FTCS 的离散是 (u_i^{n+1}=u_i^n+r(u_{i-1}^n-2u_i^n+u_{i+1}^n))，其中 (r=\alpha\Delta t/\Delta x^2)。二维泊松方程 (-\Delta u=f) 用五点差分，内部点满足 ((4u_{i,j}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1})/h^2=f_{i,j})。这两个格式足够覆盖课程实验里大部分“写差分、解矩阵、看误差”的需求。

网格先约定清楚：一维节点数nx，间距dx=(x1-x0)/(nx-1)；二维用u[j][i]，j是 y 方向行，i是 x 方向列。边界节点保留在数组里，内部未知量是nx-2或(nx-2)*(ny-2)。这个约定一旦定下，后面所有函数、输出和误差计算都按它走，不要一会儿[i][j]一会儿[j][i]。

#include <vector> #include <cmath> #include <stdexcept> struct Grid1D { int nx; // 节点数，不是单元数 double x0, x1, dx; // 区间起点、终点、空间步长 std::vector<double> x; }; Grid1D make_grid(int nx, double x0, double x1) { if (nx < 3) throw std::invalid_argument("nx must be >= 3"); Grid1D g; g.nx = nx; g.x0 = x0; g.x1 = x1; g.dx = (x1 - x0) / static_cast<double>(nx - 1); g.x.resize(nx); for (int i = 0; i < nx; ++i) { g.x[i] = x0 + i * g.dx; } return g; }

这段代码只做一件事：把一维区间均匀离散成节点坐标。nx是节点数，dx由区间长度除以nx-1得到，x向量用来赋初值、算解析解误差和输出。参数改的时候注意，nx增加会让dx变小，如果时间步dt不变，显式格式的稳定数会迅速变大，后面第 4 章会专门算这个匹配关系。

注意：nx是节点数，不是单元数。课程实验里常有人把nx当单元数，循环写成i < nx，结果最后一个节点落在区间外，边界条件怎么调都不对。

2.2 系数矩阵怎么存：稠密、三对角与按行压缩的取舍

一维隐式格式会产生三对角系统，二维五点差分会产生稀疏带状系统。课程实验的数据规模通常不大，稠密矩阵写起来直观，但容易把边界贡献塞错位置；三对角用 Thomas 算法是 (O(n))，代码短、调试快；二维如果不想手写迭代，可以用按行压缩存稀疏矩阵，但对新手来说索引更容易翻车。选型建议是：一维隐式优先三对角，二维小规模先用迭代法，规模上来再考虑稀疏存储。

三对角系统写成 (a_i u_{i-1}+b_i u_i+c_i u_{i+1}=d_i)，其中a是下对角，b是主对角，c是上对角，d是右端项。边界条件不要留在矩阵里，统一移到右端项，这样 Thomas 算法的首尾处理简单很多。

#include <vector> #include <cmath> #include <stdexcept> // Thomas 算法：求解三对角系统 a_i u_{i-1} + b_i u_i + c_i u_{i+1} = d_i std::vector<double> solve_tridiag(const std::vector<double>& a, const std::vector<double>& b, const std::vector<double>& c, const std::vector<double>& d) { int n = static_cast<int>(b.size()); std::vector<double> cp(n), dp(n), x(n); cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]; for (int i = 1; i < n; ++i) { double m = b[i] - a[i] * cp[i - 1]; if (std::abs(m) < 1e-14) { throw std::runtime_error("zero pivot in tridiagonal solve"); } cp[i] = (i + 1 < n) ? c[i] / m : 0.0; dp[i] = (d[i] - a[i] * dp[i - 1]) / m; } x[n - 1] = dp[n - 1]; for (int i = n - 2; i >= 0; --i) { x[i] = dp[i] - cp[i] * x[i + 1]; } return x; }

这段代码是追赶法，也就是 Thomas 算法。a[0]和c[n-1]应该置 0，因为它们对应边界外的未知量；b不能出现零主元，否则说明离散系数或边界处理有问题。参数d里要包含边界值移项后的贡献，比如左边界固定为u0，第一行右端项就要减去a[0]*u0。如果这里漏了，解出来的内部点看着平滑，边界附近会一直差一截。

提示：隐式格式的边界贡献要移到右端项d，不要留在矩阵里。很多“收敛阶只有一阶”的问题，源头就是边界项处理不一致。

2.3 时间推进的稳定性参数：显式、隐式与 CFL 条件

显式 FTCS 的稳定条件是 (r=\alpha\Delta t/\Delta x^2\le 0.5)（一维）。二维如果dx=dy，稳定条件是 (4r\le 1)，不是 (2r\le 1)。对流项另算 CFL 数：(|a|\Delta t/\Delta x\le 1)。显式的优点是每步只扫一遍数组，代码短，容易验证收敛阶；缺点是对dt很敏感，长时间积分或者扩散系数大时步长会被压得很小。隐式没有这个限制，但每步要解线性系统，代码量和调试成本都上去。

选型可以按实验目的来：只想快速看温度分布趋势，用显式；要跑长时间、边界复杂、或者扩散系数很大，用隐式或 Crank-Nicolson；课程实验要展示收敛阶，显式更容易把空间误差单独拆出来。下面这段是一维显式推进的核心循环，边界在外部函数里更新。

#include <vector> // 一维热传导 u_t = alpha u_xx，显式 FTCS 推进一步 void step_ftcs(std::vector<double>& u, double alpha, double dt, double dx) { int n = static_cast<int>(u.size()); std::vector<double> un = u; // 保存旧时间层 double r = alpha * dt / (dx * dx); for (int i = 1; i < n - 1; ++i) { u[i] = un[i] + r * (un[i - 1] - 2.0 * un[i] + un[i + 1]); } // u[0] 和 u[n-1] 由边界条件函数在调用后更新 }

这段代码里un是旧时间层，必须保留，否则就地更新会让u[i-1]变成新值，格式就不是 FTCS 了。r是稳定数，调用前最好打印出来；alpha是热扩散系数，dt是时间步，dx是空间步。边界点没有在循环里更新，是因为边界条件可能是 Dirichlet、Neumann 或周期，统一放到外层处理更清楚。显式格式不是“跑得动就稳定”，r接近 0.5 时误差可能先小后炸，血泪经验是先把r压到 0.4 以下再调其他参数。

注意：二维显式热传导用dx=dy时稳定条件是 (4\alpha\Delta t/\Delta x^2\le 1)。把一维的 0.5 直接套到二维，程序迟早给你颜色看。

3. 把源码跑通：编译、边界条件与收敛阶验证的三步闭环

拿到代码包后不要先改物理参数，先跑最小算例。三步是编译、边界、收敛阶，每步留一个可对比的输出。课程实验里最常见的翻车不是公式错，而是编译选项、边界覆盖顺序和误差范数不一致。下面给具体命令和代码，按顺序走一遍，能少掉很多“看起来对、一算就错”的时间。

3.1 目录结构与编译命令：g++/CMake 两种常见组织方式

单文件 Demo 适合先验证格式，多文件工程适合课程大作业。常见目录是src/放实现，include/放头文件，tests/放最小算例，data/放输出。编译时打开-Wall -Wextra，很多索引和类型问题编译器会直接提醒。优化选项-O2在跑大网格时有用，但调试阶段可以先不开，方便定位。

# 单文件直接编译，适合先跑最小算例 g++ -std=c++11 -O2 -Wall -Wextra main.cpp -o pde_demo # 多文件工程用 CMake 时，常见流程 mkdir -p build cd build cmake .. -DCMAKE_BUILD_TYPE=Release cmake --build . -j

这两条路径对应两种组织方式。单文件编译适合验证一个格式，命令短，出错信息直接；CMake 适合多文件工程，CMAKE_BUILD_TYPE=Release会打开优化，-j用多核并行编译。参数改的时候注意，调试阶段把-O2换成-O0 -g，方便用调试器看数组；发布或跑大网格再切回Release。如果出现链接错误，先看头文件声明和源文件是不是都进了 target，不要一上来怀疑算法。

提示：先编译最小main，不要一上来全工程。最小算例能跑通，再逐个加边界和输出，排查成本低很多。

3.2 边界条件与初值怎么写：Dirichlet、Neumann 与周期边界的实现差异

Dirichlet 边界直接固定端点值；Neumann 零通量用一阶差分是u[0]=u[1]，用二阶虚拟节点是u[0]=u[2]再配合边界方程；周期边界让首尾互为邻居。三种边界的循环范围不一样，周期边界内部点会跨过首尾，不能直接套 Dirichlet 的i=1..n-2。初值要和边界一致，否则第一步就会在边界附近产生一个虚假的尖峰。

#include <vector> // Dirichlet：边界值固定 void apply_dirichlet(std::vector<double>& u, double left, double right) { u.front() = left; u.back() = right; } // Neumann 零通量：一阶近似，u[0] = u[1] void apply_neumann_first_order(std::vector<double>& u) { u.front() = u[1]; u.back() = u[u.size() - 2]; } // 周期边界：首尾互为邻居 void apply_periodic(std::vector<double>& u) { u.front() = u[u.size() - 2]; u.back() = u[1]; }

这三个函数要在每次时间推进后调用，调用顺序不能乱。Dirichlet 最简单，但初值和边界不一致时，第一步会出现不连续；Neumann 一阶近似会把整体精度拉到一阶，做收敛阶验证时如果发现阶数只有 1，先检查这里；周期边界要求内部循环也按周期索引取邻居，不能只改首尾两个点。参数left、right是边界值，u是当前时间层。实际代码里建议把边界类型做成枚举，避免运行时靠if猜。

注意：Neumann 边界不要直接令u[0]=u[1]后还按 Dirichlet 的循环范围算，边界附近会少一项。周期边界同理，内部循环的索引要跟着改。

3.3 用解析解做收敛阶验证：误差表、网格加密与阶数拟合

收敛阶验证是判断差分程序是否可信的硬指标。选一个解析解，比如一维热传导取 (u(x,t)=e^{-\pi^2\alpha t}\sin(\pi x))，初值和边界都从这个解来。网格加密nx=21,41,81,161，固定最终时间，dt也按比例缩小，算离散 L2 误差 (\sqrt{\sum e_i^2\Delta x})。如果格式是二阶，误差应该大约按 (dx^2) 下降，拟合出来的观测阶接近 2。

import numpy as np # 把这里换成程序实际输出的网格数和 L2 误差 nx = np.array([21, 41, 81, 161], dtype=float) err = np.array([8.1e-3, 2.0e-3, 5.0e-4, 1.25e-4]) # 用 log(err) 对 log(dx) 做一次线性拟合，斜率取负就是观测阶 dx = 1.0 / (nx - 1) slope = np.polyfit(np.log(dx), np.log(err), 1)[0] print("observed order =", -slope)

这段脚本只做一件事：把误差随网格变化的趋势拟合成阶数。nx是节点数，dx是间距，err是对应误差。实际使用时把示例数字替换成程序输出，不要拿别人的表凑。拟合前先看误差是不是单调下降，如果某一级突然变大，先查边界和源项，不要先怀疑编译器。下面这张表是填表格式，观测阶按相邻两级误差比取对数得到。

误差表要固定最终物理时间，dt随dx同步缩小，否则时间离散误差会混进空间阶里。如果dt固定不变，网格越细，时间误差占比越大，拟合出来的阶数会偏低。另一个常见问题是误差范数少乘了dx，L2 范数变成 L∞ 趋势，看起来也在降，但阶数对不上。先把这两点对齐，再谈优化。

提示：如果误差不降，先查边界和源项，不要先怀疑编译器。差分程序里“编译器优化导致结果不对”的概率，远低于边界覆盖顺序写错。

4. 参数怎么设才不翻车：步长、库朗数与迭代容差

参数不是拍脑袋。显式看稳定数，隐式看迭代容差，输出看物理量。课程实验里最常见的参数错误是把dt设成固定值再改dx，或者把 SOR 的omega设成 1.9 但网格很小。下面按三类参数拆：空间步长和时间步长的匹配、隐式迭代的松弛因子、输出格式和切片。

4.1 空间步长与时间步长的匹配：从 CFL 数反推 dt

一维热传导显式稳定条件是 (r=\alpha\Delta t/\Delta x^2\le 0.5)，二维dx=dy时是 (4r\le 1)，一维对流是 (|a|\Delta t/\Delta x\le 1)。实际写代码时不要每次手算，封装一个函数，把安全系数留在参数里。安全系数一般取 0.8 到 0.9，留一点余量，别贴着临界值跑。

#include <cmath> // 一维显式热传导的稳定步长：alpha*dt/dx^2 <= 0.5 double stable_dt_ftcs(double alpha, double dx, double safety = 0.9) { return safety * 0.5 * dx * dx / alpha; } // 一维对流 CFL：|a|*dt/dx <= 1 double stable_dt_advection(double a, double dx, double safety = 0.9) { return safety * dx / std::abs(a); }

这两个函数分别对应扩散项和对流项。alpha是热扩散系数，a是对流速度，dx是空间步，safety是安全系数。二维热传导不要直接调一维函数，按 (4\alpha\Delta t/\Delta x^2\le 1) 自己写一行，或者把维度作为参数传进去。参数改的时候先算dt再改nx，不要反过来；先改网格再补时间步，很容易把稳定数顶到临界值以上。

注意：二维热传导用dx=dy时稳定条件是 (4r\le 1)，不是 (2r\le 1)。把一维条件直接搬过去，前面几步可能正常，后面会突然炸。

4.2 隐式格式的迭代求解：Jacobi、Gauss-Seidel 与 SOR 的松弛因子

隐式格式每步要解线性系统。一维三对角可以直接用 Thomas 算法，二维五点差分常用 Jacobi、Gauss-Seidel 或 SOR。Jacobi 每步用旧值更新，容易并行但收敛慢；Gauss-Seidel 就地更新，收敛快一些；SOR 在 Gauss-Seidel 上加松弛因子omega，omega=1就是 Gauss-Seidel，1<omega<2是超松弛。最优omega和网格规模有关，常见经验是 (2/(1+\sin(\pi/N)))，但课程实验里先跑omega=1看收敛，再逐步调到 1.5 左右比较稳。

#include <vector> #include <cmath> #include <algorithm> // 二维泊松方程 -Laplacian u = f，五点差分，dx=dy=h void sor_poisson(std::vector<std::vector<double>>& u, const std::vector<std::vector<double>>& f, double h, double omega, int max_iter, double tol) { int ny = static_cast<int>(u.size()); int nx = static_cast<int>(u[0].size()); for (int it = 0; it < max_iter; ++it) { double max_change = 0.0; for (int j = 1; j < ny - 1; ++j) { for (int i = 1; i < nx - 1; ++i) { double gs = 0.25 * (u[j][i - 1] + u[j][i + 1] + u[j - 1][i] + u[j + 1][i] + h * h * f[j][i]); double old = u[j][i]; u[j][i] = old + omega * (gs - old); max_change = std::max(max_change, std::abs(u[j][i] - old)); } } if (max_change < tol) break; } }

这段代码用 SOR 迭代解二维泊松方程。u是包含边界点的场，f是源项，h是网格间距，omega是松弛因子，max_iter是迭代上限，tol是相邻两次迭代的最大变化量。omega=1时退化成 Gauss-Seidel；omega取 1.5 到 1.8 常见，但网格很小的时候过大反而震荡。tol不要设得太小，否则白跑；也不要只看迭代次数，最好同时记录残差 (||f+\Delta u||)，残差曲线不降说明边界或源项有问题。

提示：tol是迭代变化量，不是残差。严格验证时看残差比看变化量可靠。残差平台不降，先查边界条件是不是每步都更新了。

4.3 输出与后处理：场数据格式、切片与可视化工具衔接

输出格式按用途选：CSV 简单，适合小网格和快速检查；VTK 结构化点适合 ParaView 直接看场；二进制格式快，但需要自己写读接口。课程实验里最稳的是先写文本 VTK，小网格用文本编辑器打开检查维度、间距和数据顺序，确认无误再换二进制。切片输出时注意坐标和索引的对应关系，不要只输出值不输出坐标。

#include <fstream> #include <string> #include <vector> // 输出二维结构化点 VTK，u[j][i] 中 j 是 y 行，i 是 x 列 void write <p> <a href="https://download.csdn.net/download/weixin_43652753/12377529" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>