简介《科学与工程计算》熊春光版期末复习笔记适合北理工自动化专业及相关理工科学生期末复习与考前突击也可用于预习时把握重点。压缩包内为一份PDF文档共1个文件约20.04MB便于屏读与打印标注。笔记围绕微分方程数值解法展开系统覆盖初值问题与边界条件、有限差分法中的等距划分、步长与差商选取、内部截断误差估计以及基函数法、Lagrange插值与三次样条插值等核心考点同时结合二阶常微分方程定解、矩阵形式代数方程组求解等案例帮助理解离散化思路。对高次插值的Runge现象与三次样条的连续性差异也做了辨析并梳理了B样条基函数图形对称、平移构造解空间等特点便于考前快速扫除盲区。已有2754人学习/下载是一份高效的期末复习资料。1. 期末复习这门课先把边值问题和抛物方程两条主线钉住《科学与工程计算》这门课对自动化专业的学生来说最容易踩的坑是知识点太散边界条件、差分格式、稳定性分析、插值逼近每块单独都能看懂但合在一起就是串不起来。这版期末复习笔记的好处在于它把第1章两点边值问题的数值解法和第4章抛物方程差分格式整理成了同一条主线先给定解条件再选离散格式最后判断稳定性和收敛性。把这条线理清楚考前一周足够把全课内容收进一张纸。适合平时基础一般、考前需要快速建立体系的人读也适合已经学过一遍但面对题目仍不知道用显式还是隐式格式的人。2. 从边值问题到三种差分格式定解条件、截断误差与矩阵组装2.1 初值问题和边值问题定解条件到底差在哪笔记一开篇就把定解条件分成两类与时间有关的叫初始条件与空间边界有关的叫边界条件。这个区分看起来简单但很多翻车现场都出在这里。初值问题给的是某个时间点的状态比如控制系统的状态方程给出 t0 时刻 x(0)这是初始条件而边值问题给的是空间区域两端的行为比如一根杆两端分别固定为某个温度这是边界条件。二阶常微分方程要求定解区间两个端点都有条件缺一个方程组就欠定多一个又矛盾。我一般会建议拿到题目先在草稿纸上画坐标轴把题目给的条件按时间/空间和已知值/已知导数标出来。y(0)0 和 y(0)0 是初始条件y(0)0 和 y(1)0 才是两点边值条件。把这两类混在一起后面所有格式推导都会跟着错。这版笔记把这一点放在第1章最前面复习时值得多停留两分钟。2.2 有限差分法的第一步等距剖分与中心差商有限差分法的思路是不要直接解微分方程而是把自变量区间等距切开在离散节点上近似导数。把 [a,b] 分成 N 等份节点记为 x_jajhh(b-a)/N 就是步长。对导数用有限差商逼近时一阶导通常取中心差商u(x_i) ≈ (u_{i1} − u_{i−1}) / (2h)二阶导取二阶中心差商u(x_i) ≈ (u_{i−1} − 2u_i u_{i1}) / h²为什么取中心差商而不是前向差商把 Taylor 展开写一遍就知道前向差商的截断误差是 O(h)中心差商是 O(h²)。对二阶常微分方程 y p(x)y q(x)y f(x) 来说把所有节点上的差商代进去就能整理成一个三对角线性方程组。用矩阵形式写就是笔记里的 Ag b其中 A 是三对角矩阵对角线由 q(x_i) 和 2/h² 组成次对角线来自 −1/h² 和 p(x_i)/(2h)。实际求解三对角方程组最常用的手段是追赶法我习惯把它封装成通用函数因为后面隐式格式每推进一步都要解一次三对角方程组import numpy as np def solve_tridiag(sub, diag, sup, rhs): 追赶法解三对角方程组sub/diag/sup/rhs 长度按矩阵结构对齐 n len(diag) cp np.zeros(n) # 消元后的上三角系数 dp np.zeros(n) # 消元后的右端项 cp[0] sup[0] / diag[0] dp[0] rhs[0] / diag[0] for i in range(1, n): m diag[i] - sub[i - 1] * cp[i - 1] cp[i] sup[i] / m if i n - 1 else 0.0 dp[i] (rhs[i] - sub[i - 1] * dp[i - 1]) / m x np.zeros(n) x[-1] dp[-1] for i in range(n - 2, -1, -1): x[i] dp[i] - cp[i] * x[i 1] return x # 示例-u u 1, u(0)u(1)0取 h0.1内部节点数 n9 n 9 h 0.1 diag (2.0 / h**2 1.0) * np.ones(n) # 对角线 sub sup -1.0 / h**2 * np.ones(n - 1) # 上下次对角线 y solve_tridiag(sub, diag, sup, np.ones(n)) print(y)这里 diag 是内部节点的二阶差分系数加源项系数sub 和 sup 是相邻节点耦合项。n 取 9 时步长 0.1矩阵规模很小主要用于验证追赶法实现是否正确。实际算题时N 越大截断误差 O(h²) 越小但矩阵条件数会变大不是步长越小越好。2.3 最简显式格式递推只要一行公式到了抛物方程部分笔记先给出热传导方程的定解问题u_t a u_xx f(x,t)加上初始条件 u(x,0)φ(x) 和两个边界条件。最自然的离散方式是时间方向用一阶向前差商空间方向用二阶中心差商这就是最简显式格式也叫向前差分格式(u_j^{k1} − u_j^k) / τ a(u_{j−1}^k − 2u_j^k u_{j1}^k) / h² f_j^k这个式子里的上标 k 表示第 k 个时间层下标 j 表示第 j 个空间节点。τ 是时间步长h 是空间步长。把 u_j^{k1} 单独提出来引入网格比 r aτ / h²就得到u_j^{k1} r u_{j−1}^k (1 − 2r)u_j^k r u_{j1}^k τ f_j^k这个格式最舒服的地方是递推公式直接可用当前时刻的空间值全部是已知的逐点算一遍就得到下一个时间层的值不需要解任何方程组。代价是稳定性受网格比约束具体约束是 r ≤ 1/2这个结论下一章会从传播因子的角度重新推一遍。复习时应当记住显式格式的局部截断误差是 O(τ h²)时间方向只有一阶精度。2.4 最简隐式格式同样的精度无条件稳定如果把时间方向改成向后差商就得到最简隐式格式也叫向后差分格式(u_j^{k1} − u_j^k) / τ a(u_{j−1}^{k1} − 2u_j^{k1} u_{j1}^{k1}) / h² f_j^k注意这里的空间二阶差商全部取在 k1 时间层。也就是说u_j^{k1} 不再能单独提出来它和左右两个相邻节点的未知值耦合在一起必须每推进一层就解一次三对角线性方程组。这正是上一小节那个追赶法函数的用武之地。隐式格式的局部截断误差同样是 O(τ h²)和显式格式同级但稳定性却是无条件稳定时间和空间步长怎么取都不会让误差放大。工程上如果计算域大、时间要推到很远隐式格式是更稳妥的选择。笔记里特别提了一句隐式格式形成的三对角矩阵是严格对角占优的矩阵本身性质好追赶法不需要选主元。2.5 Richardson 格式时间二阶精度为什么是陷阱前面的显式和隐式格式时间方向都只有一阶精度。想提高时间精度直觉是时间方向也换成中心差商这就是 Richardson 格式(u_j^{k1} − u_j^{k−1}) / (2τ) a(u_{j−1}^k − 2u_j^k u_{j1}^k) / h²时间方向用了中心差商之后截断误差确实变成 O(τ² h²)从精度看比前两个格式都漂亮。但它有两个问题第一它是三层格式计算 u^{k1} 需要用到 u^k 和 u^{k−1} 两层的信息启动时需要额外处理第一层第二也是最致命的这个格式无条件不稳定无论 τ 和 h 取多少误差都会随时间层增长。笔记在推导传播因子时会验证这一点复习阶段先记住结论想提高时间精度不能简单地把时间差商换成中心差商。3. 稳定性判断的手算过程传播因子、Von Neumann 条件和 r≤1/2 是怎么来的3.1 传播因子是什么把误差拆成谐波来看差分格式逐层向前推进某一层的计算误差会传到下一层。稳定性讨论的本质是问这个误差是随着层数增长还是被压制对于线性常系数差分格式Fourier 分析是标准工具。基本做法是把误差函数分解成谐波分量 e^{ijωh}看单个谐波在时间推进中放大多少倍。把近似解写成 u_j^k G^k e^{ijωh}代进齐次差分格式后指数因子可以约掉剩下的 G(ω) 就叫传播因子或者增长因子。为什么只讨论齐次格式笔记里有一段定理说明如果差分格式按范数关于初值稳定那么它也关于右端项稳定所以稳定性讨论可以归结为齐次方程组的初值稳定性。这个结论很重要它让我们不用每次都拖着 f 项做分析。3.2 显式格式稳定域的手算过程由 |G|≤1 推出 r≤1/2把 u_j^k G^k e^{ijωh} 代入最简显式格式的齐次形式。等号左边是 G^{k1}e^{ijωh}右边三项都含 G^k e^{ijωh}约掉公共因子后得到G(ω) r e^{−iωh} (1 − 2r) r e^{iωh}用欧拉公式把两个复指数合并e^{iωh} e^{−iωh} 2cos(ωh)于是G(ω) 1 − 2r 2r cos(ωh) 1 − 4r sin²(ωh/2)到这里问题就变成ω 取任意实数时G(ω) 的模是否不超过 1。sin²(ωh/2) 的取值范围是 0 到 1所以 G(ω) 落在 1−4r 到 1 之间。要让所有可能取值都不超过 1等价于要求下界不小于 −11 − 4r ≥ −1即 r ≤ 1/2注意这里不能只检查 |1−4r| ≤ 1因为 G 的上界是 1下界才是关键约束。实际算题时 r 一旦超过 0.5数值解就会开始振荡e^{ijωh} 中某些频率的谐波会被持续放大表现就是时间层越往后波形越乱。3.3 隐式格式无条件稳定的根源分母永远比 1 大最简隐式格式做同样的代入区别在于空间差商全部取在 k1 层。整理后传播因子变成G(ω) 1 / [1 4r sin²(ωh/2)]这个式子分母最小是 1最大是 14r不管 r 取多少|G(ω)| 恒小于等于 1。这就是无条件稳定的数学来源隐式格式是求解而不是递推每个新时间层的值都通过线性方程组与周围节点耦合单谐波误差天然被分母压缩。直觉上隐式格式更像在解一个稳态问题而不是把旧信息直接外推。3.4 Richardson 格式如何判断三层格式要补恒等式化成两层Richardson 格式是三层格式不能直接套两层格式的传播因子公式。笔记里的处理技巧是补一个恒等式 u_j^k u_j^k把三层格式改写成一个向量形式的二层格式u_j^{k1} u_j^{k−1} − 8r sin²(ωh/2) u_j^k配合恒等式 u_j^k u_j^k写成矩阵形式后传播矩阵的特征方程是λ² 8r sin²(ωh/2)λ − 1 0两个特征根的乘积等于常数项 −1也就是说不管 r 取多少至少有一个特征根的模不小于 1。这直接违反 Von Neumann 条件因此 Richardson 格式无条件不稳定。笔记里还特地提到如果传播矩阵的特征值有重根需要检查 Jordan 块的情况对 Richardson 格式来说特征根乘积为 −1 已经足够说明问题。3.5 相容性、稳定性与收敛性等价定理到底说了什么三个概念之间的关系很容易背混。相容性说的是差分格式的局部截断误差在 τ 和 h 趋于 0 时趋于 0即差分方程确实是微分方程的近似稳定性说的是误差在逐层推进中不失控收敛性说的是数值解趋近真解。笔记里的 Lax 等价定理把三者串在一起对一个适定的线性初值问题相容性是基础在相容的前提下稳定性是收敛性的充分必要条件。这个定理的实际指导意义是检查一个新格式先看截断误差是否是 O(τ^p h^q)再算传播因子判断稳定两者都满足收敛就不用额外验证。工程实践里很多人只做误差分析不做稳定性分析这是典型的理论挂空挡。我自己写差分程序时的习惯是先随手打印第一层的 r 值和传播因子模的最大值确认稳定域满足再做长时间推进。下面用一个小实验验证显式格式的稳定域直接对比 r0.4 和 r0.6 的表现import numpy as np L, T 1.0, 0.06 N 50 # 空间网格数 h L / N # 空间步长 a 1.0 # 扩散系数对应 u_t a u_xx x np.linspace(0, L, N 1) u0 np.sin(np.pi * x) # 满足零边界条件的初始波形 for r in [0.4, 0.6]: tau r * h * h / a M int(round(T / tau)) u np.zeros((M 1, N 1)) u[0, :] u0 for k in range(M): u[k 1, 1:-1] ( r * u[k, :-2] (1 - 2 * r) * u[k, 1:-1] r * u[k, 2:] ) exact np.exp(-np.pi**2 * T) * np.sin(np.pi * x) err np.max(np.abs(u[M] - exact)) print(fr{r}, tau{tau:.6f}, M{M}, max_err{err:.4e})这段实验限制了只做一步到的最大偏差比较。r 从 tau 计算得到tau r·h²/aM 取 T/tau 的整数值T 选 0.06 是为了让两种 r 下时间层数都能整除。解析解是 e^{−π²t}sin(πx)初始波形取热传导方程的第一特征函数空间离散本身误差很小最终误差主要反映时间推进的稳定性。r0.4 时误差保持在很小量级r0.6 时误差会指数级放大这个数字差异比任何文字描述都有说服力。4. 基函数法、配置法与最小二乘法多项式不行就换样条4.1 基函数法的总体思路把连续问题转换成有限个系数差分法是直接离散导数基函数法则换了一条路先假定近似解的形式把它表示成一组基函数的线性组合再通过某种准则确定组合系数。设近似解为u_N(x) Σ a_k φ_k(x)φ_k(x) 是基函数a_k 是待定系数。所有后续工作——无论是配置法、最小二乘法还是 Galerkin 法——本质上都是在问同一个问题怎样确定 a_k让 u_N(x) 在某种意义下最接近真解。二者的区别只在于接近的标准。笔记里专门讨论了基函数的选取对结果的影响。取不同基函数得到不同离散格式这一句话是整个第4章的枢纽。基函数选得好方程组的矩阵紧凑、条件数友好选得不好即使理论上能逼近实际算出来也是一团乱。4.2 Lagrange 插值为什么容易翻车Runge 现象最直接的想法是取多项式基函数比如 φ_k(x) x^k或者更讲究一点取 Lagrange 插值基函数让 φ_k 在节点 x_k 上取值 1在其他节点上取值 0。Lagrange 基函数构造简单插值多项式唯一看起来非常理想。翻车点在于高次多项式。节点一多插值多项式在两个端部附近会出现剧烈振荡这就是 Runge 现象。对一个光滑函数做等距节点插值次数越高端部振荡越严重近似效果反而更差。笔记里这个判断写得很直白高次多项式插值会出现 Runge 现象。所以在实际数值方法中几乎不会用高次多项式对全域做插值要么分区用低次多项式要么换样条基函数。4.3 三次 B 样条二阶连续可导、紧支撑的折中二阶常微分方程的解至少要求近似函数二阶连续可导否则差商和残差都没有意义。三次 B 样条正好满足 C² 连续而且每个基函数只在四个相邻网格区间上有非零值这个性质叫紧支撑它保证离散后的矩阵是稀疏的不会产生全局耦合。以标准三次 B 样条为例分段定义在区间 [0,4] 上步长取 1 时的表达式是0 ≤ x 1B(x) x³ / 61 ≤ x 2B(x) (−3x³ 12x² − 12x 4) / 62 ≤ x 3B(x) (3x³ − 24x² 60x − 44) / 63 ≤ x ≤ 4B(x) (−x³ 12x² − 48x 64) / 6每个基函数的图形关于中心对称最大值为 1不同节点上的基函数可以由这个基本图形平移得到。我一般会提醒自己正因为每个基函数只覆盖最多三个区间配置法最后得到的矩阵是窄带宽的而不是满矩阵。这一点对数值求解效率是决定性的。4.4 配置法与最小二乘法残差为零还是残差最小配置法的思路是把近似解代入微分方程后一般不会精确满足会剩下一个残差 R(x; a)。选取配置点 x_i强制残差在这些点上等于 0得到代数方程组 R(x_i) 0。基函数有 N 个系数边界条件通常给出 2 个约束所以配置点数量取 N−2 个即可。注意配置点个数不能乱取取多了方程组超定取少了欠定。配置法只保证有限个点上残差为零残差函数整体可能仍然很大。笔记因此引入最小二乘法目标改为让残差平方的积分极小min ∫ R(x; a)² dx对 a_k 求偏导并令其为零得到一组线性方程。和配置法相比最小二乘法需要考虑积分计算计算开销更大但整体残差被强制压到最小。实际使用中如果只是想要一个快速近似解我一般会选配置法并精心挑选配置点比如取 Gauss 点如果想要一个全局稳定的近似就走最小二乘或 Galerkin 路线。Galerkin 法取权函数与基函数相同本质上是最小二乘的变体笔记里把它归入同一个思想脉络。5. 期末复习避坑日记五个最容易翻车的细节5.1 网格比漏乘扩散系数 a现象显式格式计算热传导方程程序写出来和笔记一模一样但结果发散检查后发现 r 用的是 τ/h²不是 aτ/h²。原因最简显式格式的递推公式里系数 r aτ/h²a 是方程 u_t a u_xx 里的扩散系数。题目给 a2 时按 τ/h² 算出的网格比直接小了一倍稳定域判断全错。解决写代码前先单独打印一行 r 的值对照题目给的 a、τ、h 手动算一遍确认。稳定域边界是 r ≤ 1/2如果 a 不是 1必须体现在 r 中。5.2 初始条件和边界条件混为一谈现象遇到 y p y q y f 的题目题目给 y(0)0 和 y(0)0有人直接按两点边值问题处理硬套差分格式。原因y(0)0 和 y(0)0 都是 t0 时刻的条件属于初值问题。边值问题必须在空间区间两端分别给条件比如 y(0)0 和 y(1)0。题目给的是初值差分格式却按边值问题列三对角方程组矩阵维数都对不上。解决拿到题目先花十秒画坐标轴在图上标出每个条件所在的位置。初值问题用时间推进格式边值问题用空间离散格式这是两条完全不同的路线。5.3 Richardson 格式翻车后还在加密步长现象用 Richardson 格式算抛物方程误差随网格加密反而增大数值解出现锯齿状振荡。最开始以为步长不够小把 h 和 τ 同时减半结果更糟。原因Richardson 格式是三层格式传播矩阵的特征根乘积为 −1无条件不满足 Von Neumann 条件。无论步长多小它都不稳定这不是精度问题是格式本身的问题。解决考前记住一个结论想提高时间方向精度不能用 Richardson 格式硬冲正路是 Crank-Nicolson 这类在空间上也做中间层处理的格式。看到三层格式第一反应是检查稳定性而不是检查步长。5.4 配置点个数数错导致欠定方程组现象用三次 B 样条做配置法取了 N 个基函数边界条件 2 个配置点随手取了 N/2 个最后方程个数比未知量少线性方程组解不出来。原因配置点的数量不是随意取的。未知量总数减去边界条件约束数才是可用的配置点数。B 样条基函数数量多配置点也要对应增加多项式基函数少配置点就少。解决动手之前先数数基函数个数 m边界条件个数 b配置点个数取 m−b。笔记里那个多项式近似例子里3 个系数配 2 个边界条件只需要 1 个配置点就是这个逻辑。5.5 只会背 O(h²)不会写截断误差主项现象考试要求写出某格式的局部截断误差表达式有人只写了一句误差是 O(h²)被扣分。原因O(h²) 是阶数不是表达式。截断误差计算要把真解 u(x_j, t_k) 代入差分格式做 Taylor 展开后收集最低阶非零项主项里必须带系数和偏导项。笔记里显式格式的截断误差就写成了 R (τ/2)u_tt − (a h²/12)u_xxxx … 这样的形式。解决复习时每个格式至少手推一遍 Taylor 展开。顺序固定把 u 在展开点泰勒展开代入差分方程逐项合并找出第一个不消失的项。这个过程练三遍之后考场上不会再慌。6. 考前两小时用一张表和两道题把三种格式钉死6.1 考前对照表三种格式一句话记住三种古典格式的差异集中在这张表上考前扫一眼就能回忆全部要点格式时间精度空间精度稳定性每层求解方式最简显式O(τ)O(h²)r ≤ 1/2 才稳定直接递推无需求解方程组最简隐式O(τ)O(h²)无条件稳定每层解三对角方程组RichardsonO(τ²)O(h²)无条件不稳定三层递推无法实用显式格式省在每步计算量小但时间步长被网格比卡死隐式格式每步要解三对角方程组但时间步长可以放到很大Richardson 的存在价值就是给一张反面典型提醒我们时间精度不是靠换差商提上去的。表格里每层求解方式这一列决定了对程序实现的选择显式循环体一行公式隐式循环体里要调用三对角求解器。6.2 两道自测题第一题是稳定域计算给定 u_t 2u_xxh 0.05τ 0.001。先算网格比 r 2 × 0.001 / 0.0025 0.8大于 0.5用显式格式会不稳定。要让显式格式稳定τ 必须不大于 h²/(2a) 0.000625。如果题目要求 τ 不变就只能换隐式格式。这道题每年都会有变体核心是 r 的定义和 1/2 这个边界值。第二题是配置法的最简单演练求解 y x²边界条件 y(0)0、y(1)0设近似解 u_N(x) a·x(1−x)。这个试函数本身满足两个边界条件只有一个待定系数 a。因为只用一个配置点自然取区间中点 x0.5。残差为 R(x) u_N − x² −2a − x²令 R(0.5) −2a − 0.25 0解得 a −0.125。这道题的价值在于展示整个流程写近似解、算残差、选配置点、解系数四步走完。配置点数量不足或残差写错在这个例子里会马上暴露。从那以后我每次考前都会把三种格式的传播因子手推一遍再默写一遍边界条件分类最后用这张表复盘一遍。稳不稳定、几阶精度、要不要解方程组三句话全记住。希望对正在复习的你也有帮助。本文还有配套的精品资源点击获取