简介《科学与工程计算》熊春光版期末复习笔记面向自动化、理工科本科生及备考者用于期末系统复习与预习。笔记围绕微分方程数值求解展开系统梳理初值问题与边界条件、有限差分法含节点、步长与中心差商逼近、内部截断误差、基函数法、拉格朗日插值及龙格现象、三次样条插值、配置法等核心内容并给出将常微分方程改写为矩阵形式及引入边界条件消除虚拟未知量的完整流程有助于构建从理论到算法的完整知识框架。资源为1个PDF文件大小约20.04MB排版紧凑重点公式标注清晰适合打印或平板批注使用。目前已有2767人学习浏览是期末冲刺阶段的高频下载资料。通过这份笔记可快速掌握两点边值问题的数值解法理解不同离散格式的构造思路并领会三次样条基函数对称平移的构造特性对控制系统数学建模和科学计算能力提升均有直接帮助。1. 网格比 r 的生死线为什么同一个热传导方程会算出三种命运同一道热传导方程把时间导数分别换成向前差商、向后差商和中心差商会得到显式、隐式、Richardson 三种差分格式。三者看起来都合理命运却截然不同显式格式在网格比 r0.5 时振荡发散隐式格式任意网格比都稳定Richardson 格式则无条件不稳定。这个反直觉现象说明差分格式的稳定性不是实现细节而是格式本身的数学性质。理解它就理解了《科学与工程计算》这门课的主线——从边值问题的差分离散到抛物型方程的三套格式再到 Von Neumann 稳定性判定和 Lax 等价定理。自动化专业的控制系统仿真、状态估计都要走到“连续模型离散化”这一步这个判断标准早晚用得上。2. 有限差分法中心差商如何同时服务边值问题和初值问题2.1 定解条件初值、边值与两点边值问题笔记开篇把定解条件分成两类与时间有关的初始条件构成初值问题与空间有关的边界条件构成边值问题。对二阶常微分方程数值求解需要定解区间两端的边界条件这就是“两点边值问题”的由来。这里的“两点”不是两个孤立的点而是区间端点的配对左端点给一个条件右端点给一个条件解才能被唯一确定。为什么必须两端都给二阶方程的通解带两个任意常数一个边界条件只能确定其中一个只有两端信息齐全离散后的代数方程组才会是方阵、有唯一解。这个看起来简单的结论在下面组装矩阵时最容易翻车——程序里常常只处理内部节点把边界行漏掉结果矩阵奇异或者解在端点处翘起来。2.2 差商近似从导数定义到内部截断误差把区间等距划分节点 x_i x_0 i hh 是步长。用适当差商近似导数中心差商是这套方法的核心武器y(x_i) ≈ (y_{i1} - y_{i-1}) / (2h)y(x_i) ≈ (y_{i1} - 2y_i y_{i-1}) / h²代入二阶常微分方程 y p(x)y q(x)y f(x)得到每个内部节点上的三点格式(y_{i1} - 2y_i y_{i-1})/h² p_i (y_{i1} - y_{i-1})/(2h) q_i y_i f_i这里的每一步都不精确差商代替导数丢掉的是 Taylor 展开的高阶项。笔记把这个局部误差记为 R_i称为内部截断误差。它的量级是 O(h²)来自二阶导数的四阶项和一阶导数的二阶项。h 减半截断误差大约降到四分之一这就是加密网格管用的数学原因。# 组装两点边值问题的三点格式系数示意 import numpy as np n 10 h 1.0 / n A np.zeros((n - 1, n - 1)) F np.zeros(n - 1) for i in range(1, n): # 内部节点 x_i i*h x i * h a 1/h**2 - p(x)/(2*h) # y_{i-1} 的系数 b q(x) - 2/h**2 # y_i 的系数 c 1/h**2 p(x)/(2*h) # y_{i1} 的系数 A[i-1, i-1] b if i 1: A[i-1, i-2] a if i n-1: A[i-1, i] c F[i-1] f(x)这个循环把每个内部节点的三点方程写入矩阵行。a、b、c 三项里p/(2h) 来自一阶导数的中心差商q - 2/h² 来自二阶导数合并如果 p(x) 数值很大而 h 很小p/(2h) 会明显放大这是对流占优问题里差分解振荡的根源之一。边界条件要单独处理首末两行把 u(x_0)、u(x_n) 的已知值搬到右端项矩阵才闭合。2.3 三对角方程组从差分格式到可求解的代数结构内部节点 i1,2,...,N-1 各给出一个方程加上两端边界条件补足首尾行得到的是一张三对角矩阵。非零元只落在主对角线、上对角线和下对角线上。这个结构不是巧合三点格式里每个方程只涉及相邻节点离散化的“局部性”直接决定了矩阵的带宽。名称表达式作用节点 x_ix_0 i h空间离散位置步长 h(b-a)/N控制截断误差量级网格比 raτ/h²决定显式格式稳定性内部截断误差 R_iO(h²)差分逼近微分方程的局部误差三对角系统可以用追赶法Thomas 算法求解复杂度 O(N)而一般高斯消去是 O(N³)。对这里的问题追赶法的优势不只是快主对角占优条件下它不需要选主元数值稳定性也有保证。笔记里特意把格式写成矩阵形式 A y F目的就是让人看到这个代数结构而不是一头扎进循环里去调系数。常见误区是“边界条件是在改方程”。实际上边界条件是补方程它给首末行提供独立信息。如果把边界行删掉方程组欠定如果用一个与内部格式无关的方式去写边界行追赶法的递推关系会被破坏。正确做法是让边界行和内部行一起组装成同一个三对角结构只是右端项里带上已知边界值。2.4 配置法与最小二乘法当差商法不够用时的两个替代并非所有问题都方便直接打差分。笔记在边值问题这一章里还给了基函数法路线设 u(x) ≈ Σ a_k φ_k(x)把参数 a_k 当未知量。这样就不需要每个节点都写差分方程而是让残差在某种意义下变小。配置法要求残差在选定的配置点上严格等于零最小二乘法要求残差平方的积分最小。笔记给了一个具体算例解 u u f边界条件 u(0)u(1)0设 u(x)a₁xa₂x²。两个边界条件已经消掉一个自由度剩一个未知数。配置法取区间中点 x1/2 作配置点残差方程直接解出 a₁最小二乘法则对 ∫(uu-f)² dx 关于 a₁ 求导等于零。两个方法在这个小算例上都能算对但思维方式完全不同配置法只盯有限个点最小二乘盯整个区间上的能量。这个例子还点出一个工程事实配置点数量必须等于剩余未知量数量。如果用三次 B 样条做基函数未知量动辄几十个配置点也要同样多而多项式逼近参数少一个中点就够。选基函数和选配置点其实是同一个决策的两面。3. 基函数法与样条逼近Runge 现象之后为什么是三次 B 样条接棒3.1 基函数法的一般框架基函数法的核心是把解限制在一个有限维函数空间里u(x) ≈ Σ a_k φ_k(x)φ_k 是基函数a_k 是待定系数。笔记把基函数分为三类多项式插值型φ_k x^k即 Lagrange 插值、B 样条插值型、Fourier 逼近型。选不同基函数离散格式的形态完全不同——稠密还是稀疏、二阶光滑还是任意光滑、有没有 Runge 振荡全都由 φ_k 决定。3.2 Runge 现象高次多项式插值的反直觉崩塌笔记明确写高次多项式会出现 Runge 现象。这在等距节点、多项式次数升高时尤其明显插值误差在区间端点附近不降反升出现波纹状振荡。原因要从 Lebesgue 常数理解等距节点的高次插值算子范数随节点数指数增长端点附近的基函数出现“过冲”微小的数据扰动被放大成大幅振荡。这个现象直接宣告了“次数越高越精确”的直觉破产。实用方案是分段低次插值而三次 B 样条恰好是分段三次、整体二阶连续可导的折中选择。3.3 三次 B 样条对称、平移与虚拟未知量的消除二阶常微分方程的解需要二阶连续可导三次 B 样条刚好满足。笔记强调的两个几何性质很关键每个基函数图形对称最大值统一为 1不同节点的基函数可以由同一个模板平移得到。翻译成实现语言就是程序里只需要一份局部基函数表达式所有节点的基函数都是它搬过去的结果不用逐节点积出系数。局部支撑性质带来一个间接后果x_i 处的函数值只被相邻三个节点的样条函数影响代入微分方程后每个方程只包含 a_{i-1}、a_i、a_{i1} 三个系数。于是第一行方程会出现 a_{-1}最后一行出现 a_{N1}——这两个是区间外的“虚拟未知量”不是真正要解的变量。# 用边界条件消去虚拟未知量伪代码示意 # 边界 u(0)0 给出关系 a_{-1} 4*a_0 a_1 0 # 把 a_{-1} -4*a_0 - a_1 代回第一个内部方程 row0 assemble_row(i1) # 内含 a_{-1}, a_0, a_1 三项 A[0, 0] -4 * row0[a_{-1}] # 折算到 a_0 A[0, 1] -1 * row0[a_{-1}] # 折算到 a_1代码里实现的不是某个固定公式而是“把区间外基函数的贡献折回真实未知量”这一思想。边界条件的种类决定折回关系的形式第一类边界直接给函数值第二类边界给导数导数边界会多出一项样条导数的系数。处理完首末两行方程组才是闭合、可逆的三对角或近三对角系统。3.4 配置点怎么选三次 B 样条配置的矩阵组装配置点可以取节点 x_i也可以取相邻节点中点。取节点时每个配置点对应一个残差方程矩阵接近三对角取中点时带宽略增但某些问题的条件数会更好。笔记里在配置法一节特别点明未知量个数决定配置点个数B 样条未知量多所以需要大量配置点。系统尺寸从“一个方程、一个未知数”的多项式逼近变成“几十个方程、几十个未知数”的稀疏线性系统这就是为什么这一章会落在样条——它用更大的系统换来了稳定性和光滑性。4. 抛物型方程离散化显式、隐式与三层格式的分岔点在哪4.1 差分方法的标准步骤笔记把差分方法的标准步骤归纳为网格剖分、离散化建立差分格式、求解代数方程组、用插值从离散解恢复近似解。前两步决定格式性质后两步决定实现代价。对一维非齐次热传导方程 u_t a u_xx f(x,t)网格剖分给出空间步长 h 与时间步长 τ网点值记为 u_i^k ≈ u(x_i, t_k)。这里最核心的是引入符号 r aτ/h²称为网格比。它同时包含扩散系数 a、时间步长 τ、空间步长 h显式格式的稳定性、误差的传播速度都浓缩在这个无量纲数里。为什么是“网格比”这个组合而非单独的 τ 或 h因为在显式格式的迭代表达式里下一时刻的值是上一时刻三个邻居的加权平均权重由 r 唯一确定其他组合方式得不到这样干净的结构。4.2 最简显式格式时间用一阶向前差商空间用二阶中心差商得到最简显式格式也叫向前差分格式u^{k1}j r u^k{j1} (1 - 2r) u^k_j r u^k_{j-1} τ f^k_j它叫“显式”是因为 u^{k1}_j 只依赖 k 时刻的已知值不需要解方程组每层推进代价 O(N)。但注意递推系数当前点的权重是 1-2r两个邻居各 r。要使新值仍是旧值的凸组合必须有 1-2r ≥ 0即 r ≤ 0.5。一旦 r 超过 0.5当前点的权重变负下一时刻的值会“反向依赖”旧值空间上立刻出现锯齿。# 最简显式格式的时间推进向量化示意 r a * tau / h**2 if r 0.5: raise ValueError(f网格比 r{r:.3f} 0.5高频分量会被逐层放大) for k in range(int(T / tau)): u_new[1:-1] ( r * u[2:] (1 - 2*r) * u[1:-1] r * u[:-2] tau * f[1:-1] ) u_new[0], u_new[-1] left_bc, right_bc u, u_new u_new, u这个向量化写法里u[2:] 对应 j1u[:-2] 对应 j-1边界值直接赋值给首末点。稳定性检查放在循环外而不是等振荡出现后是最便宜的后悔药——一次 raise 的成本远低于排查一组发散数据。局部截断误差是 O(τ h²)时间方向一阶、空间方向二阶。这意味着就算 r 满足条件时间步长 τ 也不能取得太大否则误差照样被时间方向的一阶项主导。4.3 最简隐式格式时间向后差商代替向前差商得到最简隐式格式-r u^{k1}_{j-1} (1 2r) u^{k1}j - r u^{k1}{j1} u^k_j τ f^k_j每推进一步需要解三对角方程组单步代价高于显式。但它换来一个更值钱的属性稳定性与网格比无关。笔记的说法是“隐式格式稳定性与网格比无任何关系称格式无条件稳定或绝对稳定”。工程上这意味着时间步长可以取大省下的步数通常远超过解三对角矩阵的开销——对扩散问题尤其明显。注意局部截断误差仍然只到 O(τ h²)。隐式格式并不比显式更精确它只是更稳定。这是最常见的误解把“隐式”误当成“高阶”。实际它是一阶时间精度只是稳定域大允许你把 τ 拉到显式不敢碰的量级。4.4 Richardson 格式三层格式的更高阶尝试Richardson 格式把时间方向换成中心差商期望用对称性消掉一阶时间误差u^{k1}j u^{k-1}j 2r (u^k{j1} - 2u^k_j u^k{j-1})局部截断误差确实提升到 O(τ² h²)。但它是三层格式计算 u^{k1} 需要 u^k 和 u^{k-1} 两个历史层而显式、隐式都是两层格式。三层格式有两个立即后果第一第一层 u¹ 不能由该格式自己产生必须先用其他两层格式启动一步第二稳定性分析不能直接套两层格式的标量增长因子必须把相邻两层组合成状态向量再处理。笔记对 Richardson 格式的处理方式是先把它改写成等价的两层方程组再求传播矩阵的特征值。结论是特征方程的两个特征值乘积恒为 -1绝对值之积为 1必有一个特征值的模大于或等于 1对任何网格比都不满足 Von Neumann 条件。也就是说这个看起来时间精度更高的格式根本不适合实际计算。提示一个格式局部截断误差阶数再高若增长因子越界网格加密只会让数值解在逼近真解之前先被误差淹没。高阶与稳定是两件事必须分开验证。5. 稳定性分析避坑Von Neumann 条件怎么用才不会翻车5.1 增长因子把格式变成谐波分量再判断稳定性分析的思路是傅里叶分解误差在网格上是所有频率谐波的叠加格式对每个频率的放大倍数就是增长因子 G(ω)。对常系数线性格式把 u^n_j λ^n e^{iω jh} 代入齐次差分方程化简后得到的 G(ω) 是个只依赖 ωh 和 r 的标量或矩阵。Von Neumann 条件是存在常数 C使对任意实数 ω 和时间层 n都有 |G(ω)| ≤ 1 O(τ)。常系数情形下直接要求 |G(ω)| ≤ 1。为什么只讨论齐次方程因为线性格式把非齐次右端项视为逐层累加的扰动而每一层扰动又按同一增长因子传播。若齐次解不被放大带右端项的完整问题也不会发散。笔记里“稳定性的讨论可以归结为初值稳定”这句话指的就是这个约化过程。# 显式格式增长因子扫描写代码前就能完成的稳定性判断 import numpy as np def G_explicit(r, theta): return 1 - 4*r*np.sin(theta/2)**2 # theta ωh for r in (0.4, 0.5, 0.6): vals [abs(G_explicit(r, t)) for t in np.linspace(0, np.pi, 400)] flag 稳定 if max(vals) 1.0 else 不稳定 print(fr{r:.1f} max|G|{max(vals):.3f} {flag}) # 输出: r0.4 max|G|1.000 稳定 # r0.5 max|G|1.000 临界 # r0.6 max|G|1.400 不稳定这个扫描只需在 θ ∈ [0, π] 上取几百个点成本几乎为零。每拿到一个新差商组合都先跑一遍这种扫描再写推进循环比事后从数据里找振荡原因省太多时间。5.2 避坑一显式格式 r0.5 时看到的是振荡而不是误差现象用最简显式格式跑热传导r0.6数值解在空间方向上呈锯齿形振荡时间越往后振荡越剧烈最终溢出。原因增长因子 G(ω) 1 - 4r sin²(ωh/2)。当 r0.5 时高频分量sin² 接近 1对应的 G 小于 -1每个时间步振幅被乘以绝对值大于 1 的负数。相邻节点一个被放大、一个被缩小自然形成棋盘式交替振荡。这和差分方程里 1-2r 变负是同一件事的两种叙述。解决要么缩小 τ 使 r ≤ 0.5要么换隐式格式把约束交给三对角求解。显式格式不是不能用而是必须把它控制在稳定区域内使用。我的血泪教训是跑显式格式前先算 r而不是等振荡出现再翻车。5.3 避坑二三层格式不能拿两层格式的理论直接套现象把 Richardson 格式的谐波解代进去按照两层格式的方式求标量增长因子得出“可能稳定”的错误结论甚至在某些 r 下算出看似合格的 |G|。原因三层格式的真实状态是 (u^k, u^{k-1}) 的联合单变量代入丢掉了状态维数。必须先把三层格式改写成两层方程组v^k (u^k, u^{k-1})^T再求 2×2 传播矩阵的特征值。Richardson 格式的特征方程为 λ² 8r sin²(ωh/2) λ - 1 0两根之积恒为 -1故模长之积为 1至少一个模长 ≥ 1。这一点与 r 具体取值无关所以它是无条件不稳定。解决遇到任何三层格式先把状态扩维成两层方程组再对传播矩阵做 Von Neumann 检验。笔记里专门演示了这个降维过程目的就是别把两层格式的直觉带过来。5.4 避坑三样条配置时忘了消虚拟未知量现象B 样条配置法组装出的线性系统第一行和最后一行总是多出一个“不存在的变量”直接求解报奇异或边界附近数值异常。原因B 样条基函数在定解区间外仍然有支撑方程组第一行包含 a_{-1}最后一行包含 a_{N1}它们不是节点上的真实自由度。这是局部支撑函数的边界效应不是代码 bug而是模型闭合问题。解决用边界条件补两个方程把 a_{-1}、a_{N1} 分别用 a_0、a_1 和 a_{N-1}、a_N 表示后代回。第一类边界给函数值关系第二类边界给导数关系。做完这一步方程组才是方阵、可逆、可追赶的。第 3.3 节的伪代码演示的就是这个折回过程。6. 设计取舍为什么相容性加稳定性才等价于收敛性Lax 等价定理是这份笔记的思想收口对一个适定的线性初值问题如果差分格式相容那么格式收敛当且仅当格式稳定。三者的角色可以这样拆解相容性管局部——每个网格点上的差分方程在逼近原偏微分方程截断误差随步长趋零稳定性管全局——某一步产生的误差不会被后续步指数放大收敛性管最终——网格加密时数值解逼近真解。定理把它们分成两个独立验证的性质这是“等价”二字的真正含义。这个等价关系背后有一组设计取舍。显式格式每步代价低但必须锁死 r隐式格式单步贵但换来任意网格比都稳Richardson 格式局部精度高却因为全局不稳定被一票否决。换句话说一个格式是否可用不在局部误差多么漂亮而在误差能不能在时间推进中承得住。Lax 等价定理把“这个格式能不能用”从玄学变成可以在写代码前验证的判定局部看截断误差阶全局看增长因子模两个条件一凑结论就出来了。工程上做连续系统离散化也是一样的逻辑隐式欧拉和显式欧拉收敛阶相同但非线性刚度系统里显式欧拉根本推不动换隐式后大时间步就平稳跑完。我做模拟项目X的连续系统离散化时踩过这个坑当时只盯着差分公式的截断误差以为精度高就完事结果一上刚度矩阵就振荡。从那以后我每次拿到新格式都先写出增长因子并检查边界条件是否闭合然后才允许自己写推进循环。希望帮到你——把这份笔记从边值问题的三点格式一路推到 Richardson 格式的特征方程推导链在草稿纸上完整走一遍收获会比孤立背公式大得多。本文还有配套的精品资源点击获取 
长河君 · 2026/10/23 8:01:43 · 阅读 1.2k 有限差分法稳定性解析：从网格比到Von Neumann条件
简介《科学与工程计算》熊春光版期末复习笔记面向自动化、理工科本科生及备考者用于期末系统复习与预习。笔记围绕微分方程数值求解展开系统梳理初值问题与边界条件、有限差分法含节点、步长与中心差商逼近、内部截断误差、基函数法、拉格朗日插值及龙格现象、三次样条插值、配置法等核心内容并给出将常微分方程改写为矩阵形式及引入边界条件消除虚拟未知量的完整流程有助于构建从理论到算法的完整知识框架。资源为1个PDF文件大小约20.04MB排版紧凑重点公式标注清晰适合打印或平板批注使用。目前已有2767人学习浏览是期末冲刺阶段的高频下载资料。通过这份笔记可快速掌握两点边值问题的数值解法理解不同离散格式的构造思路并领会三次样条基函数对称平移的构造特性对控制系统数学建模和科学计算能力提升均有直接帮助。1. 网格比 r 的生死线为什么同一个热传导方程会算出三种命运同一道热传导方程把时间导数分别换成向前差商、向后差商和中心差商会得到显式、隐式、Richardson 三种差分格式。三者看起来都合理命运却截然不同显式格式在网格比 r0.5 时振荡发散隐式格式任意网格比都稳定Richardson 格式则无条件不稳定。这个反直觉现象说明差分格式的稳定性不是实现细节而是格式本身的数学性质。理解它就理解了《科学与工程计算》这门课的主线——从边值问题的差分离散到抛物型方程的三套格式再到 Von Neumann 稳定性判定和 Lax 等价定理。自动化专业的控制系统仿真、状态估计都要走到“连续模型离散化”这一步这个判断标准早晚用得上。2. 有限差分法中心差商如何同时服务边值问题和初值问题2.1 定解条件初值、边值与两点边值问题笔记开篇把定解条件分成两类与时间有关的初始条件构成初值问题与空间有关的边界条件构成边值问题。对二阶常微分方程数值求解需要定解区间两端的边界条件这就是“两点边值问题”的由来。这里的“两点”不是两个孤立的点而是区间端点的配对左端点给一个条件右端点给一个条件解才能被唯一确定。为什么必须两端都给二阶方程的通解带两个任意常数一个边界条件只能确定其中一个只有两端信息齐全离散后的代数方程组才会是方阵、有唯一解。这个看起来简单的结论在下面组装矩阵时最容易翻车——程序里常常只处理内部节点把边界行漏掉结果矩阵奇异或者解在端点处翘起来。2.2 差商近似从导数定义到内部截断误差把区间等距划分节点 x_i x_0 i hh 是步长。用适当差商近似导数中心差商是这套方法的核心武器y(x_i) ≈ (y_{i1} - y_{i-1}) / (2h)y(x_i) ≈ (y_{i1} - 2y_i y_{i-1}) / h²代入二阶常微分方程 y p(x)y q(x)y f(x)得到每个内部节点上的三点格式(y_{i1} - 2y_i y_{i-1})/h² p_i (y_{i1} - y_{i-1})/(2h) q_i y_i f_i这里的每一步都不精确差商代替导数丢掉的是 Taylor 展开的高阶项。笔记把这个局部误差记为 R_i称为内部截断误差。它的量级是 O(h²)来自二阶导数的四阶项和一阶导数的二阶项。h 减半截断误差大约降到四分之一这就是加密网格管用的数学原因。# 组装两点边值问题的三点格式系数示意 import numpy as np n 10 h 1.0 / n A np.zeros((n - 1, n - 1)) F np.zeros(n - 1) for i in range(1, n): # 内部节点 x_i i*h x i * h a 1/h**2 - p(x)/(2*h) # y_{i-1} 的系数 b q(x) - 2/h**2 # y_i 的系数 c 1/h**2 p(x)/(2*h) # y_{i1} 的系数 A[i-1, i-1] b if i 1: A[i-1, i-2] a if i n-1: A[i-1, i] c F[i-1] f(x)这个循环把每个内部节点的三点方程写入矩阵行。a、b、c 三项里p/(2h) 来自一阶导数的中心差商q - 2/h² 来自二阶导数合并如果 p(x) 数值很大而 h 很小p/(2h) 会明显放大这是对流占优问题里差分解振荡的根源之一。边界条件要单独处理首末两行把 u(x_0)、u(x_n) 的已知值搬到右端项矩阵才闭合。2.3 三对角方程组从差分格式到可求解的代数结构内部节点 i1,2,...,N-1 各给出一个方程加上两端边界条件补足首尾行得到的是一张三对角矩阵。非零元只落在主对角线、上对角线和下对角线上。这个结构不是巧合三点格式里每个方程只涉及相邻节点离散化的“局部性”直接决定了矩阵的带宽。名称表达式作用节点 x_ix_0 i h空间离散位置步长 h(b-a)/N控制截断误差量级网格比 raτ/h²决定显式格式稳定性内部截断误差 R_iO(h²)差分逼近微分方程的局部误差三对角系统可以用追赶法Thomas 算法求解复杂度 O(N)而一般高斯消去是 O(N³)。对这里的问题追赶法的优势不只是快主对角占优条件下它不需要选主元数值稳定性也有保证。笔记里特意把格式写成矩阵形式 A y F目的就是让人看到这个代数结构而不是一头扎进循环里去调系数。常见误区是“边界条件是在改方程”。实际上边界条件是补方程它给首末行提供独立信息。如果把边界行删掉方程组欠定如果用一个与内部格式无关的方式去写边界行追赶法的递推关系会被破坏。正确做法是让边界行和内部行一起组装成同一个三对角结构只是右端项里带上已知边界值。2.4 配置法与最小二乘法当差商法不够用时的两个替代并非所有问题都方便直接打差分。笔记在边值问题这一章里还给了基函数法路线设 u(x) ≈ Σ a_k φ_k(x)把参数 a_k 当未知量。这样就不需要每个节点都写差分方程而是让残差在某种意义下变小。配置法要求残差在选定的配置点上严格等于零最小二乘法要求残差平方的积分最小。笔记给了一个具体算例解 u u f边界条件 u(0)u(1)0设 u(x)a₁xa₂x²。两个边界条件已经消掉一个自由度剩一个未知数。配置法取区间中点 x1/2 作配置点残差方程直接解出 a₁最小二乘法则对 ∫(uu-f)² dx 关于 a₁ 求导等于零。两个方法在这个小算例上都能算对但思维方式完全不同配置法只盯有限个点最小二乘盯整个区间上的能量。这个例子还点出一个工程事实配置点数量必须等于剩余未知量数量。如果用三次 B 样条做基函数未知量动辄几十个配置点也要同样多而多项式逼近参数少一个中点就够。选基函数和选配置点其实是同一个决策的两面。3. 基函数法与样条逼近Runge 现象之后为什么是三次 B 样条接棒3.1 基函数法的一般框架基函数法的核心是把解限制在一个有限维函数空间里u(x) ≈ Σ a_k φ_k(x)φ_k 是基函数a_k 是待定系数。笔记把基函数分为三类多项式插值型φ_k x^k即 Lagrange 插值、B 样条插值型、Fourier 逼近型。选不同基函数离散格式的形态完全不同——稠密还是稀疏、二阶光滑还是任意光滑、有没有 Runge 振荡全都由 φ_k 决定。3.2 Runge 现象高次多项式插值的反直觉崩塌笔记明确写高次多项式会出现 Runge 现象。这在等距节点、多项式次数升高时尤其明显插值误差在区间端点附近不降反升出现波纹状振荡。原因要从 Lebesgue 常数理解等距节点的高次插值算子范数随节点数指数增长端点附近的基函数出现“过冲”微小的数据扰动被放大成大幅振荡。这个现象直接宣告了“次数越高越精确”的直觉破产。实用方案是分段低次插值而三次 B 样条恰好是分段三次、整体二阶连续可导的折中选择。3.3 三次 B 样条对称、平移与虚拟未知量的消除二阶常微分方程的解需要二阶连续可导三次 B 样条刚好满足。笔记强调的两个几何性质很关键每个基函数图形对称最大值统一为 1不同节点的基函数可以由同一个模板平移得到。翻译成实现语言就是程序里只需要一份局部基函数表达式所有节点的基函数都是它搬过去的结果不用逐节点积出系数。局部支撑性质带来一个间接后果x_i 处的函数值只被相邻三个节点的样条函数影响代入微分方程后每个方程只包含 a_{i-1}、a_i、a_{i1} 三个系数。于是第一行方程会出现 a_{-1}最后一行出现 a_{N1}——这两个是区间外的“虚拟未知量”不是真正要解的变量。# 用边界条件消去虚拟未知量伪代码示意 # 边界 u(0)0 给出关系 a_{-1} 4*a_0 a_1 0 # 把 a_{-1} -4*a_0 - a_1 代回第一个内部方程 row0 assemble_row(i1) # 内含 a_{-1}, a_0, a_1 三项 A[0, 0] -4 * row0[a_{-1}] # 折算到 a_0 A[0, 1] -1 * row0[a_{-1}] # 折算到 a_1代码里实现的不是某个固定公式而是“把区间外基函数的贡献折回真实未知量”这一思想。边界条件的种类决定折回关系的形式第一类边界直接给函数值第二类边界给导数导数边界会多出一项样条导数的系数。处理完首末两行方程组才是闭合、可逆的三对角或近三对角系统。3.4 配置点怎么选三次 B 样条配置的矩阵组装配置点可以取节点 x_i也可以取相邻节点中点。取节点时每个配置点对应一个残差方程矩阵接近三对角取中点时带宽略增但某些问题的条件数会更好。笔记里在配置法一节特别点明未知量个数决定配置点个数B 样条未知量多所以需要大量配置点。系统尺寸从“一个方程、一个未知数”的多项式逼近变成“几十个方程、几十个未知数”的稀疏线性系统这就是为什么这一章会落在样条——它用更大的系统换来了稳定性和光滑性。4. 抛物型方程离散化显式、隐式与三层格式的分岔点在哪4.1 差分方法的标准步骤笔记把差分方法的标准步骤归纳为网格剖分、离散化建立差分格式、求解代数方程组、用插值从离散解恢复近似解。前两步决定格式性质后两步决定实现代价。对一维非齐次热传导方程 u_t a u_xx f(x,t)网格剖分给出空间步长 h 与时间步长 τ网点值记为 u_i^k ≈ u(x_i, t_k)。这里最核心的是引入符号 r aτ/h²称为网格比。它同时包含扩散系数 a、时间步长 τ、空间步长 h显式格式的稳定性、误差的传播速度都浓缩在这个无量纲数里。为什么是“网格比”这个组合而非单独的 τ 或 h因为在显式格式的迭代表达式里下一时刻的值是上一时刻三个邻居的加权平均权重由 r 唯一确定其他组合方式得不到这样干净的结构。4.2 最简显式格式时间用一阶向前差商空间用二阶中心差商得到最简显式格式也叫向前差分格式u^{k1}j r u^k{j1} (1 - 2r) u^k_j r u^k_{j-1} τ f^k_j它叫“显式”是因为 u^{k1}_j 只依赖 k 时刻的已知值不需要解方程组每层推进代价 O(N)。但注意递推系数当前点的权重是 1-2r两个邻居各 r。要使新值仍是旧值的凸组合必须有 1-2r ≥ 0即 r ≤ 0.5。一旦 r 超过 0.5当前点的权重变负下一时刻的值会“反向依赖”旧值空间上立刻出现锯齿。# 最简显式格式的时间推进向量化示意 r a * tau / h**2 if r 0.5: raise ValueError(f网格比 r{r:.3f} 0.5高频分量会被逐层放大) for k in range(int(T / tau)): u_new[1:-1] ( r * u[2:] (1 - 2*r) * u[1:-1] r * u[:-2] tau * f[1:-1] ) u_new[0], u_new[-1] left_bc, right_bc u, u_new u_new, u这个向量化写法里u[2:] 对应 j1u[:-2] 对应 j-1边界值直接赋值给首末点。稳定性检查放在循环外而不是等振荡出现后是最便宜的后悔药——一次 raise 的成本远低于排查一组发散数据。局部截断误差是 O(τ h²)时间方向一阶、空间方向二阶。这意味着就算 r 满足条件时间步长 τ 也不能取得太大否则误差照样被时间方向的一阶项主导。4.3 最简隐式格式时间向后差商代替向前差商得到最简隐式格式-r u^{k1}_{j-1} (1 2r) u^{k1}j - r u^{k1}{j1} u^k_j τ f^k_j每推进一步需要解三对角方程组单步代价高于显式。但它换来一个更值钱的属性稳定性与网格比无关。笔记的说法是“隐式格式稳定性与网格比无任何关系称格式无条件稳定或绝对稳定”。工程上这意味着时间步长可以取大省下的步数通常远超过解三对角矩阵的开销——对扩散问题尤其明显。注意局部截断误差仍然只到 O(τ h²)。隐式格式并不比显式更精确它只是更稳定。这是最常见的误解把“隐式”误当成“高阶”。实际它是一阶时间精度只是稳定域大允许你把 τ 拉到显式不敢碰的量级。4.4 Richardson 格式三层格式的更高阶尝试Richardson 格式把时间方向换成中心差商期望用对称性消掉一阶时间误差u^{k1}j u^{k-1}j 2r (u^k{j1} - 2u^k_j u^k{j-1})局部截断误差确实提升到 O(τ² h²)。但它是三层格式计算 u^{k1} 需要 u^k 和 u^{k-1} 两个历史层而显式、隐式都是两层格式。三层格式有两个立即后果第一第一层 u¹ 不能由该格式自己产生必须先用其他两层格式启动一步第二稳定性分析不能直接套两层格式的标量增长因子必须把相邻两层组合成状态向量再处理。笔记对 Richardson 格式的处理方式是先把它改写成等价的两层方程组再求传播矩阵的特征值。结论是特征方程的两个特征值乘积恒为 -1绝对值之积为 1必有一个特征值的模大于或等于 1对任何网格比都不满足 Von Neumann 条件。也就是说这个看起来时间精度更高的格式根本不适合实际计算。提示一个格式局部截断误差阶数再高若增长因子越界网格加密只会让数值解在逼近真解之前先被误差淹没。高阶与稳定是两件事必须分开验证。5. 稳定性分析避坑Von Neumann 条件怎么用才不会翻车5.1 增长因子把格式变成谐波分量再判断稳定性分析的思路是傅里叶分解误差在网格上是所有频率谐波的叠加格式对每个频率的放大倍数就是增长因子 G(ω)。对常系数线性格式把 u^n_j λ^n e^{iω jh} 代入齐次差分方程化简后得到的 G(ω) 是个只依赖 ωh 和 r 的标量或矩阵。Von Neumann 条件是存在常数 C使对任意实数 ω 和时间层 n都有 |G(ω)| ≤ 1 O(τ)。常系数情形下直接要求 |G(ω)| ≤ 1。为什么只讨论齐次方程因为线性格式把非齐次右端项视为逐层累加的扰动而每一层扰动又按同一增长因子传播。若齐次解不被放大带右端项的完整问题也不会发散。笔记里“稳定性的讨论可以归结为初值稳定”这句话指的就是这个约化过程。# 显式格式增长因子扫描写代码前就能完成的稳定性判断 import numpy as np def G_explicit(r, theta): return 1 - 4*r*np.sin(theta/2)**2 # theta ωh for r in (0.4, 0.5, 0.6): vals [abs(G_explicit(r, t)) for t in np.linspace(0, np.pi, 400)] flag 稳定 if max(vals) 1.0 else 不稳定 print(fr{r:.1f} max|G|{max(vals):.3f} {flag}) # 输出: r0.4 max|G|1.000 稳定 # r0.5 max|G|1.000 临界 # r0.6 max|G|1.400 不稳定这个扫描只需在 θ ∈ [0, π] 上取几百个点成本几乎为零。每拿到一个新差商组合都先跑一遍这种扫描再写推进循环比事后从数据里找振荡原因省太多时间。5.2 避坑一显式格式 r0.5 时看到的是振荡而不是误差现象用最简显式格式跑热传导r0.6数值解在空间方向上呈锯齿形振荡时间越往后振荡越剧烈最终溢出。原因增长因子 G(ω) 1 - 4r sin²(ωh/2)。当 r0.5 时高频分量sin² 接近 1对应的 G 小于 -1每个时间步振幅被乘以绝对值大于 1 的负数。相邻节点一个被放大、一个被缩小自然形成棋盘式交替振荡。这和差分方程里 1-2r 变负是同一件事的两种叙述。解决要么缩小 τ 使 r ≤ 0.5要么换隐式格式把约束交给三对角求解。显式格式不是不能用而是必须把它控制在稳定区域内使用。我的血泪教训是跑显式格式前先算 r而不是等振荡出现再翻车。5.3 避坑二三层格式不能拿两层格式的理论直接套现象把 Richardson 格式的谐波解代进去按照两层格式的方式求标量增长因子得出“可能稳定”的错误结论甚至在某些 r 下算出看似合格的 |G|。原因三层格式的真实状态是 (u^k, u^{k-1}) 的联合单变量代入丢掉了状态维数。必须先把三层格式改写成两层方程组v^k (u^k, u^{k-1})^T再求 2×2 传播矩阵的特征值。Richardson 格式的特征方程为 λ² 8r sin²(ωh/2) λ - 1 0两根之积恒为 -1故模长之积为 1至少一个模长 ≥ 1。这一点与 r 具体取值无关所以它是无条件不稳定。解决遇到任何三层格式先把状态扩维成两层方程组再对传播矩阵做 Von Neumann 检验。笔记里专门演示了这个降维过程目的就是别把两层格式的直觉带过来。5.4 避坑三样条配置时忘了消虚拟未知量现象B 样条配置法组装出的线性系统第一行和最后一行总是多出一个“不存在的变量”直接求解报奇异或边界附近数值异常。原因B 样条基函数在定解区间外仍然有支撑方程组第一行包含 a_{-1}最后一行包含 a_{N1}它们不是节点上的真实自由度。这是局部支撑函数的边界效应不是代码 bug而是模型闭合问题。解决用边界条件补两个方程把 a_{-1}、a_{N1} 分别用 a_0、a_1 和 a_{N-1}、a_N 表示后代回。第一类边界给函数值关系第二类边界给导数关系。做完这一步方程组才是方阵、可逆、可追赶的。第 3.3 节的伪代码演示的就是这个折回过程。6. 设计取舍为什么相容性加稳定性才等价于收敛性Lax 等价定理是这份笔记的思想收口对一个适定的线性初值问题如果差分格式相容那么格式收敛当且仅当格式稳定。三者的角色可以这样拆解相容性管局部——每个网格点上的差分方程在逼近原偏微分方程截断误差随步长趋零稳定性管全局——某一步产生的误差不会被后续步指数放大收敛性管最终——网格加密时数值解逼近真解。定理把它们分成两个独立验证的性质这是“等价”二字的真正含义。这个等价关系背后有一组设计取舍。显式格式每步代价低但必须锁死 r隐式格式单步贵但换来任意网格比都稳Richardson 格式局部精度高却因为全局不稳定被一票否决。换句话说一个格式是否可用不在局部误差多么漂亮而在误差能不能在时间推进中承得住。Lax 等价定理把“这个格式能不能用”从玄学变成可以在写代码前验证的判定局部看截断误差阶全局看增长因子模两个条件一凑结论就出来了。工程上做连续系统离散化也是一样的逻辑隐式欧拉和显式欧拉收敛阶相同但非线性刚度系统里显式欧拉根本推不动换隐式后大时间步就平稳跑完。我做模拟项目X的连续系统离散化时踩过这个坑当时只盯着差分公式的截断误差以为精度高就完事结果一上刚度矩阵就振荡。从那以后我每次拿到新格式都先写出增长因子并检查边界条件是否闭合然后才允许自己写推进循环。希望帮到你——把这份笔记从边值问题的三点格式一路推到 Richardson 格式的特征方程推导链在草稿纸上完整走一遍收获会比孤立背公式大得多。本文还有配套的精品资源点击获取