简介偏微分方程数值解差分方法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^{n1}u_i^nr(u_{i-1}^n-2u_i^nu_{i1}^n))其中 (r\alpha\Delta t/\Delta x^2)。二维泊松方程 (-\Delta uf) 用五点差分内部点满足 ((4u_{i,j}-u_{i-1,j}-u_{i1,j}-u_{i,j-1}-u_{i,j1})/h^2f_{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::vectordouble 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_castdouble(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))代码短、调试快二维如果不想手写迭代可以用按行压缩存稀疏矩阵但对新手来说索引更容易翻车。选型建议是一维隐式优先三对角二维小规模先用迭代法规模上来再考虑稀疏存储。存储方式适用场景时间复杂度主要风险稠密二维数组教学演示、小规模(O(n^3)) 直接解边界项和内部项混在一起三对角数组一维隐式、Crank-Nicolson(O(n))首尾对角元没置零按行压缩 CSR二维大网格、稀疏系统依赖求解器行偏移和列索引对不上三对角系统写成 (a_i u_{i-1}b_i u_ic_i u_{i1}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_{i1} d_i std::vectordouble solve_tridiag(const std::vectordouble a, const std::vectordouble b, const std::vectordouble c, const std::vectordouble d) { int n static_castint(b.size()); std::vectordouble 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)一维。二维如果dxdy稳定条件是 (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::vectordouble u, double alpha, double dt, double dx) { int n static_castint(u.size()); std::vectordouble 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 以下再调其他参数。注意二维显式热传导用dxdy时稳定条件是 (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 -stdc11 -O2 -Wall -Wextra main.cpp -o pde_demo # 多文件工程用 CMake 时常见流程 mkdir -p build cd build cmake .. -DCMAKE_BUILD_TYPERelease cmake --build . -j这两条路径对应两种组织方式。单文件编译适合验证一个格式命令短出错信息直接CMake 适合多文件工程CMAKE_BUILD_TYPERelease会打开优化-j用多核并行编译。参数改的时候注意调试阶段把-O2换成-O0 -g方便用调试器看数组发布或跑大网格再切回Release。如果出现链接错误先看头文件声明和源文件是不是都进了 target不要一上来怀疑算法。提示先编译最小main不要一上来全工程。最小算例能跑通再逐个加边界和输出排查成本低很多。3.2 边界条件与初值怎么写Dirichlet、Neumann 与周期边界的实现差异Dirichlet 边界直接固定端点值Neumann 零通量用一阶差分是u[0]u[1]用二阶虚拟节点是u[0]u[2]再配合边界方程周期边界让首尾互为邻居。三种边界的循环范围不一样周期边界内部点会跨过首尾不能直接套 Dirichlet 的i1..n-2。初值要和边界一致否则第一步就会在边界附近产生一个虚假的尖峰。#include vector // Dirichlet边界值固定 void apply_dirichlet(std::vectordouble u, double left, double right) { u.front() left; u.back() right; } // Neumann 零通量一阶近似u[0] u[1] void apply_neumann_first_order(std::vectordouble u) { u.front() u[1]; u.back() u[u.size() - 2]; } // 周期边界首尾互为邻居 void apply_periodic(std::vectordouble 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))初值和边界都从这个解来。网格加密nx21,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], dtypefloat) 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是对应误差。实际使用时把示例数字替换成程序输出不要拿别人的表凑。拟合前先看误差是不是单调下降如果某一级突然变大先查边界和源项不要先怀疑编译器。下面这张表是填表格式观测阶按相邻两级误差比取对数得到。nxdxL2 误差观测阶210.058.1e-3-410.0252.0e-32.02810.01255.0e-42.001610.006251.25e-42.00误差表要固定最终物理时间dt随dx同步缩小否则时间离散误差会混进空间阶里。如果dt固定不变网格越细时间误差占比越大拟合出来的阶数会偏低。另一个常见问题是误差范数少乘了dxL2 范数变成 L∞ 趋势看起来也在降但阶数对不上。先把这两点对齐再谈优化。提示如果误差不降先查边界和源项不要先怀疑编译器。差分程序里“编译器优化导致结果不对”的概率远低于边界覆盖顺序写错。4. 参数怎么设才不翻车步长、库朗数与迭代容差参数不是拍脑袋。显式看稳定数隐式看迭代容差输出看物理量。课程实验里最常见的参数错误是把dt设成固定值再改dx或者把 SOR 的omega设成 1.9 但网格很小。下面按三类参数拆空间步长和时间步长的匹配、隐式迭代的松弛因子、输出格式和切片。4.1 空间步长与时间步长的匹配从 CFL 数反推 dt一维热传导显式稳定条件是 (r\alpha\Delta t/\Delta x^2\le 0.5)二维dxdy时是 (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不要反过来先改网格再补时间步很容易把稳定数顶到临界值以上。模型稳定条件推荐安全系数失败现象一维热传导显式(\alpha\Delta t/\Delta x^2\le 0.5)0.8 到 0.9几步后 NaN 或剧烈振荡二维热传导显式(4\alpha\Delta t/\Delta x^2\le 1)0.8 到 0.9场值出现棋盘格一维对流显式(a\Delta t/\Delta x\le 1)注意二维热传导用dxdy时稳定条件是 (4r\le 1)不是 (2r\le 1)。把一维条件直接搬过去前面几步可能正常后面会突然炸。4.2 隐式格式的迭代求解Jacobi、Gauss-Seidel 与 SOR 的松弛因子隐式格式每步要解线性系统。一维三对角可以直接用 Thomas 算法二维五点差分常用 Jacobi、Gauss-Seidel 或 SOR。Jacobi 每步用旧值更新容易并行但收敛慢Gauss-Seidel 就地更新收敛快一些SOR 在 Gauss-Seidel 上加松弛因子omegaomega1就是 Gauss-Seidel1omega2是超松弛。最优omega和网格规模有关常见经验是 (2/(1\sin(\pi/N)))但课程实验里先跑omega1看收敛再逐步调到 1.5 左右比较稳。#include vector #include cmath #include algorithm // 二维泊松方程 -Laplacian u f五点差分dxdyh void sor_poisson(std::vectorstd::vectordouble u, const std::vectorstd::vectordouble f, double h, double omega, int max_iter, double tol) { int ny static_castint(u.size()); int nx static_castint(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是相邻两次迭代的最大变化量。omega1时退化成 Gauss-Seidelomega取 1.5 到 1.8 常见但网格很小的时候过大反而震荡。tol不要设得太小否则白跑也不要只看迭代次数最好同时记录残差 (||f\Delta u||)残差曲线不降说明边界或源项有问题。提示tol是迭代变化量不是残差。严格验证时看残差比看变化量可靠。残差平台不降先查边界条件是不是每步都更新了。4.3 输出与后处理场数据格式、切片与可视化工具衔接输出格式按用途选CSV 简单适合小网格和快速检查VTK 结构化点适合 ParaView 直接看场二进制格式快但需要自己写读接口。课程实验里最稳的是先写文本 VTK小网格用文本编辑器打开检查维度、间距和数据顺序确认无误再换二进制。切片输出时注意坐标和索引的对应关系不要只输出值不输出坐标。#include fstream #include string #include vector // 输出二维结构化点 VTKu[j][i] 中 j 是 y 行i 是 x 列 void write p a hrefhttps://download.csdn.net/download/weixin_43652753/12377529 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
