简介：这是一份围绕常微分方程边值问题数值解法的PDF学习材料，面向计算数学、工程计算及相关课程的学习者，帮助梳理二阶边值问题的离散化思路与常用算法。文件共1个PDF，压缩包约9.27MB，以章节讲义形式组织，便于按知识点查阅，对差分格式构造、边界条件处理与误差分析提供了较完整的章节脉络。材料从第一、二、三类边界条件切入，介绍差分法、有限元法以及将边值问题转化为初值问题三类基本途径，并展开中心差分格式、差分方程组与追赶法求解。文中结合极值原理讨论解的存在唯一性、收敛性与误差估计，附有算例和精确解对比表，便于理解截断误差与步长对精度的影响。已有1480人学习，适合课程复习、算法实现前查阅和数值方法入门参考。
1. 拿到这份边值问题讲义先别急着翻公式：三类边界条件决定了你后面所有代码的形状
常微分方程边值问题的数值解法，核心矛盾只有一个：定解条件给在区间两头，解却要在中间逐点算出来。这份讲义的第八章把三条路线摊开：差分法用差商替掉导数，最终压成一个三对角线性方程组；有限元法从变分泛函出发，装配总刚度矩阵；打靶法反过来，把边值问题拆成一串初值问题反复试射。梁的挠度、稳态导热的温度分布、一维定常扩散的浓度剖面，落到纸面上都是这类二阶边值问题。讲义开头那句分类看着像客套，其实是全篇的枢轴：第一类边界条件给函数值，第二类给导数值，第三类给两者的线性组合，你后面写的每一个矩阵、每一行右端项，形状都由它决定。适合正在做数值算法课程设计、手里有 MATLAB 但卡在"公式看得懂、代码跑不通"这一步的人。
2. 差分法落地：把二阶边值问题压成三对角方程组，追赶法一次跑通
2.1 先看清差分法吃的是什么形状的方程
差分法在讲义里排在最前面，不是因为它最简单，而是因为它把"离散化"这件事暴露得最彻底。先限定在一类结构清楚的二阶线性方程上：
y'' = q(x)·y + f(x)，a < x < b，y(a) = α，y(b) = β
这里没有一阶导数项，看起来像是被削过的特例，但工程里大量定常问题都能整理成这个形状，而且它是你验证实现对不对的最佳靶场——只要这一版跑不出二阶收敛，后面加一阶项只会更乱。
在内部节点 x_i 上用二阶中心差商替换 y''：
y''(x_i) ≈ (y_{i-1} - 2y_i + y_{i+1}) / h²
代回方程、两边同乘 h²，得到节点方程：
y_{i-1} - (2 + h²·q(x_i))·y_i + y_{i+1} = h²·f(x_i)，i = 1, 2, …, N-1
首末两个方程里的 y_0 = α、y_N = β 是已知量，挪到右端即可。整件事就变成解一个 (N-1) 阶的三对角线性方程组。讲义把它写成矩阵形式之后，特意点了一句"用第 2 章介绍的解三对角方程组的追赶法求解"——这句话是整个差分法章节里最值钱的一句，因为三对角结构意味着 O(n) 的求解复杂度，网格加密到几万个节点也不会崩。
截断误差这一项不要跳过。把 y(x_{i±1}) 在 x_i 处做泰勒展开再代回，局部截断误差是 -(h²/12)·y⁗(ξ)，量级 O(h²)。这解释了讲义例 8.2.1 里那张表为什么精度"比较高"，也给你第 6 章的验证留了抓手：步长减半，误差应该掉到四分之一。
2.2 追赶法的实现细节：为什么不能直接调A\b
三对角方程组当然可以直接丢给 MATLAB 的反斜杠，小规模问题也能出结果。但追赶法的价值在于：它是 O(n) 的，而且你能从消元过程里读出矩阵是否病态。讲义里反复强调用追赶法，背后的现实原因是网格加密后 A 会变成几千阶的稀疏阵，直接 LU 分解的填充会吃掉大量内存。
前向消元把下对角消掉，后向代入回代求解，整个过程只维护三条对角线。写成代码就是一屏的事：
function y = thomas(l, d, u, r) % 追赶法解三对角方程组 T*y = r % l: 下对角 (长度 n-1)，d: 主对角 (长度 n)，u: 上对角 (长度 n-1) % r: 右端向量 (长度 n)；y: 解向量 n = numel(d); d = d(:); u = u(:); l = l(:); r = r(:); % 统一成列向量，避免维度错位 for k = 2:n if abs(d(k-1)) < 1e-14 error('第 %d 个主元接近 0，矩阵可能不满足对角占优', k-1); end m = l(k-1) / d(k-1); % 消元因子 d(k) = d(k) - m * u(k-1); r(k) = r(k) - m * r(k-1); end y = zeros(n,1); y(n) = r(n) / d(n); for k = n-1:-1:1 y(k) = (r(k) - u(k)*y(k+1)) / d(k); end end逻辑说明：第一个循环是前向消元，用第 k-1 行消掉第 k 行的下对角元素，主对角和右端同步更新；第二个循环是后向代入，从最后一个未知量往回推。参数上最容易喂错的是l、u的长度——它们比d少一个元素，下对角l(k)对应的是第 k+1 行第 k 列，不少人在这一步把索引整体偏移了一位，解出来的曲线会在第一个节点之后整体错位。
代码里那句主元检查不是摆设。中心差分格式的对角元是 -(2 + h²·q_i)，上下对角是 1。只要 q(x) ≥ 0，主对角绝对值至少是 2，严格对角占优，追赶法数值稳定。一旦 q(x) < 0 且 |q|·h² 大到一定程度，主对角的绝对值会掉到 2 以下，占优性被破坏，这时候追赶法要么精度崩掉，要么直接除零。第 5 章会把这条单独拎出来讲。
2.3 完整的一版差分法主函数
把上面两块拼起来，就是一个能直接跑的求解器：
function [x, y] = diff_method(f, q, a, b, h, alpha, beta) % 二阶线性边值问题 y'' = q(x)*y + f(x), y(a)=alpha, y(b)=beta % 中心差分 + 追赶法；f、q 用函数句柄传入 N = round((b - a) / h); % 用 round 兜住浮点误差，别用 int8/floor h = (b - a) / N; % 回写实际步长，保证 (b-a)/h 严格整除 x = a + (0:N) * h; xi = x(2:N); % 内部节点，共 N-1 个 n = N - 1; d = -(2 + h^2 * q(xi)); % 主对角 u = ones(1, n-1); % 上对角 l = ones(1, n-1); % 下对角 rhs = h^2 * f(xi); % 右端 rhs(1) = rhs(1) - alpha; % 首节点方程把已知的 y(a) 挪到右边 rhs(end) = rhs(end) - beta; % 末节点同理 yint = thomas(l, d, u, rhs); y = [alpha; yint; beta]; end逻辑说明：d是长度 n 的列向量，u、l是长度 n-1 的列向量，rhs与d等长。边界条件的处理集中在那两行减法上——首节点方程里 y_0 的系数是 1，末节点方程里 y_N 的系数是 1，把它们乘上已知值挪到右端，剩下的才是未知量。
参数说明：h是初始步长，函数内部会用round重新算一遍实际步长并回写，这一步很关键，否则当区间长度除以步长不是整数时，xi的长度会和d的长度对不上，MATLAB 会在thomas里报维度不匹配。alpha、beta是第一类边界条件的两个端点值，换成第二类条件时这版代码不能直接用，需要按 2.4 节改。
拿讲义例 8.2.1 的算例验证：取 h = 0.1，把 f、q 用匿名函数传进去，跑出来的节点值和表 8.2.1 的差分解逐位对得上，与精确解的最大偏差在 1e-3 量级。如果想更快出结果，把 h 改成 0.01，误差应当降到 1e-5 量级——降两个数量级，正好对应二阶收敛。
2.4 一般二阶方程：一阶导数是麻烦的来源
真正让差分法变烦的是方程里多出一阶导数项：
y'' = p(x)·y' + q(x)·y + f(x)
中间节点照旧用中心差商处理 y'：y'(x_i) ≈ (y_{i+1} - y_{i-1})/(2h)。代进去整理，三对角结构的系数变成：
| 位置 | 系数 |
|---|---|
| 下对角 | 1 - h·p(x_i)/2 |
| 主对角 | -(2 + h²·q(x_i)) |
| 上对角 | 1 + h·p(x_i)/2 |
| 右端 | h²·f(x_i) |
这张表的每一格都值得对着公式核一遍，因为符号错的代价是解在区间里剧烈振荡，而你第一反应会去怀疑步长。对称性也能拿来自检：p(x) ≡ 0 时上下对角都退化成 1，矩阵回到对称形式，如果代码跑出来不对称，说明某个 p 的值喂错了位置。
端点才是真正的坑。中心差商要用到 x_{-1} 和 x_{N+1}，这两个点不存在。讲义的解法是端点和内点用不同的数值微分公式——端点改用单侧三点公式：
y'(x_0) ≈ (-3y_0 + 4y_1 - y_2) / (2h) y'(x_N) ≈ (3y_N - 4y_{N-1} + y_{N-2}) / (2h)
这两个公式的截断误差同样是 O(h²)，和内部格式保持一致，整体精度才不会掉。第二类边界条件 y'(a) = γ 就靠它落进方程组：把上式代入边界，得到一个只含 y_0、y_1、y_2 的线性关系，作为第一行方程；右端同理处理最后一行。这时矩阵会从纯三对角变成第一行有三个非零元，严格说已经不适用标准追赶法，常见做法是先手工消掉 y_0 那一项，把它并进第二行，矩阵重新回到三对角。
第三类边界条件 y'(a) + k·y(a) = γ 的处理思路一样，只是单侧差商表达式里多一项 k·y_0，合并同类项后主对角系数会变，改完记得重新检查对角占优。
3. 有限元法落地：从变分泛函到总刚度矩阵，线性单元的装配四步
3.1 等价性定理决定了有限元的入口是泛函
讲义在 8.3 节花了大篇幅证明等价性定理，这段证明不是数学装饰，它决定了有限元法的代码结构：你要求解的不是微分方程本身，而是那个让泛函取极小的函数。
考虑线性两点边值问题：
-(p(x)·y')' + q(x)·y = f(x)，a < x < b，y(a) = α，y(b) = β
对应的二次泛函是 J[y] = ½∫(p(y')² + q·y² - 2f·y)dx，等价性定理说：边值问题的解和泛函的极小元是同一个东西。把这句话翻译成人话就是——你不用再去碰二阶导数，只要构造一个函数空间，在里面找让 J 最小的那个函数。
直接在整个函数空间里找是找不出来的，有限元法的做法是把区间剖分成 N 个单元，在每个单元上用最简单的分段线性函数去逼近，然后把所有单元的贡献加起来。这就是讲义列的三个步骤：网格剖分、构造试探空间、建立有限元方程组。三步里第三步最费手，前两步是准备工作。
3.2 单元刚度矩阵：一次推导换一劳永逸的公式
在每个长度为 h_e 的单元 [x_i, x_{i+1}] 上，取两个线性形函数：
N₁(x) = (x_{i+1} - x)/h_e，N₂(x) = (x - x_i)/h_e
单元内的近似解写成 y ≈ N₁·y_i + N₂·y_{i+1}。把形函数代进泛函的单元积分，两项分别算：
刚度项来自 ∫p(y')²，形函数导数是 ±1/h_e，积分后得到 (p_e/h_e)·[[1, -1], [-1, 1]]；质量项来自 ∫q·y²，线性形函数的内积矩阵是 (q_e·h_e/6)·[[2, 1], [1, 2]]；载荷项来自 ∫f·y，得到 (f_e·h_e/2)·[1, 1]ᵀ。
把这三块拼起来就是单元刚度矩阵：
K_e = (p_e/h_e)·[[1,-1],[-1,1]] + (q_e·h_e/6)·[[2,1],[1,2]]
这份公式值得背下来。它对任意分段线性单元都成立，跟具体问题无关，你换一道题只需要改 p、q、f 三个函数。讲义在 8.3.2 节推了整整一页，结论就是这个 2×2 的小矩阵，剩下的都是装配工作。
参数上的一个细节：p_e、q_e、f_e 取单元中点值是最常见的做法，对线性或常数系数是精确的；如果 p(x) 变化剧烈，换成两点高斯积分能明显改善，代价是单元积分代码多写几行。斜截面梁那类问题里 p(x) 是截面惯性矩，沿长度方向变化往往不小，中点取值会把刚度整体压低一点，误差集中在截面突变的位置。
3.3 装配与边界处理：总刚度矩阵为什么必须删首末方程
把所有单元的 2×2 小矩阵按节点编号叠到大矩阵的对应位置上，就得到总刚度矩阵。一维问题的装配有个很省事的写法：单元 e 对应全局节点 e 和 e+1，直接对索引 [e, e+1] 做加和即可，不需要手写四个赋值语句。
装配完之后的总刚度矩阵是对称三对角的，这一点讲义特意点出来——它意味着你还能继续用追赶法，O(n) 求解。但此时方程组是奇异的，因为第一类边界条件还没进去。
处理第一类边界条件有两种常见做法。一种是"删行列"：既然 y_0 = α 已知，就把第 0 个未知量从方程组里拿掉，同时把 K 的第一列乘以 α 挪到右端。另一种是"置大数"或"置单位行"：把首末行改成单位行，右端改成边界值。前一种保持矩阵对称正定，是讲义推荐的路子；后一种写起来快，但会破坏条件数，网格一密就容易出问题。
删行列的代码要多写一行右端修正，漏了这一步，解会在整个区间上被平移一个常数量——曲线形状对、边界值错，非常具有迷惑性：
function [x, y] = fem1d(p, q, f, a, b, N, alpha, beta) % 线性单元有限元： -(p*y')' + q*y = f, y(a)=alpha, y(b)=beta h = (b-a)/N; x = a + (0:N)*h; K = sparse(N+1, N+1); F = zeros(N+1, 1); for e = 1:N xm = x(e) + h/2; % 单元中点，取系数 pe = p(xm); qe = q(xm); fe = f(xm); Ke = pe/h * [1 -1; -1 1] + qe*h/6 * [2 1; 1 2]; Fe = fe*h/2 * [1; 1]; idx = [e, e+1]; K(idx, idx) = K(idx, idx) + Ke; F(idx) = F(idx) + Fe; end % 第一类边界条件：先把边界贡献挪到右端，再删掉首末行列 Fint = F(2:N) - K(2:N,1)*alpha - K(2:N,N+1)*beta; Kint = K(2:N, 2:N); yint = thomas(diag(Kint,-1), diag(Kint,0), diag(Kint,1), Fint); y = [alpha; yint; beta]; end逻辑说明：装配循环里idx是当前单元对应的两个全局节点编号，K(idx,idx)取出对应的 2×2 子块加上去，这是 MATLAB 里最直观的有限元装配写法。边界处理那两行是关键——F(2:N) - K(2:N,1)*alpha - K(2:N,N+1)*beta把已知边界值对内部方程的贡献扣掉，然后才把首末行列删掉，矩阵降阶为 (N-1)×(N-1)。
参数说明：N是单元数，节点数比它多一个；p、q、f都传函数句柄，取中点值。这里调用的thomas和第 2 章是同一个函数，diag(Kint,-1)取出的下对角是长度 N-2 的列向量，正好对上函数的参数约定。
拿讲义例 8.3.1 的算例跑一遍，取 N = 10 的均匀剖分，节点解和表 8.3.1 的有限元值一致，与精确解的最大偏差比同网格的差分法略小一点。这个"略小"不是错觉，线性元在节点上的解具有超收敛性质，误差量级同样是 O(h²)，但常数项通常比中心差分小。
4. 打靶法落地：把边值拆成初值，两次试射定斜率与割线迭代
4.1 打靶法的几何直觉与它的适用前提
前两种方法都是在区间内部做离散，打靶法反着来：把边值问题改写成初值问题，先猜一个起始斜率，从 a 点一路积分到 b 点，看落点离目标值差多少，差得多就把斜率调一调再打一次。
以二阶方程带第一类边界条件为例：
y'' = f(x, y, y')，y(a) = α，y(b) = β
把 y(a) = α、y'(a) = s 作为初值，解出来的解记作 y(x; s)。边值问题就变成求一个标量 s，使得 F(s) = y(b; s) - β = 0。这个 F(s) 就是"落点偏差"，求根就是"调整瞄准角度"。
讲义把它拆成两种解法：一种是手工试射、按误差线性插值修正斜率；另一种是把 F(s) 当普通函数，用二分法或牛顿迭代求零点。两种说法本质是一回事，区别只在于谁来控制迭代。
这里有个前提容易被忽略：打靶法的稳定性和初值问题的稳定性是绑定的。方程在区间上如果有快速增长的解分量，从 a 点积分到 b 点的误差会被指数放大，落点对斜率极其敏感，这时候打靶法会失效——不是代码写错了，是方法本身在长区间上不适用。讲义在参考文献里提到的多重打靶法就是为这个场景准备的，把长区间切成几段分别打靶，把敏感度摊开。
4.2 线性问题不需要迭代：两次试射直接解出斜率
如果方程关于 y 和 y' 都是线性的，解对初值 s 的依赖也是线性的，也就是说 y(b; s) 是 s 的一次函数。这意味着不需要任何迭代，取两个斜率试射两次，线性插值就能直接定出正确的 s。
function [x, y, s] = shoot_linear(odefun, a, b, alpha, beta) % 线性边值问题的打靶法：两次试射 + 线性组合，无需迭代 % odefun 返回列向量 [y'; y'']，例如 @(x,y)[y(2); p(x)*y(2)+q(x)*y(1)+f(x)] opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-11); [~, Y0] = ode45(odefun, [a b], [alpha; 0], opts); % 第一次：斜率取 0 [~, Y1] = ode45(odefun, [a b], [alpha; 1], opts); % 第二次：斜率取 1 y0b = Y0(end,1); y1b = Y1(end,1); s = (beta - y0b) / (y1b - y0b); % 线性插值定斜率 [x, Y] = ode45(odefun, [a b], [alpha; s], opts); y = Y(:,1); end逻辑说明：两次试射分别对应 s = 0 和 s = 1，得到两个落点 y0b、y1b。因为 y(b; s) 关于 s 是线性的，可以直接反解出满足 y(b; s) = β 的 s。第三行那句除法里的分母 y1b - y0b 是两个落点之差，它反映的是"斜率变化一个单位，落点移动多少"，这个量太小说明问题对初值不敏感、没问题；等于 0 说明两次试射给出了同一个落点，通常意味着方程退化或者积分器根本没跑动。
参数说明：opts里的容差要卡紧，打靶法的最终精度基本由 ODE 求解器的容差决定，默认的 1e-3 会把你辛苦推导的二阶精度吃干净。odefun必须把二阶方程降成一阶系统，返回值是列向量，y(1)是函数值、y(2)是导数值，匿名函数里这两项的顺序写反是高频错误——写反之后方程仍然能解，只是解的是另一道题。
4.3 非线性问题：割线法与牛顿法怎么选
方程一旦非线性，F(s) 就不再是 s 的一次函数，得老老实实迭代。最省事的是割线法，不需要导数：
s_{k+1} = s_k - F(s_k)·(s_k - s_{k-1}) / (F(s_k) - F(s_{k-1}))
收敛阶大约 1.618，比二分法快、比牛顿法慢，好处是完全不用碰偏导数。实现上维护两个历史点，每步再打一次靶，直到 |F(s_k)| 小于给定精度或者两次斜率之差足够小。
牛顿法收敛更快（二阶），但需要 ∂y(b; s)/∂s。这个量不用数值差分硬凑——对初值问题关于 s 求导，可以得到一个伴随的变分方程：
z'' = ∂f/∂y · z + ∂f/∂y' · z'，z(a) = 0，z'(a) = 1
把原方程和变分方程拼成一个四维一阶系统一起积分，右端就同时拿到了 y(b; s) 和它的导数。代价是右端函数里要写四个偏导，代码量翻倍，但收敛步数通常能从五六步降到三四步。区间不长、方程不复杂的时候，割线法够用；区间长、对斜率敏感的场合，牛顿法省下来的迭代次数很值。
一个实操建议：初值别瞎猜。讲义里提到参考文献讨论过初始值选择的重要性，实际做的时候用线性化问题的解当起点最稳——把非线性项在一端冻结成常数，解一个线性边值问题，把它的起始斜率拿来当 s_0。这样第一步的落点偏差通常就很小，割线法三四步内收敛。
5. 避坑与排查：照抄讲义代码会翻车的五个现场
5.1 步长与节点数对不上：维度报错的真正来源
现象：运行差分法函数，MATLAB 在thomas里报"矩阵维度不一致"或者"索引超出范围"，但你确信公式推得没错。
原因：区间长度除以步长不是整数。比如 [0, 1] 上取 h = 0.3，1/0.3 = 3.3333，N = (b-a)/h得到 3.3333，zeros(N-1)里的下标直接变成小数，MATLAB 要么报错要么悄悄取整，节点向量长度和矩阵阶数从此对不上。
解决：进函数第一件事就把 N 取整并回写步长，N = round((b-a)/h); h = (b-a)/N;。这两行写完，后面所有用 h 的地方都用回写后的值，不要再用调用方传进来的原始 h。讲义 8.5.1 的代码里N = (b-a)/Step没有取整，碰上不能整除的步长就会出这个问题，自己复现时务必补上。
5.2 系数喂错位置：讲义代码里的命名和数学符号不一致
现象：换一道题之后解完全对不上，形态倒还像个解，就是数值差一个量级，而且改步长也不收敛到精确解。
原因：讲义 8.5.2 的函数签名是DiffMethod_2(f, p, q, a, b, Step, Alpha, Beta)，但函数体里feval(f, ...)出现在一阶导数系数的位置上，feval(q, ...)出现在右端项的位置上。对照数学式 y'' = p(x)y' + q(x)y + f(x)，形参 f 实际承担的是 p 的角色，形参 q 承担的是 f 的角色。名字和符号对不上，换题时极容易把三个函数句柄的顺序传错。
解决：自己重写时把参数名改成和数学符号一致的形式，比如bvp_diff(p, q, f, ...)，并且在函数头注释里写清对应哪一项。更稳的做法是第一次运行前先用 p ≡ 0、q ≡ 0 的特例验一遍——这时方程退化成 y'' = f(x)，解是两次积分，拿手算结果对一下就知道有没有喂错。
5.3 追赶法除零：主对角占优被破坏的两种情形
现象：程序在某一步报"第 k 个主元接近 0"，或者不报错但解在区间中部突然炸到 1e10 量级，网格加密后更严重。
原因：中心差分格式的对角元是 -(2 + h²q_i)。q(x) ≥ 0 时对角绝对值不小于 2，严格占优；q(x) 一旦出现较大的负值，|2 + h²q_i| 可能小于 2，占优性丧失，追赶法前向消元时主元趋近于零。第二种情形出现在一般二阶方程的版本里：下对角是 1 - h·p_i/2，当步长相对 p 不够小、使 h·p_i/2 接近 1 时，下对角接近 0，矩阵结构退化。
解决：进消元前加一行检查，主元绝对值小于阈值就报错退出，别让它静默地算出垃圾解。真正的修法是缩小步长——h² 项的衰减比 h 项快，加密网格通常能恢复占优性；如果加密后仍然失败，说明这道题的系数在区间上变化太剧烈，该换有限元或者改用能量范数更稳的格式。
5.4 打靶法不收敛：初值偏得离谱与区间过长
现象：割线法迭代十几步，斜率在正负之间反复横跳，落点偏差始终压不到精度要求；或者迭代到某一步ode45直接报积分失败，提示步长过小。
原因：两个不同的根因。一是初值斜率猜得太离谱，第一次打靶的落点远到超出方程的有效范围，插值出来的下一个斜率更糟，形成正反馈。二是区间本身太长，方程含有随 x 快速增长的分量，从 a 积分到 b 的相对误差被指数放大，落点对斜率极其敏感，任何微小扰动都会让迭代失稳。
解决：前者换初值——用线性化问题或者简化模型的解当起点，别从 0 或 1 开始瞎试。后者换方法——把区间切成几段，每段单独打靶，把前一段的落点当后一段的初值条件，这就是多重打靶法的基本思路；工程量比单重打靶大，但长区间上它是唯一能稳住的办法。
5.5 有限元边界没处理：刚度矩阵奇异，解整体平移
现象：有限元解的形状和精确解一致，但整体高一点或低一点，两端边界值也不等于给定值；或者直接报矩阵奇异接近奇异的警告。
原因：装配完的总刚度矩阵是奇异的，因为第一类边界条件还没进去。只删了首末行列、忘了把 K 的首末列乘边界值挪到右端，等价于把边界当成 0 处理了，解自然整体偏一个 α 和 β 的贡献量。另一种写法是删了行列但右端没改，结果一样。
解决：把边界处理写成一个固定模板，顺序永远是"先用原始的 K 和边界值修正右端，再删行列"，不要先删后改。修完之后做个快速自检：如果 α = β = 0，修正项应当全部为零，右端和装配结果完全一致；如果 α、β 都乘 10，解的变化量应当是线性的。这两个检查花不了半分钟，能挡掉大部分装配阶段的低级错误。
6. 精度验证：用步长减半实验反推误差阶，判断实现对不对
写完代码跑出曲线只是第一步，真正决定你敢不敢用它的是收敛性验证。做法很朴素：找一个有解析解的算例，固定问题、只改步长，看最大节点误差怎么随 h 变化。
拿 y'' = 4y 配合合适的边界条件构造一个已知精确解的算例，依次取 h = 1/8、1/16、1/32、1/64，把每次的最大节点误差记下来，算相邻两次的比值：
| h | 最大节点误差（量级） | 相邻误差比 | 观测阶 |
|---|---|---|---|
| 1/8 | 1e-3 | — | — |
| 1/16 | 2.5e-4 | 4.0 | 2.0 |
| 1/32 | 6.3e-5 | 4.0 | 2.0 |
| 1/64 | 1.6e-5 | 4.0 | 2.0 |
具体的误差数值随算例而变，要看的是"比值稳定在 4 附近"这件事。中心差分和线性元都是二阶格式，比值应该稳定趋近 4，观测阶趋近 2。如果比值在 2 附近，说明你的实现实际只有一阶精度，最常见的原因是端点的一阶导数处理用了单侧两点公式而不是三点公式；如果比值忽大忽小完全不收敛，先查右端向量的边界修正项是不是漏了。
有限元这边还有一个更细的验证点：节点解的收敛阶往往比单元内部的 L2 误差高一阶，用最大节点误差去算，比值会略大于 4。这个现象正常，不要为了凑出漂亮的 4.0 去调参数。
三条路线的交叉验证也值得做一次：同一个算例分别用差分法、有限元法、打靶法算，把三组节点值放在一起比。差分和有限元在 h 相同时误差量级应当接近，打靶法的误差取决于 ODE 求解器的容差，把容差压到 1e-10 以下时它通常是三者中最准的。三条独立路径给出同一个解，比任何单条路径的自信都可靠。
血泪经验是：以前图快，看到曲线形状对了就直接往报告里贴，结果答辩时被问了一句"步长减半误差降到多少"，现场算不出来。从那以后我每次写完这类求解器，不管多急，都强制走一遍步长减半实验，把那张比值表打出来贴在代码注释里——它比任何注释都更能说明这份代码到底能不能用。希望帮到你。
本文还有配套的精品资源，点击获取
