简介：《科学与工程计算》熊春光版期末复习笔记，面向自动化、理工科本科生及备考者，用于期末系统复习与预习。笔记围绕微分方程数值求解展开，系统梳理初值问题与边界条件、有限差分法（含节点、步长与中心差商逼近）、内部截断误差、基函数法、拉格朗日插值及龙格现象、三次样条插值、配置法等核心内容，并给出将常微分方程改写为矩阵形式及引入边界条件消除虚拟未知量的完整流程，有助于构建从理论到算法的完整知识框架。资源为1个PDF文件，大小约20.04MB，排版紧凑，重点公式标注清晰，适合打印或平板批注使用。目前已有2767人学习浏览，是期末冲刺阶段的高频下载资料。通过这份笔记，可快速掌握两点边值问题的数值解法，理解不同离散格式的构造思路，并领会三次样条基函数对称平移的构造特性，对控制系统数学建模和科学计算能力提升均有直接帮助。
1. 网格比 r 的生死线：为什么同一个热传导方程会算出三种命运
同一道热传导方程，把时间导数分别换成向前差商、向后差商和中心差商，会得到显式、隐式、Richardson 三种差分格式。三者看起来都合理，命运却截然不同：显式格式在网格比 r>0.5 时振荡发散，隐式格式任意网格比都稳定，Richardson 格式则无条件不稳定。这个反直觉现象说明，差分格式的稳定性不是实现细节，而是格式本身的数学性质。理解它，就理解了《科学与工程计算》这门课的主线——从边值问题的差分离散，到抛物型方程的三套格式，再到 Von Neumann 稳定性判定和 Lax 等价定理。自动化专业的控制系统仿真、状态估计都要走到“连续模型离散化”这一步，这个判断标准早晚用得上。
2. 有限差分法：中心差商如何同时服务边值问题和初值问题
2.1 定解条件：初值、边值与两点边值问题
笔记开篇把定解条件分成两类：与时间有关的初始条件构成初值问题，与空间有关的边界条件构成边值问题。对二阶常微分方程，数值求解需要定解区间两端的边界条件，这就是“两点边值问题”的由来。这里的“两点”不是两个孤立的点，而是区间端点的配对：左端点给一个条件，右端点给一个条件，解才能被唯一确定。
为什么必须两端都给？二阶方程的通解带两个任意常数，一个边界条件只能确定其中一个；只有两端信息齐全，离散后的代数方程组才会是方阵、有唯一解。这个看起来简单的结论，在下面组装矩阵时最容易翻车——程序里常常只处理内部节点，把边界行漏掉，结果矩阵奇异或者解在端点处翘起来。
2.2 差商近似：从导数定义到内部截断误差
把区间等距划分，节点 x_i = x_0 + i h，h 是步长。用适当差商近似导数，中心差商是这套方法的核心武器：
y'(x_i) ≈ (y_{i+1} - y_{i-1}) / (2h)
y''(x_i) ≈ (y_{i+1} - 2y_i + y_{i-1}) / h²
代入二阶常微分方程 y'' + p(x)y' + q(x)y = f(x)，得到每个内部节点上的三点格式：
(y_{i+1} - 2y_i + y_{i-1})/h² + p_i (y_{i+1} - 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_{i+1} 的系数 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 三对角方程组：从差分格式到可求解的代数结构
内部节点 i=1,2,...,N-1 各给出一个方程，加上两端边界条件补足首尾行，得到的是一张三对角矩阵。非零元只落在主对角线、上对角线和下对角线上。这个结构不是巧合：三点格式里每个方程只涉及相邻节点，离散化的“局部性”直接决定了矩阵的带宽。
| 名称 | 表达式 | 作用 |
|---|---|---|
| 节点 x_i | x_0 + i h | 空间离散位置 |
| 步长 h | (b-a)/N | 控制截断误差量级 |
| 网格比 r | aτ/h² | 决定显式格式稳定性 |
| 内部截断误差 R_i | O(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₁x+a₂x²。两个边界条件已经消掉一个自由度，剩一个未知数。配置法取区间中点 x=1/2 作配置点，残差方程直接解出 a₁；最小二乘法则对 ∫(u''+u'-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_{i+1} 三个系数。于是第一行方程会出现 a_{-1}，最后一行出现 a_{N+1}——这两个是区间外的“虚拟未知量”，不是真正要解的变量。
# 用边界条件消去虚拟未知量（伪代码示意） # 边界 u(0)=0 给出关系 a_{-1} + 4*a_0 + a_1 = 0 # 把 a_{-1} = -4*a_0 - a_1 代回第一个内部方程 row0 = assemble_row(i=1) # 内含 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^{k+1}j = r u^k{j+1} + (1 - 2r) u^k_j + r u^k_{j-1} + τ f^k_j
它叫“显式”是因为 u^{k+1}_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:] 对应 j+1，u[:-2] 对应 j-1，边界值直接赋值给首末点。稳定性检查放在循环外而不是等振荡出现后，是最便宜的后悔药——一次 raise 的成本远低于排查一组发散数据。
局部截断误差是 O(τ + h²)：时间方向一阶、空间方向二阶。这意味着就算 r 满足条件，时间步长 τ 也不能取得太大，否则误差照样被时间方向的一阶项主导。
4.3 最简隐式格式
时间向后差商代替向前差商，得到最简隐式格式：
-r u^{k+1}_{j-1} + (1 + 2r) u^{k+1}j - r u^{k+1}{j+1} = u^k_j + τ f^k_j
每推进一步需要解三对角方程组，单步代价高于显式。但它换来一个更值钱的属性：稳定性与网格比无关。笔记的说法是“隐式格式稳定性与网格比无任何关系，称格式无条件稳定或绝对稳定”。工程上这意味着时间步长可以取大，省下的步数通常远超过解三对角矩阵的开销——对扩散问题尤其明显。
注意局部截断误差仍然只到 O(τ + h²)。隐式格式并不比显式更精确，它只是更稳定。这是最常见的误解：把“隐式”误当成“高阶”。实际它是一阶时间精度，只是稳定域大，允许你把 τ 拉到显式不敢碰的量级。
4.4 Richardson 格式：三层格式的更高阶尝试
Richardson 格式把时间方向换成中心差商，期望用对称性消掉一阶时间误差：
u^{k+1}j = u^{k-1}j + 2r (u^k{j+1} - 2u^k_j + u^k{j-1})
局部截断误差确实提升到 O(τ² + h²)。但它是三层格式：计算 u^{k+1} 需要 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(f"r={r:.1f} max|G|={max(vals):.3f} {flag}") # 输出: r=0.4 max|G|=1.000 稳定 # r=0.5 max|G|=1.000 临界 # r=0.6 max|G|=1.400 不稳定这个扫描只需在 θ ∈ [0, π] 上取几百个点，成本几乎为零。每拿到一个新差商组合，都先跑一遍这种扫描再写推进循环，比事后从数据里找振荡原因省太多时间。
5.2 避坑一：显式格式 r>0.5 时看到的是振荡而不是误差
现象：用最简显式格式跑热传导，r=0.6，数值解在空间方向上呈锯齿形振荡，时间越往后振荡越剧烈，最终溢出。
原因：增长因子 G(ω) = 1 - 4r sin²(ωh/2)。当 r>0.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_{N+1}，它们不是节点上的真实自由度。这是局部支撑函数的边界效应，不是代码 bug，而是模型闭合问题。
解决：用边界条件补两个方程，把 a_{-1}、a_{N+1} 分别用 a_0、a_1 和 a_{N-1}、a_N 表示后代回。第一类边界给函数值关系，第二类边界给导数关系。做完这一步，方程组才是方阵、可逆、可追赶的。第 3.3 节的伪代码演示的就是这个折回过程。
6. 设计取舍：为什么相容性加稳定性才等价于收敛性
Lax 等价定理是这份笔记的思想收口：对一个适定的线性初值问题，如果差分格式相容，那么格式收敛当且仅当格式稳定。三者的角色可以这样拆解：相容性管局部——每个网格点上的差分方程在逼近原偏微分方程，截断误差随步长趋零；稳定性管全局——某一步产生的误差不会被后续步指数放大；收敛性管最终——网格加密时数值解逼近真解。定理把它们分成两个独立验证的性质，这是“等价”二字的真正含义。
这个等价关系背后有一组设计取舍。显式格式每步代价低但必须锁死 r，隐式格式单步贵但换来任意网格比都稳，Richardson 格式局部精度高却因为全局不稳定被一票否决。换句话说，一个格式是否可用，不在局部误差多么漂亮，而在误差能不能在时间推进中承得住。Lax 等价定理把“这个格式能不能用”从玄学变成可以在写代码前验证的判定：局部看截断误差阶，全局看增长因子模，两个条件一凑，结论就出来了。工程上做连续系统离散化也是一样的逻辑：隐式欧拉和显式欧拉收敛阶相同，但非线性刚度系统里显式欧拉根本推不动，换隐式后大时间步就平稳跑完。
我做模拟项目X的连续系统离散化时踩过这个坑：当时只盯着差分公式的截断误差，以为精度高就完事，结果一上刚度矩阵就振荡。从那以后我每次拿到新格式，都先写出增长因子并检查边界条件是否闭合，然后才允许自己写推进循环。希望帮到你——把这份笔记从边值问题的三点格式一路推到 Richardson 格式的特征方程，推导链在草稿纸上完整走一遍，收获会比孤立背公式大得多。
本文还有配套的精品资源，点击获取
