简介面向学习数值分析与科学计算的高校学生这份东南大学出版配套上机题文档围绕MATLAB编程实践展开聚焦三次样条插值函数、重积分计算和常微分方程初值问题数值解三类典型题目。文档不仅给出第一型三次样条插值、复化梯形公式与外推、RK4、AB4及AB4-AM4预测校正等方法的算法思路还结合车门曲线型值点、二重积分和给定初值问题展示计算过程与结果对比并附有精度排序和编程反思便于读者理解公式如何落地为程序、误差从何而来。文档同时强调主程序与子函数分离的组织方式提到逐次二分步长、外推提高精度以及公式输入错误引发的排错经历对课程作业自查和MATLAB实现复盘较有参考价值。资源包仅含1个docx文件约566KB篇幅紧凑适合上机前预习与课后对照。目前已有904人学习可作为数值分析上机训练的实用样例。1. 三道上机题卡人的从来不是公式数值分析上机题最容易给人错觉三次样条的三弯矩方程、复化梯形加外推、RK4 和 AB4-AM4 的系数公式都在书上印着抄下来似乎就完事。真坐到 MATLAB 前面翻车的往往不是推导能力而是把 μ 和 λ 的位置写反、把h/24抄成h/12、把内层外层积分变量搞混——这些地方不报错只给你一条看着挺光滑、其实整体偏了的曲线。这份数值分析上机题资料里连着做了三道题三次样条插值车门曲线、二重积分复化梯形逐次二分加外推、常微分方程初值问题RK4、AB4、AB4-AM4 和改进版对比。下面按这三道题的顺序拆把每道题的算法落点、可抄的 MATLAB 骨架、以及我实测过会踩的坑都摆出来。新手可以照着把程序跑通熟手能对一下判断依据做课程作业的、临时要补 MATLAB 数值程序的都能直接拿走。2. 三次样条插值三弯矩方程怎么补边界、怎么解三对角三次样条这道题的价值不在结果在于它把“公式推导 → 线性方程组 → 程序实现 → 结果回代”四件事完整走了一遍。第一型给定两端一阶导数的三次样条本质就是解一个三对角线性方程组求二阶导数 M再用 M 把每个小区间上的三次多项式写出来。思路清楚了代码不到五十行思路没理清改一天也调不出来。2.1 型值点先核对一遍再谈公式车门曲线的型值点一共 11 个x 从 0 到 10 等距步长 h1。资料里那行数据在复制、截图、转 Word 的过程中很容易粘连2.51和3.30挨在一起一不留神就变成13.30而后面 4.04、4.70、5.22 一路平滑上去唯独第二个点跳到 13插值出来必然是一段剧烈振荡。我一般会先把数据摆成表肉眼扫一遍量级再往程序里灌。i012345678910x_i012345678910y_i2.513.304.044.705.225.545.785.405.575.705.80边界条件是两端一阶导数已知左端 y(0)0.8右端 y(10)0.2。这两个数是第一型Clamped条件的全部输入少了任何一个方程组就欠两行M 解不出来。等距情况下 h_j 全等于 1这件事后面能用来手算校验别浪费。2.2 三弯矩方程的三个系数缺一个都解不出对内部节点 j1,…,n-1三弯矩方程是μ_j·M_{j-1} 2·M_j λ_j·M_{j1} d_j其中 h_j x_{j1} - x_jμ_j h_{j-1}/(h_{j-1}h_j)λ_j h_j/(h_{j-1}h_j) 1 - μ_j右端 d_j 6·f[x_{j-1}, x_j, x_{j1}]也就是 6 倍二阶差商。资料里把 λ 写成了 α这是符号习惯问题不影响计算但要清楚它对应的是 M_{j1} 的系数不是 M_{j-1} 的。μ 和 λ 一个挂在左邻、一个挂在右邻写反了方程组依然能解、程序依然不报错只是解出来的 M 序列会整体错位画出来的曲线在节点处一阶导不连续——肉眼不细看根本发现不了这是这道题最阴的地方。二阶差商我习惯显式写出来避免嵌套调用出错f[x_{j-1}, x_j, x_{j1}] [ (y_{j1}-y_j)/h_j − (y_j-y_{j-1})/h_{j-1} ] / (h_{j-1}h_j)等距、h1 的时候这个式子退化成 (y_{j1} - 2y_j y_{j-1})/2可以直接口算几个 d_j 来对程序输出。这是我最常用的一招先手算 j1 那一行的 μ、λ、d 三个数跟程序打印的第一行对一遍对不上就别往下跑了越跑越乱。2.3 第一型边界条件怎么补两行方程内部的 n-1 个方程不够一共 n1 个未知量 M_0 到 M_n还差两行由端点导数条件补。第一型的补法是2·M_0 M_1 (6/h_0)·( f[x_0,x_1] − y_0 )M_{n-1} 2·M_n (6/h_{n-1})·( yn − f[x{n-1},x_n] )注意两行里差商和导数的先后顺序是反的左端是“差商减导数”右端是“导数减差商”。这个顺序反了误差不会暴露成崩溃只会让整条曲线两端微微翘起。第二型给定二阶导数的补法不一样是 2M_0 M_1 6·f[x_0,x_1,x_2] 这种形式第三型周期性更不一样。三类型别在同一个程序里混用是最常见的翻车点之一做题前先确认题目给的是哪一型。拼起来就是一个 (n1)×(n1) 的三对角方程组 A·M d直接用 MATLAB 的反斜杠A\d能解但既然题目要求“编制通用程序”把追赶法Thomas 算法写出来更符合要求而且 O(n) 的复杂度在节点多的时候优势明显。2.4 通用程序的骨架构造—追赶—回代先写主函数负责组装三对角矩阵、调用追赶法、返回 Mfunction [M, h] spline3_clamped(x, y, dy0, dyn) % 第一型三次样条输入型值点 x,y 与两端一阶导数 dy0,dyn % 输出各节点二阶导数 M长度与 x 相同和区间步长 h n length(x) - 1; % 子区间个数 h diff(x); % h(j) x(j1) - x(j) A zeros(n1); d zeros(n1,1); % 左边界2*M0 M1 6/h0 * (f[x0,x1] - y0) A(1,1) 2; A(1,2) 1; d(1) 6/h(1) * ((y(2)-y(1))/h(1) - dy0); % 内部节点mu*M(j-1) 2*M(j) lam*M(j1) 6*f[x(j-1),x(j),x(j1)] for j 2:n mu h(j-1) / (h(j-1) h(j)); lam 1 - mu; A(j,j-1) mu; A(j,j) 2; A(j,j1) lam; d(j) 6 * ((y(j1)-y(j))/h(j) - (y(j)-y(j-1))/h(j-1)) ... / (h(j-1) h(j)); end % 右边界M(n-1) 2*Mn 6/h(n) * (yn - f[x(n-1),x(n)]) A(n1,n) 1; A(n1,n1) 2; d(n1) 6/h(n) * (dyn - (y(n1)-y(n))/h(n)); M thomas(A, d); % 追赶法求解见下 end追赶法把三对角矩阵压成三个向量前代消元、回代求解注意第一行和最后一行的消元边界function x thomas(A, b) % 三对角方程组追赶法A 为方阵只取三条对角线 n length(b); a diag(A,-1); bb diag(A,0); c diag(A,1); for i 2:n m a(i-1) / bb(i-1); % 消元因子 bb(i) bb(i) - m * c(i-1); b(i) b(i) - m * b(i-1); end x zeros(n,1); x(n) b(n) / bb(n); % 回代 for i n-1:-1:1 x(i) (b(i) - c(i)*x(i1)) / bb(i); end enddiag(A,-1)取的是下标 (i, i-1)对应 a(i-1)diag(A,1)取的是 (i, i1)对应 c(i)。这两个索引偏移错一位追赶法解出来的结果会差得很离谱但矩阵本身仍然可解不报错。所以写完先拿max(abs(A*x - b))验残差残差在 1e-12 量级才算通过。2.5 打印 S(i0.5) 的定位技巧题目要打印 S(0.5)、S(1.5)……S(9.5)。插值函数是分段定义的得先判断查询点落在哪个子区间。我用find(x xq, 1, last)定位左端点下标再做上下界钳位防止 xq 恰好等于最后一个节点时越界function v eval_spline(x, y, M, xq) % 用二阶导数 M 回代计算三次样条在 xq 处的值 h diff(x); v zeros(size(xq)); for k 1:numel(xq) j find(x xq(k), 1, last); j min(max(j, 1), length(x)-1); % 钳位防止越界 hj h(j); t1 (x(j1) - xq(k)) / hj; % 左基函数 t2 (xq(k) - x(j)) / hj; % 右基函数 v(k) t1*y(j) t2*y(j1) ... ((t1^3 - t1)*M(j) (t2^3 - t2)*M(j1)) * hj^2 / 6; end end这个式子是三次样条的标准回代公式(t^3 - t)那一项乘的是 M 和 h²/6不是 h/6量纲才对得上。写完之后做两个自检一是eval_spline在节点处的值必须等于 y误差 1e-12 以内二是两端的一阶导数应该逼近 0.8 和 0.2差分数值近似即可。这两条过了S(i0.5) 的结果基本不会错。3. 重积分外推从 Tn 到 Sn 的那一步才是精度来源二重积分这题看着比样条简单其实考点全在“逐次二分 外推”上。题目给的是 I(f) ∫_c^d ( ∫_a^b f(x,y) dx ) dy要求用复化梯形公式、逐次二分步长、外推算到相邻两次近似值之差的绝对值不超过 ε。如果只是把两次梯形值一减就交那外推等于没做精度停在 O(h²)二分十几次也追不上要求。3.1 为什么是 4/3 和 -1/3一次外推把误差阶抬上去复化梯形公式的误差展开是 T(h) I C·h² O(h⁴)。步长减半后 T(h/2) I C·h²/4 O(h⁴)。两式一消把 h² 项干掉S (4/3)·T(h/2) − (1/3)·T(h) I O(h⁴)这就是资料里写的 S_n (4/3)T_{2n} − (1/3)T_n。系数 4/3 和 -1/3 不是拍脑袋来的是消 h² 项消出来的4/3 − 1/3 1 保证常数项对4/3·(1/4) − 1/3 0 保证 h² 项对。写代码时这两个系数写反、或者写成 3/4 和 1/4结果会稳稳地偏离真值而且二分越多次偏离越“稳定”看起来特别像收敛了。理解这一点之后判断依据就很清楚外推后的 S 序列收敛速度应该明显快于 T 序列。如果打印出来发现 S 和 T 收敛得一样慢八成是系数或者 T_{2n}/T_n 的顺序搞反了。3.2 内层 x、外层 y二重循环的顺序不能乱题目写的是 ∫_c^d (∫_a^b f(x,y) dx) dy内层对 x 积、外层对 y 积。程序里两层循环必须严格对应外层循环变量是 y从 c 到 d内层是 x从 a 到 b。我见过不少把两层换过来的写法程序照样出数、照样收敛只是算出来的是另一个积分的值——因为积分区域是矩形的时候交换次序数学上等价但如果被积函数里的 x、y 在调用时传反了比如f(x,y)写成了f(y,x)结果就悄悄错了。被积函数 f(x,y) tan(x² y²) 里 x、y 是平方和交换了看不出差别这类题最容易掩盖参数顺序错误。我的习惯是给临时测试函数f(x,y) x*10 y跑一遍看外层积分变量是不是真的在外面两分钟就能确认比事后调半天强。积分限这一题是 a0、bπ/3、c0、dπ/6上界小、下界大这种写法在题目里也常见(b-a)一减如果是负的整个积分符号就反了。写完先把积分限代进hx (b-a)/nx打一行看看正负。3.3 三个文件的拆分myfun / getT / 主程序资料里把被积函数、复化梯形子函数、主程序分成三个文件这个拆法是对的——被积函数一改就能换题梯形子函数可以单独测。被积函数function z myfun(x, y) % 被积函数注意 x 是内层积分变量、y 是外层 z tan(x.^2 y.^2); end复化梯形二重积分注意权重角点乘 1/4边上乘 1/2内部乘 1。写代码时把权重拆成wx * wy最不容易错function T getT(f, a, b, c, d, nx, ny) % 矩形区域上的复化梯形二重积分nx/ny 为 x/y 方向等分数 hx (b-a)/nx; hy (d-c)/ny; T 0; for i 0:ny y c i*hy; wy 1; if i0 || iny, wy 0.5; end % y 方向端点权重 inner 0; for j 0:nx x a j*hx; wx 1; if j0 || jnx, wx 0.5; end % x 方向端点权重 inner inner wx * f(x, y); end T T wy * inner; end T T * hx * hy; % 最后统一乘步长 endhx*hy放在循环外面统一乘比每次内层都乘一遍快也避免累加过程中浮点误差被放大。这个写法在 nx、ny 都到几千的时候差别能看出来。3.4 二分停止条件写在哪儿决定了会不会死循环题目要求“相邻两次近似值之差的绝对值不超过 ε”ε 取 0.5e-5。关键是这个“近似值”指的是外推后的 S不是原始的 T。判据写成abs(T2 - T1) eps就偏了题意写成abs(S2 - S1) eps才对。eps0 0.5e-5; nx 2; ny 2; Tn getT(myfun, 0, pi/3, 0, pi/6, nx, ny); % 粗网格 nx 2*nx; ny 2*ny; T2n getT(myfun, 0, pi/3, 0, pi/6, nx, ny); Sn 4/3*T2n - 1/3*Tn; % 第一次外推 it 0; while it 20 % 上限兜底防死循环 Tn T2n; nx 2*nx; ny 2*ny; T2n getT(myfun, 0, pi/3, 0, pi/6, nx, ny); S2n 4/3*T2n - 1/3*Tn; if abs(S2n - Sn) eps0 break; end Sn S2n; it it 1; end fprintf(积分值 %.10f二分次数 %d\n, S2n, it);it 20这个上限不是可选项是必需品。浮点数在 1e-5 附近会进入震荡区两个相邻的 S 值可能来回跳永远不满足小于 ε没上限就是死循环。而且二分到 nx、ny 上千以后单次 getT 的耗时是平方级增长的卡死一次能等十几分钟。3.5 运行慢的真因每次二分都在重算旧节点资料里反思“运行速度慢是很大的不足”这个判断是准确的但原因不是“MATLAB 慢”而是每次getT都从头把所有节点重算了一遍。实际上复化梯形公式有递推关系T_{2n} T_n / 2 (h/2) · Σ f(新增节点)新增节点的横坐标是 a (2k-1)h/2数量正好是 n 个。用这个递推每次二分的计算量从 O(4n²) 降到 O(n²) 的一小部分尤其两个方向同时二分的时候节省非常明显。改法是把getT拆成“首算”和“加细”两个函数主循环里只算新增点。做作业不一定要求但把这段写进反思里比写“还需要更深入学习 MATLAB”有说服力得多。4. 常微分方程初值RK4、AB4、AB4-AM4 和改进版的误差差在哪第三道题最有意思因为它把四种方法放在同一个初值问题上横向比。初值问题是 y −x²y²0 ≤ x ≤ 1.5y(0) 3步长 h 0.1精确解 y 3/(1x³)。四种方法、同一组数据、同一个误差范数排序出来一目了然。这道题真正考的不是“会不会写 RK4”而是“知不知道为什么 AB4 会比 RK4 差两个数量级”。4.1 先把精确解和右端函数准备好右端函数和精确解各写一个文件后面所有方法共用function f fxy(x, y) % 一阶导数右端项 f -x.^2 .* y.^2; end function y fx(x) % 精确解用于误差对比 y 3 ./ (1 x.^3); end精确解这个文件很容易被当成“可有可无”实际它是整道题的裁判。没有它四种方法的误差没法排序题目第6问的“能得到哪些结论”就只能凭感觉写。我一般会把精确解和数值解画在同一张图上再加一张误差随 x 变化的图误差曲线的形状比单个数字能说明问题——RK4 的误差曲线几乎是平的AB4 的误差曲线一路往上爬这个差别一眼就能看出来。验证一下精确解和右端项一致对 y 3/(1x³) 求导得 −9x²/(1x³)²而 −x²y² −9x²/(1x³)²两边相等。这一步花三十秒能避免后面所有对比都建立在一个写错的 fx 上。4.2 RK4 是单步自启动AB4 不是RK4 是单步法只依赖 y_i 就能算出 y_{i1}所以第一个点起步毫无障碍function [x, y] rk4(f, x0, xend, y0, h) % 经典四阶 Runge-Kutta x (x0:h:xend); N length(x); y zeros(N,1); y(1) y0; for i 1:N-1 xi x(i); yi y(i); k1 f(xi, yi); k2 f(xi h/2, yi h/2*k1); k3 f(xi h/2, yi h/2*k2); k4 f(xi h, yi h*k3); y(i1) yi h/6 * (k1 2*k2 2*k3 k4); end end四个 k 的权重是 1:2:2:1除以 6这是四阶精度的来源。k2、k3 都是在中点取值区别只在于用 k1 还是 k2 去估斜率这一点写反了阶数会掉到二阶误差大约变 h² 倍——h0.1 时差 100 倍左右但程序不报错只能靠和精确解比才看得出来。AB4 是四阶显式 Adams 方法用前面四个点的 f 值外推y_{i1} y_i (h/24)·[55f_i − 59f_{i-1} 37f_{i-2} − 9f_{i-3}]这里有个极易抄错的地方系数是 55、−59、37、−9前面的因子是 h/24不是 h/12。资料正文里这一行写的是 h/12标准四阶 AB4 的权系数分母是 24对 55、−59、37、−9 来说24 才能让它们加起来等于 24保证常数函数精确成立55−5937−9 24。这是个非常实用的自检四个系数加起来必须等于分母不等就说明抄错了。AB4 自身不能启动因为 i1 时需要 f_0、f_{-1}、f_{-2}前三个点得由别的单步法提供。题目要求用 RK4 提供初值那就先跑 RK4 得到 y_1、y_2、y_3从 i4 开始用 AB4 递推。这里有个坑RK4 提供的初值必须用同一个步长 h0.1不能图省事用更小的步长先算前三个点否则后面 AB4 用的 f 序列步长不一致整个递推全废。4.3 AB4-AM4 的预测校正先猜后修AB4 是显式方法用已知点外推稳定性差AM4 是隐式方法公式右边含 y_{i1}得迭代解。把两者接起来就是预测校正先用 AB4 猜一个 y_{i1}^{(p)}代进 AM4 的右端算 f_{i1}再修正一次得到 y_{i1}^{(c)}。AM4 的公式是y_{i1}^{(c)} y_i (h/24)·[9f(x_{i1}, y_{i1}^{(p)}) 19f_i − 5f_{i-1} f_{i-2}]注意这里只做一次校正不迭代到收敛。迭代到收敛当然更准但计算量上去了而且四阶方法的截断误差摆在那儿多迭代几次收益有限。只做一次校正是 PECE 模式Predict-Evaluate-Correct-Evaluate工程上最常用的折中。校正这一步把 AB4 的误差压下去大约一个量级因为隐式公式的稳定域比显式大得多。这也是为什么资料里 AB4 的误差 2 范数是 0.004809而 AB4-AM4 降到 0.000920差距五倍多。4.4 改进版多出来的那一步251/270 从哪来改进的 AB4-AM4 在预测值和校正值之间再做一次加权组合y_{i1} (251/270)·y_{i1}^{(c)} (19/270)·y_{i1}^{(p)}这两个系数不是凑的是把预测值和校正值的误差主项做抵消推出来的。251 19 270权重归一保证常值解不被破坏。校正值权重远大于预测值说明这一步是“以校正值为主、拿预测值做小幅修正”。function [x, y] ab4am4_improved(f, x0, xend, y0, h) % 改进的 AB4-AM4 预测校正前三个点由 RK4 提供 x (x0:h:xend); N length(x); y zeros(N,1); y(1) y0; % 用 RK4 起步步长必须一致 y(1:4) rk4(f, x0, x03*h, y0, h); F f(x, y); % 预先算好各点 f 值避免重复调用 for i 4:N-1 yp y(i) h/24 * (55*F(i) - 59*F(i-1) 37*F(i-2) - 9*F(i-3)); fp f(x(i1), yp); % 预测点处的斜率 yc y(i) h/24 * (9*fp 19*F(i) - 5*F(i-1) F(i-2)); y(i1) 251/270 * yc 19/270 * yp; % 加权组合 F(i1) f(x(i1), y(i1)); % 更新 f 序列供下一步用 end endF这个数组是关键优化AB4 和 AM4 都要回头取前几个点的 f 值每次都重新调f是浪费。预先存成一个向量循环里直接索引代码也短。注意F(i1)要在算完 y(i1) 之后更新顺序反了就会用到旧的 f 值。y(1:4) rk4(...)这一行返回的是 4 个点正好覆盖 i0 到 i3从 i4 开始 AB4 才有足够的四个历史点可用。这里如果返回 3 个点后面索引会从 4 开始但 F(1) 是空的MATLAB 不报错、直接拿 0 参与运算结果全错。4.5 四种方法跑同一组数据误差 2 范数排序统一步长 h0.1x 从 0 到 1.5共 16 个点。误差用和精确解之差的 2 范数xe (0:0.1:1.5); ye fx(xe); [~, y1] rk4(fxy, 0, 1.5, 3, 0.1); [~, y2] ab4(fxy, 0, 1.5, 3, 0.1); % 内部用 RK4 起步 [~, y3] ab4am4(fxy, 0, 1.5, 3, 0.1); [~, y4] ab4am4_improved(fxy, 0, 1.5, 3, 0.1); fprintf(RK4%.6f AB4%.6f AB4-AM4%.6f 改进%.6f\n, ... norm(y1-ye), norm(y2-ye), norm(y3-ye), norm(y4-ye));结果对照如下方法误差2 范数精度排名RK40.0000381AB40.0048094AB4-AM40.0009203改进的 AB4-AM40.0005312精度从高到低RK4 改进 AB4-AM4 AB4-AM4 AB4。这个排序和理论完全吻合。RK4 虽然是单步法但每一步只依赖当前点误差不累积历史信息而且四阶精度在 h0.1 这种不算小的步长下非常稳。AB4 是显式多步法用到前四个点的 f 值外推误差主项系数大55/24 量级加上稳定性区域窄误差一路累积。AB4-AM4 靠隐式校正把误差压下来五倍改进版再靠加权组合压到 RK4 的十四倍以内。值得注意的一点AB4-AM4 比 AB4 准但计算量是两倍每次迭代多算一次 f。如果只看“精度/计算量”RK4 在这道题上是完胜的——这也是为什么工程上 RK4 用得远比 Adams 系列广。Adams 系列的优势在于右端函数 f 特别贵、且步长可以取得很小时多步法每步只算一次新的 f总调用次数少于 RK4。5. 上机避坑五个反复出现的现象、原因和改法下面这五条都是上机过程中真实出现过、而且在新一批作业里还会再出现的。每一条都按“看到什么现象 → 为什么会这样 → 怎么改”写照着核对能省不少时间。5.1 公式与数据层面的坑坑一样条曲线整体平滑但节点处一阶导明显断裂。现象把 S(x) 在 [0,10] 上密集取点画图曲线看着连续但放大到节点附近能看出“折角”和理论上的 C² 连续对不上。原因三弯矩方程里 μ 和 λ 的挂靠位置写反了。μ_j 乘的是 M_{j-1}λ_j 乘的是 M_{j1}两者数值上互补和为 1写成同一行时肉眼看不出区别。方程组照样有唯一解解出来的 M 序列整体错位。解决等距情况下 μλ0.5这个特例下写反了程序也能跑对所以不能只拿等距数据测。造一组不等距的 x比如 [0 1 3 6 10]手算第一行的 μ、λ、d和程序打印的中间变量逐项对。另外在回代后加一条检查max(abs(diff(eval_spline(x,y,M,x))))应该逼近两端给定的导数值不逼近就是方程组错了。坑二三次样条第一个区间结果离谱后面全部正常。现象S(0.5) 的值大得离谱S(1.5) 之后基本合理。原因型值点录入时2.51和3.30被粘成了13.30第二个点跳高。三弯矩方程的第一行右端项直接由这个点算出来误差被放大进 M_0 和 M_1。解决数据录入后先算一阶差分检查有没有异常跳变。等距情况下差分应该缓慢变化突然出现一个 9 的跳变百分之百是录入错误。另外把 y 的均值、极差打出来看一眼车门曲线这种物理量的极差应该在合理范围内。坑三AB4 的误差大得离谱比理论值差好几个量级。现象RK4 误差 3.8e-5AB4 误差 4.8e-3差了两个数量级比预期的大。原因AB4 公式那一行的因子抄成了 h/12。分母 24 是四阶 AB4 的标配分母写成 12 相当于把整步长放大了一倍。还有一个隐藏原因前三个点如果不用同一步长的 RK4 提供f 序列的步长不一致递推公式不成立。解决用系数和检查——55 − 59 37 − 9 24等于分母。这一条对任何 Adams 类公式都适用写之前先加一遍。前三个点严格用同一个 RK4 函数、同一个 h 算不做任何“先小步长跑准一点”的优化。坑四二重积分结果比预期大一倍或小一半。现象算出来的积分值总是和参考答案差个整数倍二分次数多了也收敛不到正确答案。原因复化梯形的权重分配错了。二重积分的四个角点权重是 1/4四条边上的点是 1/2内部点是 1。常见错误是只在外层考虑权重、内层全用 1或者反过来结果正好差一个因子。解决拿 f(x,y)1 做单元测试积分值应该严格等于 (b-a)(d-c)。这个测试五秒钟能立刻暴露权重错误。然后再拿 f(x,y)xy 测积分值应该等于 (b-a)(d-c)(abcd)/2能把内层外层的变量顺序也一并验出来。5.2 程序结构与调用层面的坑坑五函数文件明明存在调用却报“未定义函数”。现象myfun.m就在当前目录主程序里myfun(0,0)报错找不到。原因MATLAB 要求文件名和函数名严格一致且函数文件不能放在脚本内部除非用局部函数语法。另一种常见情况是主程序本身是个脚本函数定义写在了脚本中间旧版本 MATLAB 会报错。解决文件名、函数名、调用名三者统一大小写也要一致。检查当前工作目录是不是文件所在目录用which myfun能直接看到 MATLAB 实际找到的是哪个文件——这个命令我每次都先用一遍比猜快得多。坑六改进的 AB4-AM4 结果和 AB4-AM4 完全一样。现象两个方法输出的误差一模一样改进似乎没生效。原因加权组合那一行写成了y(i1) yc或者系数写成了251/270 19/270 1但把 yp 和 yc 的位置写反了。位置写反的情况下结果会明显变差改成“几乎一样”通常是系数只写了一半。解决在加权那一行后面加一句if i5, fprintf(yp%.6f yc%.6f y%.6f\n, yp, yc, y(i1)); end把第一次迭代的三个数打出来手算一遍 251/270yc 19/270yp 对不对。两三行调试代码能定位的问题不要靠盯着公式找。6. 误差排序之后按精度选方法按习惯保结果三道士机题走完最有价值的结论来自第三题那张误差表RK4 0.000038、改进 AB4-AM4 0.000531、AB4-AM4 0.000920、AB4 0.004809。这个排序不是偶然它把“单步 vs 多步”“显式 vs 隐式”“是否校正”三组设计选择的代价和收益摆在了明面上。单步法误差不累积历史、显式法稳定性差、隐式校正能压误差但计算量翻倍——这三条判断依据比记住任何一个公式都管用。设计维度选择代价收益单步 / 多步单步RK4每步要算 4 次 f误差不累积起步无依赖显式 / 隐式显式AB4稳定域窄误差主项系数大每步只算 1 次 f是否校正预测校正AB4-AM4每步多算 1 次 f误差降约五倍是否加权改进版加权组合几乎为零再降约 40%真正要带走的一个技巧是用水滴法验证阶数。把步长从 h 减到 h/2看误差是不是按 2⁴ 16 倍下降。RK4 从 h0.1 换到 h0.05误差应该降到 1/16 左右如果只降到 1/4说明实际是二阶方法代码里某个地方阶数掉了。这个检查对任何数值方法都适用不依赖参考答案% 阶数验证误差比应接近 2^pp 为方法阶数 for h [0.1 0.05 0.025] xe (0:h:1.5); [~, yn] rk4(fxy, 0, 1.5, 3, h); err norm(yn - fx(xe)) / sqrt(length(xe)); % 均方根误差消除点数影响 fprintf(h%.3f rmse%.3e\n, h, err); end注意这里用的是均方根误差而不是原始 2 范数因为不同步长下点数不一样直接比 2 范数会混进点数的影响。连着三个 h 打印出来误差比稳定在 16 附近才敢说这段 RK4 写对了。回过头看我在这三道题上栽的跟头一半是公式抄写一半是数据结构没想清楚就开始敲。后来养成的习惯是任何数值程序写之前先拿一个能手工验算的小规模输入等距、两三个点、简单函数把中间变量打印出来逐项对对上了再换真题数据。这条习惯救过我不止一次——三次样条那次就是因为等距特例下 μ 和 λ 都是 0.5写反了也看不出来才拖到别人帮忙跑的时候才发现。从那以后我每次写三对角方程组都强制先用不等距的 5 个点跑一遍手算对照再也不敢拿等距数据当验收标准。希望帮到你。本文还有配套的精品资源点击获取
简介面向学习数值分析与科学计算的高校学生这份东南大学出版配套上机题文档围绕MATLAB编程实践展开聚焦三次样条插值函数、重积分计算和常微分方程初值问题数值解三类典型题目。文档不仅给出第一型三次样条插值、复化梯形公式与外推、RK4、AB4及AB4-AM4预测校正等方法的算法思路还结合车门曲线型值点、二重积分和给定初值问题展示计算过程与结果对比并附有精度排序和编程反思便于读者理解公式如何落地为程序、误差从何而来。文档同时强调主程序与子函数分离的组织方式提到逐次二分步长、外推提高精度以及公式输入错误引发的排错经历对课程作业自查和MATLAB实现复盘较有参考价值。资源包仅含1个docx文件约566KB篇幅紧凑适合上机前预习与课后对照。目前已有904人学习可作为数值分析上机训练的实用样例。1. 三道上机题卡人的从来不是公式数值分析上机题最容易给人错觉三次样条的三弯矩方程、复化梯形加外推、RK4 和 AB4-AM4 的系数公式都在书上印着抄下来似乎就完事。真坐到 MATLAB 前面翻车的往往不是推导能力而是把 μ 和 λ 的位置写反、把h/24抄成h/12、把内层外层积分变量搞混——这些地方不报错只给你一条看着挺光滑、其实整体偏了的曲线。这份数值分析上机题资料里连着做了三道题三次样条插值车门曲线、二重积分复化梯形逐次二分加外推、常微分方程初值问题RK4、AB4、AB4-AM4 和改进版对比。下面按这三道题的顺序拆把每道题的算法落点、可抄的 MATLAB 骨架、以及我实测过会踩的坑都摆出来。新手可以照着把程序跑通熟手能对一下判断依据做课程作业的、临时要补 MATLAB 数值程序的都能直接拿走。2. 三次样条插值三弯矩方程怎么补边界、怎么解三对角三次样条这道题的价值不在结果在于它把“公式推导 → 线性方程组 → 程序实现 → 结果回代”四件事完整走了一遍。第一型给定两端一阶导数的三次样条本质就是解一个三对角线性方程组求二阶导数 M再用 M 把每个小区间上的三次多项式写出来。思路清楚了代码不到五十行思路没理清改一天也调不出来。2.1 型值点先核对一遍再谈公式车门曲线的型值点一共 11 个x 从 0 到 10 等距步长 h1。资料里那行数据在复制、截图、转 Word 的过程中很容易粘连2.51和3.30挨在一起一不留神就变成13.30而后面 4.04、4.70、5.22 一路平滑上去唯独第二个点跳到 13插值出来必然是一段剧烈振荡。我一般会先把数据摆成表肉眼扫一遍量级再往程序里灌。i012345678910x_i012345678910y_i2.513.304.044.705.225.545.785.405.575.705.80边界条件是两端一阶导数已知左端 y(0)0.8右端 y(10)0.2。这两个数是第一型Clamped条件的全部输入少了任何一个方程组就欠两行M 解不出来。等距情况下 h_j 全等于 1这件事后面能用来手算校验别浪费。2.2 三弯矩方程的三个系数缺一个都解不出对内部节点 j1,…,n-1三弯矩方程是μ_j·M_{j-1} 2·M_j λ_j·M_{j1} d_j其中 h_j x_{j1} - x_jμ_j h_{j-1}/(h_{j-1}h_j)λ_j h_j/(h_{j-1}h_j) 1 - μ_j右端 d_j 6·f[x_{j-1}, x_j, x_{j1}]也就是 6 倍二阶差商。资料里把 λ 写成了 α这是符号习惯问题不影响计算但要清楚它对应的是 M_{j1} 的系数不是 M_{j-1} 的。μ 和 λ 一个挂在左邻、一个挂在右邻写反了方程组依然能解、程序依然不报错只是解出来的 M 序列会整体错位画出来的曲线在节点处一阶导不连续——肉眼不细看根本发现不了这是这道题最阴的地方。二阶差商我习惯显式写出来避免嵌套调用出错f[x_{j-1}, x_j, x_{j1}] [ (y_{j1}-y_j)/h_j − (y_j-y_{j-1})/h_{j-1} ] / (h_{j-1}h_j)等距、h1 的时候这个式子退化成 (y_{j1} - 2y_j y_{j-1})/2可以直接口算几个 d_j 来对程序输出。这是我最常用的一招先手算 j1 那一行的 μ、λ、d 三个数跟程序打印的第一行对一遍对不上就别往下跑了越跑越乱。2.3 第一型边界条件怎么补两行方程内部的 n-1 个方程不够一共 n1 个未知量 M_0 到 M_n还差两行由端点导数条件补。第一型的补法是2·M_0 M_1 (6/h_0)·( f[x_0,x_1] − y_0 )M_{n-1} 2·M_n (6/h_{n-1})·( yn − f[x{n-1},x_n] )注意两行里差商和导数的先后顺序是反的左端是“差商减导数”右端是“导数减差商”。这个顺序反了误差不会暴露成崩溃只会让整条曲线两端微微翘起。第二型给定二阶导数的补法不一样是 2M_0 M_1 6·f[x_0,x_1,x_2] 这种形式第三型周期性更不一样。三类型别在同一个程序里混用是最常见的翻车点之一做题前先确认题目给的是哪一型。拼起来就是一个 (n1)×(n1) 的三对角方程组 A·M d直接用 MATLAB 的反斜杠A\d能解但既然题目要求“编制通用程序”把追赶法Thomas 算法写出来更符合要求而且 O(n) 的复杂度在节点多的时候优势明显。2.4 通用程序的骨架构造—追赶—回代先写主函数负责组装三对角矩阵、调用追赶法、返回 Mfunction [M, h] spline3_clamped(x, y, dy0, dyn) % 第一型三次样条输入型值点 x,y 与两端一阶导数 dy0,dyn % 输出各节点二阶导数 M长度与 x 相同和区间步长 h n length(x) - 1; % 子区间个数 h diff(x); % h(j) x(j1) - x(j) A zeros(n1); d zeros(n1,1); % 左边界2*M0 M1 6/h0 * (f[x0,x1] - y0) A(1,1) 2; A(1,2) 1; d(1) 6/h(1) * ((y(2)-y(1))/h(1) - dy0); % 内部节点mu*M(j-1) 2*M(j) lam*M(j1) 6*f[x(j-1),x(j),x(j1)] for j 2:n mu h(j-1) / (h(j-1) h(j)); lam 1 - mu; A(j,j-1) mu; A(j,j) 2; A(j,j1) lam; d(j) 6 * ((y(j1)-y(j))/h(j) - (y(j)-y(j-1))/h(j-1)) ... / (h(j-1) h(j)); end % 右边界M(n-1) 2*Mn 6/h(n) * (yn - f[x(n-1),x(n)]) A(n1,n) 1; A(n1,n1) 2; d(n1) 6/h(n) * (dyn - (y(n1)-y(n))/h(n)); M thomas(A, d); % 追赶法求解见下 end追赶法把三对角矩阵压成三个向量前代消元、回代求解注意第一行和最后一行的消元边界function x thomas(A, b) % 三对角方程组追赶法A 为方阵只取三条对角线 n length(b); a diag(A,-1); bb diag(A,0); c diag(A,1); for i 2:n m a(i-1) / bb(i-1); % 消元因子 bb(i) bb(i) - m * c(i-1); b(i) b(i) - m * b(i-1); end x zeros(n,1); x(n) b(n) / bb(n); % 回代 for i n-1:-1:1 x(i) (b(i) - c(i)*x(i1)) / bb(i); end enddiag(A,-1)取的是下标 (i, i-1)对应 a(i-1)diag(A,1)取的是 (i, i1)对应 c(i)。这两个索引偏移错一位追赶法解出来的结果会差得很离谱但矩阵本身仍然可解不报错。所以写完先拿max(abs(A*x - b))验残差残差在 1e-12 量级才算通过。2.5 打印 S(i0.5) 的定位技巧题目要打印 S(0.5)、S(1.5)……S(9.5)。插值函数是分段定义的得先判断查询点落在哪个子区间。我用find(x xq, 1, last)定位左端点下标再做上下界钳位防止 xq 恰好等于最后一个节点时越界function v eval_spline(x, y, M, xq) % 用二阶导数 M 回代计算三次样条在 xq 处的值 h diff(x); v zeros(size(xq)); for k 1:numel(xq) j find(x xq(k), 1, last); j min(max(j, 1), length(x)-1); % 钳位防止越界 hj h(j); t1 (x(j1) - xq(k)) / hj; % 左基函数 t2 (xq(k) - x(j)) / hj; % 右基函数 v(k) t1*y(j) t2*y(j1) ... ((t1^3 - t1)*M(j) (t2^3 - t2)*M(j1)) * hj^2 / 6; end end这个式子是三次样条的标准回代公式(t^3 - t)那一项乘的是 M 和 h²/6不是 h/6量纲才对得上。写完之后做两个自检一是eval_spline在节点处的值必须等于 y误差 1e-12 以内二是两端的一阶导数应该逼近 0.8 和 0.2差分数值近似即可。这两条过了S(i0.5) 的结果基本不会错。3. 重积分外推从 Tn 到 Sn 的那一步才是精度来源二重积分这题看着比样条简单其实考点全在“逐次二分 外推”上。题目给的是 I(f) ∫_c^d ( ∫_a^b f(x,y) dx ) dy要求用复化梯形公式、逐次二分步长、外推算到相邻两次近似值之差的绝对值不超过 ε。如果只是把两次梯形值一减就交那外推等于没做精度停在 O(h²)二分十几次也追不上要求。3.1 为什么是 4/3 和 -1/3一次外推把误差阶抬上去复化梯形公式的误差展开是 T(h) I C·h² O(h⁴)。步长减半后 T(h/2) I C·h²/4 O(h⁴)。两式一消把 h² 项干掉S (4/3)·T(h/2) − (1/3)·T(h) I O(h⁴)这就是资料里写的 S_n (4/3)T_{2n} − (1/3)T_n。系数 4/3 和 -1/3 不是拍脑袋来的是消 h² 项消出来的4/3 − 1/3 1 保证常数项对4/3·(1/4) − 1/3 0 保证 h² 项对。写代码时这两个系数写反、或者写成 3/4 和 1/4结果会稳稳地偏离真值而且二分越多次偏离越“稳定”看起来特别像收敛了。理解这一点之后判断依据就很清楚外推后的 S 序列收敛速度应该明显快于 T 序列。如果打印出来发现 S 和 T 收敛得一样慢八成是系数或者 T_{2n}/T_n 的顺序搞反了。3.2 内层 x、外层 y二重循环的顺序不能乱题目写的是 ∫_c^d (∫_a^b f(x,y) dx) dy内层对 x 积、外层对 y 积。程序里两层循环必须严格对应外层循环变量是 y从 c 到 d内层是 x从 a 到 b。我见过不少把两层换过来的写法程序照样出数、照样收敛只是算出来的是另一个积分的值——因为积分区域是矩形的时候交换次序数学上等价但如果被积函数里的 x、y 在调用时传反了比如f(x,y)写成了f(y,x)结果就悄悄错了。被积函数 f(x,y) tan(x² y²) 里 x、y 是平方和交换了看不出差别这类题最容易掩盖参数顺序错误。我的习惯是给临时测试函数f(x,y) x*10 y跑一遍看外层积分变量是不是真的在外面两分钟就能确认比事后调半天强。积分限这一题是 a0、bπ/3、c0、dπ/6上界小、下界大这种写法在题目里也常见(b-a)一减如果是负的整个积分符号就反了。写完先把积分限代进hx (b-a)/nx打一行看看正负。3.3 三个文件的拆分myfun / getT / 主程序资料里把被积函数、复化梯形子函数、主程序分成三个文件这个拆法是对的——被积函数一改就能换题梯形子函数可以单独测。被积函数function z myfun(x, y) % 被积函数注意 x 是内层积分变量、y 是外层 z tan(x.^2 y.^2); end复化梯形二重积分注意权重角点乘 1/4边上乘 1/2内部乘 1。写代码时把权重拆成wx * wy最不容易错function T getT(f, a, b, c, d, nx, ny) % 矩形区域上的复化梯形二重积分nx/ny 为 x/y 方向等分数 hx (b-a)/nx; hy (d-c)/ny; T 0; for i 0:ny y c i*hy; wy 1; if i0 || iny, wy 0.5; end % y 方向端点权重 inner 0; for j 0:nx x a j*hx; wx 1; if j0 || jnx, wx 0.5; end % x 方向端点权重 inner inner wx * f(x, y); end T T wy * inner; end T T * hx * hy; % 最后统一乘步长 endhx*hy放在循环外面统一乘比每次内层都乘一遍快也避免累加过程中浮点误差被放大。这个写法在 nx、ny 都到几千的时候差别能看出来。3.4 二分停止条件写在哪儿决定了会不会死循环题目要求“相邻两次近似值之差的绝对值不超过 ε”ε 取 0.5e-5。关键是这个“近似值”指的是外推后的 S不是原始的 T。判据写成abs(T2 - T1) eps就偏了题意写成abs(S2 - S1) eps才对。eps0 0.5e-5; nx 2; ny 2; Tn getT(myfun, 0, pi/3, 0, pi/6, nx, ny); % 粗网格 nx 2*nx; ny 2*ny; T2n getT(myfun, 0, pi/3, 0, pi/6, nx, ny); Sn 4/3*T2n - 1/3*Tn; % 第一次外推 it 0; while it 20 % 上限兜底防死循环 Tn T2n; nx 2*nx; ny 2*ny; T2n getT(myfun, 0, pi/3, 0, pi/6, nx, ny); S2n 4/3*T2n - 1/3*Tn; if abs(S2n - Sn) eps0 break; end Sn S2n; it it 1; end fprintf(积分值 %.10f二分次数 %d\n, S2n, it);it 20这个上限不是可选项是必需品。浮点数在 1e-5 附近会进入震荡区两个相邻的 S 值可能来回跳永远不满足小于 ε没上限就是死循环。而且二分到 nx、ny 上千以后单次 getT 的耗时是平方级增长的卡死一次能等十几分钟。3.5 运行慢的真因每次二分都在重算旧节点资料里反思“运行速度慢是很大的不足”这个判断是准确的但原因不是“MATLAB 慢”而是每次getT都从头把所有节点重算了一遍。实际上复化梯形公式有递推关系T_{2n} T_n / 2 (h/2) · Σ f(新增节点)新增节点的横坐标是 a (2k-1)h/2数量正好是 n 个。用这个递推每次二分的计算量从 O(4n²) 降到 O(n²) 的一小部分尤其两个方向同时二分的时候节省非常明显。改法是把getT拆成“首算”和“加细”两个函数主循环里只算新增点。做作业不一定要求但把这段写进反思里比写“还需要更深入学习 MATLAB”有说服力得多。4. 常微分方程初值RK4、AB4、AB4-AM4 和改进版的误差差在哪第三道题最有意思因为它把四种方法放在同一个初值问题上横向比。初值问题是 y −x²y²0 ≤ x ≤ 1.5y(0) 3步长 h 0.1精确解 y 3/(1x³)。四种方法、同一组数据、同一个误差范数排序出来一目了然。这道题真正考的不是“会不会写 RK4”而是“知不知道为什么 AB4 会比 RK4 差两个数量级”。4.1 先把精确解和右端函数准备好右端函数和精确解各写一个文件后面所有方法共用function f fxy(x, y) % 一阶导数右端项 f -x.^2 .* y.^2; end function y fx(x) % 精确解用于误差对比 y 3 ./ (1 x.^3); end精确解这个文件很容易被当成“可有可无”实际它是整道题的裁判。没有它四种方法的误差没法排序题目第6问的“能得到哪些结论”就只能凭感觉写。我一般会把精确解和数值解画在同一张图上再加一张误差随 x 变化的图误差曲线的形状比单个数字能说明问题——RK4 的误差曲线几乎是平的AB4 的误差曲线一路往上爬这个差别一眼就能看出来。验证一下精确解和右端项一致对 y 3/(1x³) 求导得 −9x²/(1x³)²而 −x²y² −9x²/(1x³)²两边相等。这一步花三十秒能避免后面所有对比都建立在一个写错的 fx 上。4.2 RK4 是单步自启动AB4 不是RK4 是单步法只依赖 y_i 就能算出 y_{i1}所以第一个点起步毫无障碍function [x, y] rk4(f, x0, xend, y0, h) % 经典四阶 Runge-Kutta x (x0:h:xend); N length(x); y zeros(N,1); y(1) y0; for i 1:N-1 xi x(i); yi y(i); k1 f(xi, yi); k2 f(xi h/2, yi h/2*k1); k3 f(xi h/2, yi h/2*k2); k4 f(xi h, yi h*k3); y(i1) yi h/6 * (k1 2*k2 2*k3 k4); end end四个 k 的权重是 1:2:2:1除以 6这是四阶精度的来源。k2、k3 都是在中点取值区别只在于用 k1 还是 k2 去估斜率这一点写反了阶数会掉到二阶误差大约变 h² 倍——h0.1 时差 100 倍左右但程序不报错只能靠和精确解比才看得出来。AB4 是四阶显式 Adams 方法用前面四个点的 f 值外推y_{i1} y_i (h/24)·[55f_i − 59f_{i-1} 37f_{i-2} − 9f_{i-3}]这里有个极易抄错的地方系数是 55、−59、37、−9前面的因子是 h/24不是 h/12。资料正文里这一行写的是 h/12标准四阶 AB4 的权系数分母是 24对 55、−59、37、−9 来说24 才能让它们加起来等于 24保证常数函数精确成立55−5937−9 24。这是个非常实用的自检四个系数加起来必须等于分母不等就说明抄错了。AB4 自身不能启动因为 i1 时需要 f_0、f_{-1}、f_{-2}前三个点得由别的单步法提供。题目要求用 RK4 提供初值那就先跑 RK4 得到 y_1、y_2、y_3从 i4 开始用 AB4 递推。这里有个坑RK4 提供的初值必须用同一个步长 h0.1不能图省事用更小的步长先算前三个点否则后面 AB4 用的 f 序列步长不一致整个递推全废。4.3 AB4-AM4 的预测校正先猜后修AB4 是显式方法用已知点外推稳定性差AM4 是隐式方法公式右边含 y_{i1}得迭代解。把两者接起来就是预测校正先用 AB4 猜一个 y_{i1}^{(p)}代进 AM4 的右端算 f_{i1}再修正一次得到 y_{i1}^{(c)}。AM4 的公式是y_{i1}^{(c)} y_i (h/24)·[9f(x_{i1}, y_{i1}^{(p)}) 19f_i − 5f_{i-1} f_{i-2}]注意这里只做一次校正不迭代到收敛。迭代到收敛当然更准但计算量上去了而且四阶方法的截断误差摆在那儿多迭代几次收益有限。只做一次校正是 PECE 模式Predict-Evaluate-Correct-Evaluate工程上最常用的折中。校正这一步把 AB4 的误差压下去大约一个量级因为隐式公式的稳定域比显式大得多。这也是为什么资料里 AB4 的误差 2 范数是 0.004809而 AB4-AM4 降到 0.000920差距五倍多。4.4 改进版多出来的那一步251/270 从哪来改进的 AB4-AM4 在预测值和校正值之间再做一次加权组合y_{i1} (251/270)·y_{i1}^{(c)} (19/270)·y_{i1}^{(p)}这两个系数不是凑的是把预测值和校正值的误差主项做抵消推出来的。251 19 270权重归一保证常值解不被破坏。校正值权重远大于预测值说明这一步是“以校正值为主、拿预测值做小幅修正”。function [x, y] ab4am4_improved(f, x0, xend, y0, h) % 改进的 AB4-AM4 预测校正前三个点由 RK4 提供 x (x0:h:xend); N length(x); y zeros(N,1); y(1) y0; % 用 RK4 起步步长必须一致 y(1:4) rk4(f, x0, x03*h, y0, h); F f(x, y); % 预先算好各点 f 值避免重复调用 for i 4:N-1 yp y(i) h/24 * (55*F(i) - 59*F(i-1) 37*F(i-2) - 9*F(i-3)); fp f(x(i1), yp); % 预测点处的斜率 yc y(i) h/24 * (9*fp 19*F(i) - 5*F(i-1) F(i-2)); y(i1) 251/270 * yc 19/270 * yp; % 加权组合 F(i1) f(x(i1), y(i1)); % 更新 f 序列供下一步用 end endF这个数组是关键优化AB4 和 AM4 都要回头取前几个点的 f 值每次都重新调f是浪费。预先存成一个向量循环里直接索引代码也短。注意F(i1)要在算完 y(i1) 之后更新顺序反了就会用到旧的 f 值。y(1:4) rk4(...)这一行返回的是 4 个点正好覆盖 i0 到 i3从 i4 开始 AB4 才有足够的四个历史点可用。这里如果返回 3 个点后面索引会从 4 开始但 F(1) 是空的MATLAB 不报错、直接拿 0 参与运算结果全错。4.5 四种方法跑同一组数据误差 2 范数排序统一步长 h0.1x 从 0 到 1.5共 16 个点。误差用和精确解之差的 2 范数xe (0:0.1:1.5); ye fx(xe); [~, y1] rk4(fxy, 0, 1.5, 3, 0.1); [~, y2] ab4(fxy, 0, 1.5, 3, 0.1); % 内部用 RK4 起步 [~, y3] ab4am4(fxy, 0, 1.5, 3, 0.1); [~, y4] ab4am4_improved(fxy, 0, 1.5, 3, 0.1); fprintf(RK4%.6f AB4%.6f AB4-AM4%.6f 改进%.6f\n, ... norm(y1-ye), norm(y2-ye), norm(y3-ye), norm(y4-ye));结果对照如下方法误差2 范数精度排名RK40.0000381AB40.0048094AB4-AM40.0009203改进的 AB4-AM40.0005312精度从高到低RK4 改进 AB4-AM4 AB4-AM4 AB4。这个排序和理论完全吻合。RK4 虽然是单步法但每一步只依赖当前点误差不累积历史信息而且四阶精度在 h0.1 这种不算小的步长下非常稳。AB4 是显式多步法用到前四个点的 f 值外推误差主项系数大55/24 量级加上稳定性区域窄误差一路累积。AB4-AM4 靠隐式校正把误差压下来五倍改进版再靠加权组合压到 RK4 的十四倍以内。值得注意的一点AB4-AM4 比 AB4 准但计算量是两倍每次迭代多算一次 f。如果只看“精度/计算量”RK4 在这道题上是完胜的——这也是为什么工程上 RK4 用得远比 Adams 系列广。Adams 系列的优势在于右端函数 f 特别贵、且步长可以取得很小时多步法每步只算一次新的 f总调用次数少于 RK4。5. 上机避坑五个反复出现的现象、原因和改法下面这五条都是上机过程中真实出现过、而且在新一批作业里还会再出现的。每一条都按“看到什么现象 → 为什么会这样 → 怎么改”写照着核对能省不少时间。5.1 公式与数据层面的坑坑一样条曲线整体平滑但节点处一阶导明显断裂。现象把 S(x) 在 [0,10] 上密集取点画图曲线看着连续但放大到节点附近能看出“折角”和理论上的 C² 连续对不上。原因三弯矩方程里 μ 和 λ 的挂靠位置写反了。μ_j 乘的是 M_{j-1}λ_j 乘的是 M_{j1}两者数值上互补和为 1写成同一行时肉眼看不出区别。方程组照样有唯一解解出来的 M 序列整体错位。解决等距情况下 μλ0.5这个特例下写反了程序也能跑对所以不能只拿等距数据测。造一组不等距的 x比如 [0 1 3 6 10]手算第一行的 μ、λ、d和程序打印的中间变量逐项对。另外在回代后加一条检查max(abs(diff(eval_spline(x,y,M,x))))应该逼近两端给定的导数值不逼近就是方程组错了。坑二三次样条第一个区间结果离谱后面全部正常。现象S(0.5) 的值大得离谱S(1.5) 之后基本合理。原因型值点录入时2.51和3.30被粘成了13.30第二个点跳高。三弯矩方程的第一行右端项直接由这个点算出来误差被放大进 M_0 和 M_1。解决数据录入后先算一阶差分检查有没有异常跳变。等距情况下差分应该缓慢变化突然出现一个 9 的跳变百分之百是录入错误。另外把 y 的均值、极差打出来看一眼车门曲线这种物理量的极差应该在合理范围内。坑三AB4 的误差大得离谱比理论值差好几个量级。现象RK4 误差 3.8e-5AB4 误差 4.8e-3差了两个数量级比预期的大。原因AB4 公式那一行的因子抄成了 h/12。分母 24 是四阶 AB4 的标配分母写成 12 相当于把整步长放大了一倍。还有一个隐藏原因前三个点如果不用同一步长的 RK4 提供f 序列的步长不一致递推公式不成立。解决用系数和检查——55 − 59 37 − 9 24等于分母。这一条对任何 Adams 类公式都适用写之前先加一遍。前三个点严格用同一个 RK4 函数、同一个 h 算不做任何“先小步长跑准一点”的优化。坑四二重积分结果比预期大一倍或小一半。现象算出来的积分值总是和参考答案差个整数倍二分次数多了也收敛不到正确答案。原因复化梯形的权重分配错了。二重积分的四个角点权重是 1/4四条边上的点是 1/2内部点是 1。常见错误是只在外层考虑权重、内层全用 1或者反过来结果正好差一个因子。解决拿 f(x,y)1 做单元测试积分值应该严格等于 (b-a)(d-c)。这个测试五秒钟能立刻暴露权重错误。然后再拿 f(x,y)xy 测积分值应该等于 (b-a)(d-c)(abcd)/2能把内层外层的变量顺序也一并验出来。5.2 程序结构与调用层面的坑坑五函数文件明明存在调用却报“未定义函数”。现象myfun.m就在当前目录主程序里myfun(0,0)报错找不到。原因MATLAB 要求文件名和函数名严格一致且函数文件不能放在脚本内部除非用局部函数语法。另一种常见情况是主程序本身是个脚本函数定义写在了脚本中间旧版本 MATLAB 会报错。解决文件名、函数名、调用名三者统一大小写也要一致。检查当前工作目录是不是文件所在目录用which myfun能直接看到 MATLAB 实际找到的是哪个文件——这个命令我每次都先用一遍比猜快得多。坑六改进的 AB4-AM4 结果和 AB4-AM4 完全一样。现象两个方法输出的误差一模一样改进似乎没生效。原因加权组合那一行写成了y(i1) yc或者系数写成了251/270 19/270 1但把 yp 和 yc 的位置写反了。位置写反的情况下结果会明显变差改成“几乎一样”通常是系数只写了一半。解决在加权那一行后面加一句if i5, fprintf(yp%.6f yc%.6f y%.6f\n, yp, yc, y(i1)); end把第一次迭代的三个数打出来手算一遍 251/270yc 19/270yp 对不对。两三行调试代码能定位的问题不要靠盯着公式找。6. 误差排序之后按精度选方法按习惯保结果三道士机题走完最有价值的结论来自第三题那张误差表RK4 0.000038、改进 AB4-AM4 0.000531、AB4-AM4 0.000920、AB4 0.004809。这个排序不是偶然它把“单步 vs 多步”“显式 vs 隐式”“是否校正”三组设计选择的代价和收益摆在了明面上。单步法误差不累积历史、显式法稳定性差、隐式校正能压误差但计算量翻倍——这三条判断依据比记住任何一个公式都管用。设计维度选择代价收益单步 / 多步单步RK4每步要算 4 次 f误差不累积起步无依赖显式 / 隐式显式AB4稳定域窄误差主项系数大每步只算 1 次 f是否校正预测校正AB4-AM4每步多算 1 次 f误差降约五倍是否加权改进版加权组合几乎为零再降约 40%真正要带走的一个技巧是用水滴法验证阶数。把步长从 h 减到 h/2看误差是不是按 2⁴ 16 倍下降。RK4 从 h0.1 换到 h0.05误差应该降到 1/16 左右如果只降到 1/4说明实际是二阶方法代码里某个地方阶数掉了。这个检查对任何数值方法都适用不依赖参考答案% 阶数验证误差比应接近 2^pp 为方法阶数 for h [0.1 0.05 0.025] xe (0:h:1.5); [~, yn] rk4(fxy, 0, 1.5, 3, h); err norm(yn - fx(xe)) / sqrt(length(xe)); % 均方根误差消除点数影响 fprintf(h%.3f rmse%.3e\n, h, err); end注意这里用的是均方根误差而不是原始 2 范数因为不同步长下点数不一样直接比 2 范数会混进点数的影响。连着三个 h 打印出来误差比稳定在 16 附近才敢说这段 RK4 写对了。回过头看我在这三道题上栽的跟头一半是公式抄写一半是数据结构没想清楚就开始敲。后来养成的习惯是任何数值程序写之前先拿一个能手工验算的小规模输入等距、两三个点、简单函数把中间变量打印出来逐项对对上了再换真题数据。这条习惯救过我不止一次——三次样条那次就是因为等距特例下 μ 和 λ 都是 0.5写反了也看不出来才拖到别人帮忙跑的时候才发现。从那以后我每次写三对角方程组都强制先用不等距的 5 个点跑一遍手算对照再也不敢拿等距数据当验收标准。希望帮到你。本文还有配套的精品资源点击获取
