简介《科学与工程计算》期末复习笔记依据熊春光版教材整理面向北理工及相关院校自动化、理工科专业学生适用于期末考前冲刺或课程预习梳理。资源包含1个PDF文件约20.04MB已有2767人学习。笔记覆盖两点边值问题、初值条件与边界条件、有限差分法及内部截断误差等核心内容结合节点、步长、差商选择等细节整理出中点差商减小误差、忽略高阶无穷小项列写代数方程组等关键思路同时梳理基函数法、Lagrange插值的Runge现象、三次样条插值函数的连续光滑性与基函数平移对称性质并补充配置法的基本思想。整份笔记将教材重要定义、经典数值方法与易混概念提炼为便于查阅的知识框架适合考前快速建立章节脉络、把握微分方程数值解法的推导思路与应用前提。1. 《科学与工程计算》期末复习笔记能顶什么用从定解条件到差分格式的全链路纲要自动化专业的学生学《科学与工程计算》最头疼的不是某个公式看不懂而是整门课的知识点散落在各个章节里平时作业会做、一到期末就发现前后对不上。这份手写整理的复习笔记把二阶常微分方程数值解、抛物型方程的三种古典差分格式、稳定性与收敛性分析串成了一条线核心脉络是先分清定解条件再选离散方法最后用 Fourier 分析判断格式能不能用。笔记里公式推导和矩阵形式都写得比较完整适合两类人——一类是考前突击需要快速过一遍数值方法主线的另一类是上课没听懂想对照推导细节补课的。尤其对自动化方向的学生来说控制系统中大量涉及微分方程的数值求解这份资料的价值不在背公式而在帮你建立“给定一个方程怎么选格式、怎么判断结果靠不靠谱”的完整思路。2. 二阶常微分方程数值解定解条件、中心差商与三对角矩阵的构造2.1 初值问题与边值条件先分清定解条件再选格式常微分方程数值解的第一步不是写差分格式而是先搞清楚问题属于哪一类。笔记第一章明确区分了初值问题与边值问题初值问题给定的是某个时间点的函数值和导数值信息集中在起点边值问题给定的是求解区间两个端点上的约束比如二阶常微分方程 $\frac{d^2y}{dx^2} p(x)\frac{dy}{dx} q(x)y f(x)$ 在区间 $[a,b]$ 上求解就需要 $xa$ 和 $xb$ 两端的边界条件。这个区分直接决定你用什么样的离散策略。初值问题一般从起点逐步推进用单步法或多步法边值问题则要同时考虑整个区间上的所有节点离散后得到的是一个代数方程组而不是一条递推关系。笔记里讨论的有限差分法和基函数法本质上都是针对边值问题做全区间离散这也是为什么最后会落到矩阵形式的求解上。2.2 中心差商近似与内部截断误差的推导在区间 $[a,b]$ 上做等距剖分记步长为 $h$节点为 $x_i a ih$。针对微分方程在节点 $x_i$ 处的取值用差商代替导数。笔记用的是中心差商$$ y(x_i) \approx \frac{y_{i1} - y_{i-1}}{2h}, \quad y(x_i) \approx \frac{y_{i1} - 2y_i y_{i-1}}{h^2} $$把这两个式子代回原方程就得到一个差分方程。这里容易忽略的是截断误差的分析用 Taylor 展开可以验证中心差商的截断误差是 $O(h^2)$而向前差商 $y(x_i) \approx \frac{y_{i1} - y_i}{h}$ 的截断误差只有 $O(h)$。笔记里反复强调“忽略高阶无穷小项用近似值代替精确值”这一步正是有限差分法精度损失的来源。内部截断误差 $R_i$ 是理解差分格式精度等级的钥匙。比如中心差商格式的局部截断误差主项是 $-\frac{h^2}{12} y^{(4)}(\xi)$这意味着步长减半时误差大约缩小到原来的四分之一。实际排查数值结果不精确时先看截断误差阶数再看网格比设置这个顺序不要搞反。2.3 从差分方程到矩阵形式三对角结构的稀疏性引入记号 $u_i \approx y(x_i)$把差分方程整理成如下形式$$ a_i u_{i-1} b_i u_i c_i u_{i1} h^2 f_i, \quad i 1, 2, \dots, N-1 $$其中系数 $a_i, b_i, c_i$ 由 $p(x_i), q(x_i)$ 和步长 $h$ 决定。加上边界条件后全部方程写成一个三对角线性方程组 $A\mathbf{u} \mathbf{g}$。三对角结构的优势在于求解复杂度是 $O(N)$用 Thomas 追赶法就能高效求解不需要动用稠密矩阵的分解算法。import numpy as np def solve_tridiagonal(a, b, c, g): 求解三对角方程组 A u g a: 下对角线长度 N-1 b: 主对角线长度 N c: 上对角线长度 N-1 n len(b) cp np.zeros(n) gp np.zeros(n) cp[0] c[0] / b[0] gp[0] g[0] / b[0] for i in range(1, n): m b[i] - a[i-1] * cp[i-1] cp[i] c[i] / m if i n-1 else 0.0 gp[i] (g[i] - a[i-1] * gp[i-1]) / m u np.zeros(n) u[-1] gp[-1] for i in range(n-2, -1, -1): u[i] gp[i] - cp[i] * u[i1] return u这段代码实现了 Thomas 追赶法也就是三对角方程组的专用消元过程。参数a、b、c分别对应系数矩阵的下对角线、主对角线和上对角线g是右端项。追赶法分两步第一步从前往后消去下对角线元素得到上三角形式的等效方程组第二步从后往前回代求出 $u$。边界条件是通过修改b[0]、b[-1]和g[0]、g[-1]来实现的——这是最容易出错的地方。提示Python 的numpy.linalg.solve也可以解三对角方程组但通用稠密解法的复杂度是 $O(N^3)$节点数超过几千后速度差异会非常明显。3. 基函数法与配置法样条逼近、Runge 现象与最小二乘的取舍3.1 多项式插值的局限Runge 现象为什么让高次插值翻车基函数法的思路是假设解可以表示为一组基函数的线性组合 $\tilde{y}(x) \sum_{k} a_k \varphi_k(x)$然后通过某种准则确定系数 $a_k$。笔记专门提醒了一个坑如果取 $\varphi_k(x) x^k$也就是多项式插值型基函数当次数升高时会出现 Runge 现象——插值多项式在区间端点附近剧烈振荡而且次数越高振荡越严重。这不是数值实现的问题而是多项式插值本身的数学性质决定的。等距节点下的高次多项式插值不收敛于光滑函数拟合精度反而随次数提升而恶化。所以实际计算中笔记建议用三次 B 样条作为基函数。二阶常微分方程的数值解要求逼近函数具有二阶连续导数三次 B 样条恰好满足 $C^2$ 光滑性这是它被选中的关键原因。3.2 三次 B 样条基函数的构成与代入三次 B 样条基函数 $\varphi_k(x)$ 在四个相邻区间上有定义每个基函数在一个局部支集上非零。B 样条基函数的基本性质是关于中心点对称且不同节点的基函数可以通过平移得到。这意味着构造全局逼近时每个节点的值只受相邻三个基函数的影响。把 $\tilde{y}(x) \sum_{k} a_k B_k(x)$ 代入方程后在节点 $x_i$ 处取值可得到如下关系$$ \tilde{y}(x_i) \frac{1}{6}(a_{i-1} 4a_i a_{i1}) $$这个式子描述了样条系数与节点函数值之间的对应关系。代入微分方程后会得到一个关于 $a_i$ 的代数方程组。边界条件的作用是消除首尾两个虚拟未知量 $a_{-1}$ 和 $a_{N1}$——这两个量不在实际区间内部但 B 样条的支集结构使得它们进入了节点 $x_0$ 和 $x_N$ 处的方程必须通过边界条件消去。这正是基函数法实操里最容易被忽略的步骤。3.3 配置法与最小二乘法的对比残差为零点 vs 残差积分最小两种确定系数的方法是配置法与最小二乘法。配置法要求残差在若干个配置点上为零残差方程数量等于未知量数量直接得到一个代数方程组。最小二乘法则是让残差的平方积分在整个区间上最小$$ \min \int_\Omega R^2(x) dx $$笔记给出的算例很典型假设逼近函数 $\tilde{y} a_0 a_1 x a_2 x^2$两个边界条件消去两个未知量剩下一个自由参数。用配置法时在区间中点选一个配置点就够用最小二乘法时则令 $\frac{\partial}{\partial a_i} \int R^2 dx 0$ 得到方程。两种方法都能求但适用场景不同——配置法简单直接、控制点灵活但“只保证配置点上精确”最小二乘法全局误差分布更均匀适合对整体精度有要求的场景。import numpy as np from scipy.integrate import quad # 配置法选区间中点作为配置点残差为零 x_c 0.5 R_c lambda a2: (a2 * x_c**2 (1 - a2) * x_c 1) # 假设边界条件已消元 # 最小二乘法对残差平方在区间上积分并求导 def residual(a2, x): return a2 * x**2 (1 - a2) * x 1 def objective(a2): val, _ quad(lambda x: residual(a2, x)**2, 0, 1) return val # 用 scipy 最小化目标函数 from scipy.optimize import minimize_scalar res minimize_scalar(objective) print(配置法解 a2 , minimize_scalar(lambda a2: abs(residual(a2, x_c))).x) print(最小二乘解 a2 , res.x)这段代码用配置法和最小二乘法分别求解同一个含待定参数 $a_2$ 的逼近问题。residual函数返回给定 $a_2$ 时逼近函数的残差配置法直接找使中点和边界残差为零的 $a_2$最小二乘法通过quad计算残差平方积分后再用minimize_scalar找极小值点。注意配置法本质上是一个求根问题最小二乘法是一个优化问题两者数学结构完全不同不要混用求解器。提示最小二乘法对边界条件的处理没有配置法那么直观。配置法可以直接把边界条件代入消元最小二乘法则要在目标函数里加入边界残差惩罚项或者用约束优化。笔记中的算例因为自由度少所以能做显式消元实际项目里自由度多时建议用带约束的优化器。4. 抛物型方程三种古典差分格式显式、隐式与 Richardson 的精度与稳定4.1 最简显式格式向前差商时间、中心差商空间一维非齐次热传导方程 $\frac{\partial u}{\partial t} a\frac{\partial^2 u}{\partial x^2} f(x,t)$ 的数值解笔记给出了三种古典格式。最简显式格式在时间方向用一阶向前差商空间方向用二阶中心差商$$ \frac{u_j^{k1} - u_j^k}{k} a\frac{u_{j1}^k - 2u_j^k u_{j-1}^k}{h^2} f_j^k $$引入网格比 $r \frac{ak}{h^2}$整理后得到递推公式 $u_j^{k1} (1 - 2r)u_j^k r(u_{j-1}^k u_{j1}^k) k f_j^k$。这个格式的局部截断误差是 $O(k h^2)$时间方向一阶精度、空间方向二阶精度。递推计算非常简单——每一层新值只依赖上一层相邻三个点的值用向量化操作就能一次算出整层。import numpy as np import matplotlib.pyplot as plt # 参数设置 L 1.0 J 101 # 空间节点数 h L / (J - 1) a 1.0 k 0.0004 # 时间步长 r a * k / h**2 print(网格比 r , r) T 0.1 N int(T / k) x np.linspace(0, L, J) u np.sin(np.pi * x) # 初始条件 u[0] 0.0; u[-1] 0.0 # 齐次边界条件 for n in range(N): u_new np.zeros_like(u) u_new[1:-1] (1 - 2*r) * u[1:-1] r * (u[:-2] u[2:]) u_new[0] 0.0 u_new[-1] 0.0 u u_new plt.plot(x, u) plt.title(ft {T}, r {r}) plt.show()这段代码实现了最简显式格式的完整推进过程。核心参数是网格比r它由时间步长k、空间步长h和扩散系数a共同决定。u[:-2]和u[2:]分别取左邻居和右邻居序列实现向量化的空间差分。边界条件在每个时间步后强制写回保证边界上的值是齐次条件。如果r大于 0.5运行结果会出现振荡这就是后面稳定性分析要解决的问题。4.2 最简隐式格式无条件稳定的代价是解方程组最简隐式格式把空间二阶中心差商放到 $k1$ 时间层$$ \frac{u_j^{k1} - u_j^k}{k} a\frac{u_{j1}^{k1} - 2u_j^{k1} u_{j-1}^{k1}}{h^2} f_j^{k1} $$整理后每个时间步都需要求解一个三对角方程组。显式格式一次向量化计算就完成推进隐式格式却要解 $N-1$ 元的线性方程组。用 Thomas 追赶法求解三对角方程组复杂度仍然是 $O(N)$不会成为性能瓶颈。这笔开销换来的收益是稳定性——Fourier 分析得到的传播因子 $G(\theta) \frac{1}{1 4r \sin^2(\theta/2)}$对任意网格比 $r 0$ 都有 $|G(\theta)| \le 1$因此隐式格式是绝对稳定的也就是常说的无条件稳定。实际工程计算里显式格式受限于 $r \le 0.5$ 的条件空间网格加密时时间步长必须跟着缩小计算成本迅速上升。隐式格式可以独立选择时间步长精度允许的情况下用大步长推进整体效率往往更高。4.3 Richardson 格式中心差商带来的三层结构Richardson 格式在时间方向用中心差商局部截断误差提升到 $O(k^2 h^2)$$$ \frac{u_j^{k1} - u_j^{k-1}}{2k} a\frac{u_{j1}^k - 2u_j^k u_{j-1}^k}{h^2} f_j^k $$这个格式形式上比隐式简单精度也更高但它是一个三层格式——计算 $k1$ 层需要知道 $k$ 层和 $k-1$ 层的数据第一层需要特殊处理。更麻烦的是它的传播因子 $G(\theta) -4ar\sin^2(\theta/2) \pm \sqrt{1 16a^2 r^2 \sin^4(\theta/2)}$当 $\theta \to \pi$ 时 $|G(\theta)| 1$对任意步长都不满足 Von Neumann 条件因此是绝对不稳定的。笔记里用了转化技巧补充恒等式 $u^{k} u^{k}$把三层格式改写为等价的二层方程组然后用传播矩阵的特征值分析判断稳定性。这里的教训是三层格式不能直接套用二层格式的稳定性结论必须先降维。实际工程中 Richardson 格式很少直接用但它的分析过程是理解 Fourier 稳定性方法的好素材。5. 稳定性与收敛性避坑Von Neumann 条件、网格比和三个踩坑记录5.1 相容性、稳定性与收敛性Lax 等价定理划定的边界差分格式的三个基本概念是相容性、稳定性和收敛性。相容性说的是差分方程在步长趋于零时逼近微分方程也就是局部截断误差 $|R_i| \to 0$。稳定性说的是计算过程中误差不会随时间层放大约束即 $|U^k| \le M|U^0|$。收敛性说的是差分格式的数值解在步长趋零时逼近真解。三者之间的关系由 Lax 等价定理给出对一个适定的线性初值问题差分格式收敛的充要条件是相容性加稳定性。这意味着稳定性和收敛性不能分开讨论——格式是相容但不稳定数值解照样不收敛格式稳定但不相容收敛到的是错误方程的解。期末复习时这三者的定义必须能一字不差地默写简答题和判断题都爱考。5.2 传播因子推导用 Fourier 分析判别稳定性的套路Fourier 稳定性分析的标准流程是将 $u_j^k V^k e^{i\theta j}$ 代入齐次差分格式消去公因子后得到 $V^{k1} G(\theta) V^k$其中 $G(\theta)$ 就是传播因子。然后检查 $|G(\theta)| \le 1$ 是否对所有 $\theta$ 成立。最简显式格式代人后得到 $G(\theta) 1 - 2r 2r\cos\theta 1 - 4r\sin^2(\theta/2)$要求 $|G(\theta)| \le 1$ 等价于 $r \le \frac{1}{2}$。这个推导过程建议完整手推一遍不要只记结论——考试的时候给一个具体格式判断稳定性的第一步就是写传播因子第二步才是判断模长。5.3 踩坑记录三个实际计算中反复出现的问题坑一显式格式计算出现振荡现象用最简显式格式求解热传导方程某些时间层上数值解出现锯齿状振荡且振荡位置固定在某些空间节点附近。原因网格比 $r \frac{ak}{h^2}$ 超过了 0.5例如 $h0.1, k0.01, a1$ 时 $r1$已经远超稳定限。传播因子在部分频率上大于 1误差逐层放大。解决先打印网格比检查再决定是否调整步长。保持 $h$ 不变时取 $k \le \frac{h^2}{2a}$或者改用隐式格式让时间步长不再受稳定性约束。坑二边值问题的矩阵方程解出来是错的但格式本身没问题现象有限差分格式推导无误代入边界条件后求解结果与解析解偏差很大检查边界附近误差尤其明显。原因边界条件的离散方式有误。最常见的错误是把边界条件 $y(x_0)$ 用一阶向前差商近似精度只有 $O(h)$整体精度被边界降了一阶。正确做法是用中心差商近似导数边界条件或者引入虚拟节点 $x_{-1}$ 并用中心差商消去。解决边界条件单独构造离散方程不要让边界点掉出离散体系。遇到导数边界条件时先补虚拟节点再消元。坑三三层格式直接套用两层稳定判据现象用 Richardson 格式计算时不管怎么缩步长结果总是发散。原因三层格式不能直接套用二层格式的传播因子判据。三层格式需要补充恒等式转化为二层方程组然后分析传播矩阵的特征值。Richardson 格式的传播矩阵特征值模长大于 1任何步长下都不稳定这是格式本身的数学性质和参数选择无关。解决遇到三层格式先做降维处理再分析稳定性。如果只是需要二阶时间精度换成 Crank-Nicolson 格式——它是二层的精度 $O(k^2 h^2)$且无条件稳定。提示判断一个差分格式的稳定性习惯上先用 Von Neumann 条件快速筛一遍再用数值实验验证取几组不同的网格比跑几步看有没有振荡。两个方法都过才算靠谱。6. 期末复习自查把笔记用起来的三个验证步骤复习这份笔记最有效的方式不是重读一遍而是按步骤动手推。我的习惯是考前三天每天做一遍以下三步自查能发现不少以为自己会了其实没掌握的盲区。第一步拿一张白纸从零推导二阶常微分方程的有限差分格式。写出中心差商公式、代入方程、合并系数、整理成三对角矩阵形式最后解释边界条件如何进入矩阵。这一步能检验你对第二章内容的掌握程度推不出来的地方就是需要补的地方。特别留意截断误差的主项考试常考“这个格式的精度是几阶”这类题。第二步用显式格式求解一个最简单的热传导方程 $u_t u_{xx}$取 $h 0.1$分别试 $k 0.004$对应 $r 0.4$和 $k 0.006$对应 $r 0.6$对比数值解在 $t 0.1$ 时刻的表现。这一步的目的不是算题而是直观感受稳定性边界在哪。跑完这两组参数后你再去看 Von Neumann 条件的推导理解完全不一样。第三步用配置法和最小二乘法解同一个二阶边值问题对比两种方法得到的逼近函数曲线。这能帮你理解“残差为零”和“残差积分最小”的差别——前者在控制点上精确但控制点之间有偏差后者整体误差分布更均匀。实际做的时候可以用笔记里的算例把逼近函数换成带一个自由参数的多项式手算几步对比结果。有一年期末复习时我跳过第一步直接刷题结果考试碰到一道“写出差分格式并说明边界条件如何离散”的题目当场卡壳——平时作业都是电脑算的没仔细推过矩阵是怎么构造出来的。从那以后每次复习数值方法相关的课程我都强制自己先过一遍“推格式、算网格比、对比两种逼近方法”这三步再去做题。这套流程对这份笔记同样适用希望帮到你。本文还有配套的精品资源点击获取
