简介《常微分方程边值问题的数值解法.pdf》面向学习数值分析、计算数学及工程建模的读者系统梳理二阶常微分方程边值问题的数值求解思路。内容从边值问题与初值问题的差异入手围绕第一、二、三类边界条件展开重点讲解差分法、有限元法以及将边值问题转化为初值问题三类基本途径并涉及中心差分格式、网格剖分、截断误差、三对角方程组追赶法求解、极值原理、解的存在唯一性、收敛性与误差估计等知识点配有算例和结果对照表。压缩包仅含1个PDF文件大小约9.27MB便于集中阅读和按章节查阅。已有1480人学习适合需要理解边值问题离散化过程、掌握差分格式构造与精度分析并在课程学习或工程计算中查漏补缺的读者使用。1. 两端都钉住的曲线常微分方程边值问题到底在算什么把一根钢丝的两端钉在墙上中间挂一个小球它会垂成一条弧线。工程师想知道这条弧线上每一点的挠度于是写出描述它的常微分方程。方程只规定曲线该怎样弯曲却不告诉你它从哪里出发、到哪里结束两端的位置必须另外给定。只含边界条件作为定解条件的那类求解问题就叫边值问题而把纸面上的方程算成一张数值表的过程就是边值问题的数值解法。它值得单独花时间搞懂原因很实在现实里的边界条件几乎都凑不出解析解能算出来的那部分全靠一套成熟的数值算法顶着。2. 为什么边值问题比初值问题难缠定解条件的位置决定了一切同样的二阶微分方程把条件放在同一个点和放在两个端点表面上只差一个位置实际却是两种完全不同的数学对象。位置一变能用的求解思路、解是否存在、数值方法会不会翻车全都跟着变。搞清楚这一点后面三种方法各自的取舍才有落脚点。2.1 初值问题像推多米诺边值问题像拉一根两端固定的绳初值问题给的是某个起点上的函数值和导数值比如 y(a)y₀、y′(a)y₀′。有了这两条信息方程就变成了一台往前推的机器从 a 点迈出一小步算出下一点的值再迈一步一路推到 b 点。这个画面像推倒一排多米诺骨牌第一块往哪边倒后面的走向就全定了方向是单向的。边值问题不是这样。它给的是一条边界上一个条件、另一条边界上另一个条件比如 y(a)α、y(b)β。两端都被钉住你没有起步的斜率可用单向往前推的机器直接哑火。更麻烦的是这种钉法不保证一定有解。举一个经典的对照方程 y″π²y0如果边界是 y(0)0、y(1)0那么 yC·sin(πx) 对任意常数 C 都成立解有无穷多个把边界换成 y(0)0、y(1)1则一个解都没有。这两个例子说明的是同一件事边值问题的难不是数值方法没写好造成的而是问题本身就可能没有唯一答案。所以教材在讲任何算法之前必须先花一节讨论解的存在唯一性——没有这个前提算出来的曲线到底是不是那个解都说不清。2.2 三类边界条件给值、给斜率、还是给两者的组合边界条件不是只有端点值已知一种写法。教材按边界上约束的形式把它分成三类这三类的区别直接决定了数值方法里哪些方程要改写、哪些未知量要保留。第一类边界条件给定函数值本身y(a)A、y(b)B。物理上相当于端点位置被死死固定住。这类条件最省事因为端点的值直接就是已知量代入方程组即可。第二类边界条件给定导数值y′(a)C、y′(b)D。物理上相当于端点的倾斜角被规定但高度可以自由浮动。一根杆的稳态导热如果某一端包了保温层、热流为零写出来就是温度梯度为零属于这一类。数值上的麻烦在于端点值本身成了未知量得靠导数条件额外补一个方程。第三类边界条件给定函数值与导数的线性组合比如 αyβy′γ。它对应的是端点既不固定也不完全自由而是被弹性支撑拉住、或者与外界流体发生对流换热的情形。换热那个例子里边界上的热流正比于表面温度与外界温度之差展开后正好是这种混合形式。三类条件在离散化后都会变成线性方程但改写的部位不同第一类改写右端项第二类改写系数矩阵的首末行第三类两者都动。写程序时最容易出错的地方往往不是中间的差分格式而是端点那两行的符号。2.3 解一定存在吗那个容易被跳过的唯一性前提教材给出的存在唯一性定理条件大致可以这样理解方程右端函数及其偏导数在区域内连续并且它对未知函数 y 的变化率有界不能出现y 稍微一动右端就剧烈跳变的情况。满足这些条件边值问题就有唯一解。对线性问题还有一个更好用的推论只要系数函数连续且零阶项系数满足一定的符号条件唯一性就能保证。这里要强调充分条件三个字。定理成立解唯一定理不成立解未必不存在只是没法保证。实际计算中如果发现数值解随步长剧烈震荡、怎么缩小 h 都不收敛第一件该做的事不是换算法而是回头检查原问题是不是踩在了唯一性条件的边界上。唯一性还有一层实际意义它是误差估计的分母。只有当真正的解确实存在且唯一时数值解离真解有多远这句话才有意义。后面差分法里用极值原理证明收敛本质上就是在用唯一性做支点。3. 差分法把导数换成差商让微分方程变成三对角方程组在所有边值问题的数值方法里差分法的思路最直白导数不好处理那就用差商替掉它。替完之后微分方程变成一组代数方程剩下的工作交给线性代数。直白不等于简单真正要交代清楚的是三件事——差商为什么这么取、方程组为什么恰好是三对角的、以及凭什么相信这个解会收敛到真解。3.1 网格剖分与中心差商二阶导数最朴素的一次替换第一步是把区间切成小段。把 [a,b] 等分成 n 份得到一串节点 x₀a, x₁, …, xₙb相邻节点的距离 h(b−a)/n 叫步长。第 i 个内部节点上的二阶导数用左右两个邻居来近似y″(xᵢ) ≈ (y_{i1} − 2yᵢ y_{i−1}) / h²这个式子叫二阶中心差商。它为什么长这样用泰勒展开一看就明白把 y(xh) 和 y(x−h) 分别展开再相加一阶导项因为符号相反直接抵消三阶导项同样抵消剩下的主项就是 h²y″除以 h² 就还原出二阶导数。被扔掉的是 h² 量级的项所以这个近似的截断误差是 O(h²)——步长缩小一半误差大约降到四分之一。把近似代回方程 y″f(x)两边同乘 −h²整理成−y_{i−1} 2yᵢ − y_{i1} −h²f(xᵢ)每个内部节点写一个这样的方程一共 n−1 个。两端的 y₀α、yₙβ 是已知的挪到等式右边未知量正好也是 n−1 个方程数和未知数对上可以解了。3.2 差分格式为什么是三对角的追赶法又省在哪里观察上面那个方程它只牵扯 y_{i−1}、yᵢ、y_{i1} 三个量。这意味着系数矩阵的每一行只有主对角线和左右两条副对角线非零其余位置全是零。这种矩阵叫三对角矩阵。矩阵结构不是数学上的巧合而是局部近似的直接后果每个方程只描述一个节点和它左右邻居的关系远处的节点压根不出现。于是存储上只需要三条对角线计算上也就有了专门的解法。普通高斯消元解一个 N 阶稠密方程组运算量是 O(N³)。三对角矩阵用追赶法先从上往下消元、再从下往上回代每一步只处理三个数总运算量降到 O(N)。当网格加密到几千、几万个节点时这个差距不是几倍而是天壤之别。下面是一个可直接运行的 Python 实现。import numpy as np def difference_method(f, a, b, n, alpha, beta): 用中心差分法解 y f(x)边界条件 y(a)alpha, y(b)beta f : 右端函数接收 numpy 数组返回同形状数组 a, b : 求解区间 n : 区间等分数步长 h (b-a)/n alpha : 左端点的函数值 beta : 右端点的函数值 h (b - a) / n x np.linspace(a, b, n 1) # 内部节点共 n-1 个未知量矩阵按 -y_{i-1} 2y_i - y_{i1} 组装 A np.zeros((n - 1, n - 1)) for i in range(n - 1): A[i, i] 2.0 if i 0: A[i, i - 1] -1.0 if i n - 2: A[i, i 1] -1.0 rhs -h * h * f(x[1:-1]) # 右端 -h^2 f_i rhs[0] alpha # 左边界已知值移到右端 rhs[-1] beta # 右边界已知值移到右端 y_inner np.linalg.solve(A, rhs) return x, np.concatenate(([alpha], y_inner, [beta]))代码里最关键的是矩阵 A 的符号。方程是 −y_{i−1}2yᵢ−y_{i1}−h²fᵢ所以主对角线是 2两条副对角线是 −1右端项带负号。边界项之所以是加而不是减是因为 y₀ 和 yₙ 本来在等式左边乘上系数 −1 之后挪到右边就变了正号。符号一错整条曲线会朝反方向弯但程序不会报错只是结果悄悄错掉属于典型的静默失败。实际工程里不会真的用 np.linalg.solve 解三对角矩阵而是换成专门的追赶法或带状矩阵求解器。上面用稠密解法是为了让结构看得清楚。n 取到一万时稠密解法会明显变慢这时候换 scipy 的 solve_banded 是更务实的做法。3.3 极值原理给差分方程组的存在唯一和收敛找一个凭据方程列出来了代码也跑通了但还有两个问题没回答这个差分方程组一定有唯一解吗网格不断加密时差分解真的会收敛到微分方程的解吗教材用极值原理回答这两个问题。它的内容相当朴素假设有一组不全相等的数如果每个内部点都满足 −y_{i−1}2yᵢ−y_{i1} ≥ 0那么这组数里非负的最大值只能出现在两端。用大白话讲如果一串数在每个内部点都向上弯那它最高的地方不会藏在中间只可能在端点。证明用的是反证法假设最大值出现在某个内部点 xₖ那么它的左右邻居都不大于它于是 −y_{k−1}2yₖ−y_{k1} 这个组合一定是小于等于零的而且因为邻居里至少有一个严格更小它其实严格小于零这就和假设的 ≥0 直接冲突。有了极值原理唯一性就顺手了考虑齐次方程组右端全是零、边界也全是零它的解既要满足非负最大值只能是端点值 0又要满足非正最小值只能是端点值 0两头一夹整个解只能是零。齐次方程组只有零解说明系数矩阵可逆非齐次方程组自然有唯一解。收敛性的证法思路是比较函数。构造一个辅助的二次函数让它对应到离散网格上的值恰好能控制住误差再对误差项套用极值原理最后得到形如 |yᵢ−y(xᵢ)| ≤ M·h² 的估计其中 M 和真解的四阶导数有关。这个不等式说明h 趋于零时误差以 h² 的速度趋于零差分法确实收敛。它也解释了为什么网格加密是提升精度的正道——误差是步长的平方级而不是碰运气。3.4 一般二阶方程一阶导数要用三点公式端点还要特殊处理上面处理的是 y″f(x) 这种右端不含 y 和 y′ 的简单形式。一般的二阶线性方程长这样y″ p(x)y′ q(x)y f(x)。多出来的一阶导数项怎么离散是有讲究的。如果还用中心差商 (y_{i1}−y_{i−1})/(2h) 近似 y′截断误差是 O(h²)和二阶导数项的精度匹配整体精度保得住这是正确的选择。如果图省事改用单侧差商 (y_{i1}−yᵢ)/h误差只有 O(h)整个格式立刻降成一阶精度前面的功夫白费。精度匹配这件事是差分格式设计里的一条硬规矩。麻烦出在端点。内部的节点左右都有邻居中心差商直接可用端点 x₀ 左边没有节点只能构造单侧的三点公式比如用 x₀、x₁、x₂ 三点表示 x₀ 处的导数系数是 (−3y₀4y₁−y₂)/(2h)。这个公式的截断误差同样是 O(h²)代价是端点上要引入三个点而不是两个装配矩阵时首末行会多出一个非零元素。把每一项都代入、同类项合并之后得到的依然是一个三对角方程组只是系数比简单情形复杂一些主对角线和副对角线上的元素都变成了 x 的函数。解这个方程组仍然用追赶法。整个过程最考验耐心的地方不在推导而在检查——每一步移项、每一个符号、每一个端点公式都要对一遍。手推三页草稿错一个负号程序跑出来光滑漂亮但完全错误的结果这种坑在数值计算里太常见了。4. 有限元与打靶同一个边值问题另外两条完全不同的路差分法是从微分方程出发直接对着导数下手。另外两条路走的完全是别的方向一条把方程改写成求能量最小的变分问题再从变分形式出发做离散另一条干脆绕开边值问题的框架把它变回初值问题去解。两条路的出发点不同适用的场景也明显不同。4.1 变分形式把解方程换成找一个让能量最小的函数先看一个物理事实一根两端固定的弹性绳在自重下会自然形成一个形状而这个形状恰好使系统的总势能取最小值。有限元法就是从这个角度切入的。对一类线性两点边值问题可以构造一个泛函——注意泛函这个词它不是普通函数而是输入一个函数、输出一个数的映射。这里输入的是一条曲线 y(x)输出的是一个积分值J[y] ∫ₐᵇ (p·y′² q·y² − 2f·y) dx然后提出一个变分问题在所有满足边界条件的函数里找出那个让 J 取最小值的 y。教材的等价性定理说得很明确y 是原边值问题的解当且仅当它使这个泛函取极小。两条路通向同一个答案。必要性方向的证明比较直观如果 y 已经是原方程的解把它代入泛函并和任意一个扰动后的函数比较用分部积分把导数项挪走剩下的项正好是原方程的残差乘以扰动量残差为零泛函的增量就非负说明 y 确实是最小点。充分性反过来从泛函取极小出发对扰动参数求导并令其为零分部积分之后得到的就是原微分方程的弱形式再配合边界条件就还原成原问题。走这条路的好处不只是在纸面上好看。微分方程要求解有二阶导数条件苛刻变分形式里只出现一阶导数对解的光滑性要求降低了一个档次。这意味着一些在经典意义下不可导、但物理上合理的解也能被纳入进来。更实际的好处是变分形式里的积分可以拆成一个个小区间上的积分再相加这个拆开再相加的动作就是有限元法装配矩阵的全部来源。4.2 单元刚度矩阵与总刚度矩阵有限元法的装配逻辑有限元法的执行过程分三步每一步都很规范。第一步剖分。把区间切成若干个小单元单元长度记为 hᵢ。这一步和差分法的网格剖分看起来一样但意义不同差分法的节点是用来替换导数的取样点有限元法的单元是用来计算积分的积分域。第二步选试探函数。在每个单元上用最简单的分段线性函数去近似真解单元内部是一条直线两个端点上的值作为待定系数。这种函数像一个个小帐篷在每个节点上有一个帐篷尖离开该节点就线性降到零。用它们拼起来整条近似曲线就是连续的折线。第三步装配。把近似函数代回泛函每个单元上的积分算出来是一个关于该单元两个端点值的小二次型写成矩阵就是一个 2×2 的小矩阵习惯上叫单元刚度矩阵。各个单元的贡献按节点编号叠加起来就得到整体的总刚度矩阵和总右端向量。对最简单的 −y″f 加上第一类边界条件单元刚度矩阵的形式是 (1/h) 乘上一个 2×2 的矩阵元素为 1、−1、−1、1。装配时相邻两个单元在共享节点上的贡献直接相加因此总刚度矩阵每行只和相邻节点有关依然是三对角的。而且它对称、正定这两个性质让求解变得非常省心。边界条件的处理方式和差分法类似第一类条件把已知值挪到右端同时把对应的首末方程去掉。这是因为边界上的值不是自由变量不能对它提泛函取极小的要求。处理完之后方程组用追赶法解。有限元法真正的优势在于复杂区域。一维区间上它和差分法打平但换到二维的不规则区域——比如一块带孔的板、一个形状古怪的零件截面——差分法那种规整的矩形网格就很难贴合边界而有限元法可以把区域剖成任意形状的三角形单元边界贴合得好这就是它在结构力学和固体力学里长期占主导的原因。代价是程序复杂单元刚度矩阵里的积分一般算不出解析式得用高斯数值积分近似网格生成、节点编号、装配顺序也都是工程量。4.3 打靶法把未知的初始斜率当成一门要调准的炮第三条路最接近直觉。既然边值问题难在两端都被钉住、没法单向往前推那就先松开一端假设左端的斜率是某个未知数 t配上已知的 y(a)α这就凑成了一个完整的初值问题。用现成的初值问题解法一路推到 b得到一个末端值 y(b;t)。这个末端值显然依赖于 t。想要它等于目标值 β就把问题转化成一个求根问题找 t使 F(t)y(b;t)−β0。这个画面很像打靶射手站在固定的起点目标在另一端的固定位置能调的只有枪口角度。先凭经验打一发看落点偏了多少调整角度再打逐步逼近。教材给的第一种做法就是这种试射加线性插值的策略——打两发用两点的连线估计下一个该用的角度再打一发如此往复直到落点误差落在给定精度以内。第二种做法更规范把 F(t) 当成一个普通函数用二分法、牛顿迭代这类第三章讲过的求根方法来处理甚至可以先粗扫一遍 t 的取值范围找到 F 变号的两个端点再用二分法稳稳地收进去。import numpy as np from scipy.integrate import solve_ivp from scipy.optimize import brentq def shooting_method(f, a, b, alpha, beta, t_guess): 打靶法解 y f(x, y, y)边界 y(a)alpha, y(b)beta f : 接收 (x, y, yp) 返回 y 的函数 t_guess : 初始斜率的试探区间 (t_low, t_high)需保证两端 F 变号 # 把二阶方程降阶成一阶方程组s [y, y] def rhs(x, s): return [s[1], f(x, s[0], s[1])] def end_error(t): # 用当前斜率打一发返回末端值与目标值之差 sol solve_ivp(rhs, [a, b], [alpha, t], rtol1e-10, atol1e-12) return sol.y[0, -1] - beta t_star brentq(end_error, t_guess[0], t_guess[1]) sol solve_ivp(rhs, [a, b], [alpha, t_star], dense_outputTrue) return sol代码里 end_error 就是那个 F(t)brentq 负责在给定区间里找零点它要求区间两端函数值异号这也是为什么参数里要先给一个试探区间。这个区间不是随便给的斜率的量级往往和问题的尺度、系数大小有关凭经验估计一个粗略范围先用几个点扫一遍看 F 有没有变号比盲目设一个很大的区间可靠得多。打靶法的优点是省事不用自己推导差分格式也不用装配矩阵一个现成的初值问题求解器加一个求根器就搭起来了特别适合快速验证和原型阶段。它的软肋也很明确——对初值的敏感性。有些边值问题里起点斜率稍微变一点点末端值就会跑出去很远F(t) 的曲线陡得像悬崖插值和二分都会被数值误差折磨得很难受。这种情况下单次打靶会崩得换多重打靶把区间切成若干段每段各自设初始猜测把所有段的连续性条件联立起来一起解。代价是未知量从 1 个变成几十上百个程序复杂度也上来了。5. 五个最容易混的地方边值问题数值解法里的常见误区方法讲完回过头看初学者最容易绊倒的几个点。这些混淆大多不是因为概念本身多难而是因为名字相近、步骤相似脑子里自动把它们归成了一类。5.1 误区一把边值问题当成带边界的初值问题易混点在于两者都叫定解问题都给方程外加条件看起来只是条件的个数和位置不同。之所以容易混是因为教材讲初值问题在前很多人学完一套推进式的解法见到任何微分方程都想从起点往前推。区分的关键是看条件是否集中在同一点。全部条件压在同一个点上那是初值问题可以单向前推解对初值连续依赖条件分散在两端那是边值问题没有起步方向可用而且可能无解、可能多解。判断标准很简单拿到题目先看条件写在几个点上一个点就是初值两个点就是边值别急着套公式。5.2 误区二三类边界条件只记住了名字易混点在于第一类、第二类、第三类这组编号本身不带任何提示信息背下来也记不住谁是谁。之所以容易混是因为很多人只在纸面上见过它们没和具体的物理场景挂上钩。区分靠的是边界上到底规定了什么量。规定函数值本身是第一类端点被钉死规定导数值是第二类端点高度自由但倾斜角固定规定函数值与导数的线性组合是第三类端点被弹性地拉住。落到程序上还有一个更实用的判别法看边界方程里有没有未知的端点值。有就得把它当未知量处理系数矩阵的首末行要改没有直接挪到右端就行。5.3 误区三以为差分法和有限元法只是网格分得不一样易混点在于两种方法都要剖分网格、都会得到一个三对角方程组、最后都用追赶法解流程看上去高度重合。之所以容易混是因为表面步骤相似掩盖了推导起点的根本差异。区分要往上游看差分法直接从微分方程出发用差商替换导数逼近的对象是方程本身有限元法从变分形式出发先改写成一个泛函极小问题再用分段函数去逼近那个极小点逼近的对象是解。这个差异决定了两者的能力边界差分法在规则区域上实现简单、精度分析清楚有限元法能贴合复杂边界、能处理变系数和复杂区域代价是程序量和前置准备工作大得多。5.4 误区四把收敛阶当成误差一定很小易混点在于O(h²) 这个记号看起来像是在承诺精度很多人默认只要格式是二阶的结果就一定靠谱。之所以容易混是因为收敛阶描述的是误差随步长变化的趋势而实际误差还取决于一个常数因子。区分要盯住那句完整的话误差不超过 M·h²M 和真解的高阶导数有关。真解变化剧烈、导数很大的时候M 可能大到让二阶格式在小步长下依然不准。更别忘了舍入误差步长小到一定程度h² 的截断误差继续下降但浮点运算的舍入误差开始累积上升两条曲线交叉之后再加密网格反而更差。实践中的做法是取几个不同步长各算一遍看结果是否稳定收敛而不是一路把网格加密到底。5.5 误区五以为打靶法是一种新的微分方程解法易混点在于打靶法出现在边值问题这一章名字也新很多人下意识把它当成和差分法、有限元法并列的第三种离散格式。之所以容易混是因为教材把它放在同一层级介绍。区分要看它内部到底调用了什么打靶法自己没有发明任何新的离散手段它把边值问题转成初值问题之后用的还是初值问题那一章的老办法外面再套一个求根器。它更像是组装而不是发明。明白了这一点遇到打靶法不收敛时就知道该往哪查——要么是初值问题求解器的精度不够要么是求根区间没选对而不是去怀疑某个不存在的打靶格式。6. 什么时候该用哪种方法一份边界清单和我的踩坑记录三种方法没有绝对优劣只有场景匹配。把选择依据压缩成一张表比记一堆形容词有用。方法适合的场景边界条件处理主要代价差分法区间规则、系数光滑、需要快速得到二阶精度结果第一类最省事第二三类需改首末行复杂区域难贴合变步长要重新推导有限元法复杂区域、变系数、需要处理不规则边界三类条件都能自然纳入弱形式程序量大需要网格生成和数值积分打靶法快速原型验证、已有成熟的初值问题求解器只适合能稳定前推的问题对初值敏感病态问题时需要多重打靶具体到判断上可以先问三个问题。区间是不是一维且规则是的话差分法通常最划算代码短、精度分析清楚、调试成本低。区域是不是复杂、边界是不是不规则、系数是不是随位置剧烈变化是的话有限元法的前期投入能换回后期的稳定。只是想知道大致形状、想快速验证一个模型跑不跑得通打靶法半小时就能搭出原型但它给出的结果最好只当参考别直接拿去写结论。反过来说什么情况下别用也值得说清楚。差分法不适合硬套到复杂几何上把不规则区域强行塞进矩形网格边界附近的误差会大到让整个结果失去意义。有限元法不适合就算一条曲线这种小任务为了一个一维问题搭一套网格生成加装配的框架投入产出比很难看。打靶法不适合本质不稳定的问题末端值对初始斜率敏感到一定程度单次打靶会反复横跳怎么调都收不住这时应该老老实实换方法而不是继续在求根器上打补丁。我自己在这上面吃过一次亏。有个一维边值问题右端项变化很剧烈我图快用差分法步长随手取了个看起来不大的值程序跑出来曲线光滑、端点也对得上看着特别像那么回事。拿去和加密网格的结果一比中间段差了将近一个数量级。问题出在真解在那一带弯得太急M 这个常数被撑得很大二阶格式的名义精度完全没兑现。从那以后我每次写完差分程序都强制做一遍步长减半的对照两次结果的差异落在预期量级上才敢相信屏幕上的那条线。这个习惯救过我很多次也让我对格式阶数和实际精度之间的差距有了切身体会。这份资料把差分法、有限元法、打靶法三条路线连同极值原理、等价性定理、单元刚度矩阵的推导和对应的程序实现都收在了一起从公式到代码是一条完整的链。遇到边值问题不确定从哪下手的时候翻一翻对照着看比零散地在网上拼凑要省力得多。希望帮到你。本文还有配套的精品资源点击获取
