★ FEATURED ARTICLE 简介一份围绕数值分析课程上机作业整理的 docx 文档面向正在学习数值分析、需要完成 MATLAB 编程实验的本科生与研究生。内容以三次样条插值函数、重积分计算和常微分方程初值问题数值解三条主线展开既给出第一型三次样条通用程序的构造思路与车门曲线算例也梳理复化梯形公式配合逐次二分与外推求重积分的实现方式并对比 RK4、AB4、AB4-AM4 预测校正及改进格式在初值问题中的计算精度。压缩包共 1 个 docx 文件约 566KB篇幅紧凑适合作为课程作业对照与算法复习材料。已有 905 人学习读者可从中获得公式推导、程序组织方式、结果验证与误差排序等参考便于理解东南大学版教材相关上机题的求解流程并据此搭建自己的 MATLAB 实验脚本。1. 从一份上机作业说起三道题串起数值分析最容易翻车的三条链路数值分析的上机题文档看着像三份互不相关的作业实际是三条会在真实计算里反复出现的链路插值、积分、微分方程初值问题。这份资料用 MATLAB 把三道题做完了——三次样条插值求车门曲线、复化梯形加外推算二重积分、RK4 与 AB4-AM4 解常微分方程——每个题都给公式推导、程序结构和运行结果截图末尾还附了误差对比和方法排序。适合两类人正在啃数值分析课程、需要照着把作业跑出正确数字的同学以及手上有一堆 MATLAB 老脚本、想看清楚「主程序 子函数」这套写法怎么组织的人。它最大的价值不是答案本身而是把公式怎么落到索引、边界条件怎么塞进三对角方程组、精度不够时该怀疑哪一行代码这些没人写在教材里的事摆了出来。2. 运行链路拆解exp5.m 主脚本、getT.m 子函数与重积分的外推循环2.1 目录与文件职责划分这份资料的目录组织没有统一工程模板是典型的课程作业式布局按章分文件夹每道题一个主脚本通用逻辑抽成子函数。我按原文提到的 M 文件梳理了一下合理的结构大致是这样路径文件职责ch4_spline/exp4.m三次样条主程序读入型值点、调用求解、打印 S(i0.5)ch4_spline/spline3.m第一型三次样条通用函数组装并解三对角方程组ch4_spline/thomas.m追赶法解三对角方程组也可内联ch5_integral/exp5.m重积分主程序控制二分与外推循环ch5_integral/getT.m双层复化梯形返回给定步长下的积分近似ch5_integral/myfun.m被积函数单独成文件ch6_ode/exp6.m常微分方程主程序四种方法依次计算并对比ch6_ode/fxy.m右端函数 f(x,y)ch6_ode/fx.m精确解用于误差计算这种拆法的好处在改题的时候特别明显换一个被积函数只动 myfun.m换一个右端函数只动 fxy.m求解框架一行不用碰。原文里写「对于不同的定积分的形式也写作单独的子函数保证其修改的灵活性」说的就是这个意思。反过来如果所有东西都塞进一个脚本第二次做同类题基本只能复制粘贴再改改错了还找不出来。提示MATLAB 要求函数文件名与函数名一致myfun.m 里第一行必须是 function z myfun(x,y)否则调用会报「未定义函数或变量」。这一点新手最容易忽略尤其是从别人那里拷来文件后随手重命名。2.2 主脚本怎么调子函数以重积分为例重积分这道题的要求是给定 I(f)∫_c^d (∫_a^b f(x,y)dx)dy取初始步长 h、k 和精度 ε用复化梯形公式做逐次二分并应用外推思想算到相邻两次近似值之差不超过 ε。核心思路是先用步长 (h,k) 算一个 T_n再把两个方向步长都减半算 T_2n用线性组合得到 Simpson 级别的近似S_n (4/3)·T_2n − (1/3)·T_n主脚本负责的是循环和判据数值计算全部下沉到 getT.m% exp5.m 重积分主程序二分 外推 a 0; b pi/3; % 内层 x 的积分区间 c 0; d pi/6; % 外层 y 的积分区间 eps 0.5e-5; % 题面精度要求 maxit 20; % 二分次数上限防止死循环 h (b-a)/2; k (d-c)/2; % 初始步长先粗分两段 T1 getT(a,b,c,d,h,k); % 粗网格的复化梯形值 Sold T1; % 保存上一次外推结果 for it 1:maxit h h/2; k k/2; % 两个方向同时二分 T2 getT(a,b,c,d,h,k); % 细网格的复化梯形值 Snew 4/3*T2 - 1/3*T1; % 外推组合得到 Simpson 级精度 fprintf(iter%2d T2n%.10f S%.10f diff%.2e\n, ... it, T2, Snew, abs(Snew-Sold)); if abs(Snew-Sold) eps % 相邻两次近似值之差满足精度 break; end Sold Snew; T1 T2; % 滚动前进注意 T1 要更新成 T2 end fprintf(最终结果 I %.10f二分次数 %d\n, Snew, it);这段代码里有两个容易写错的地方。一是 T1 的滚动更新外推需要的是「二分前后两个梯形值」所以每轮结束必须把 T2 赋给 T1漏了这一句第二轮开始外推用的还是最初那个粗网格值结果会稳定地偏。二是 maxit 上限题面只说「算到满足精度为止」但没说不满足怎么办加上限之后至少能看到每轮的 diff 是在收敛还是在震荡。getT.m 是双层复化梯形外层对 y、内层对 x各自端点权重减半function T getT(a,b,c,d,h,k) % 双层复化梯形外层 y 步长 k内层 x 步长 h ny round((d-c)/k); % 外层分段数 T 0; for iy 0:ny y c iy*k; if iy 0 || iy ny wy 0.5; % 复化梯形端点权重 1/2 else wy 1; end T T wy * innerTrap(a,b,h,y); end T T * k; % 外层乘步长 k end function s innerTrap(a,b,h,y) % 固定 y 之后对内层 x 做复化梯形 nx round((b-a)/h); s 0; for ix 0:nx x a ix*h; if ix 0 || ix nx wx 0.5; % 内层端点同样减半 else wx 1; end s s wx * myfun(x,y); end s s * h; % 内层乘步长 h end参数含义直白a、b 是内层积分上下限c、d 是外层h、k 是两个方向的当前步长。分段数用 round((b-a)/h) 而不是直接 (b-a)/h是为了避免浮点除法带来的 3.9999999 被当成 3 段这个细节在步长反复二分之后特别重要血泪经验。2.3 被积函数独立成文件与运行方式myfun.m 只有一行有效代码function z myfun(x,y) % 被积函数tan(x^2 y^2) z tan(x.^2 y.^2); end写成独立函数文件而不是在主脚本末尾用匿名函数 (x,y) tan(x.^2y.^2)主要是为了调试时能单独测在命令行敲 myfun(0.1,0.2)看返回值是否合理比在整个二重循环里打断点快得多。跑起来的步骤不复杂但顺序要对打开 MATLAB把当前工作目录切到 ch5_integral 文件夹或者在「设置路径」里把该文件夹加进去命令行输入 exp5回车观察每轮打印的 diff 是否单调下降最后一行会给出最终积分值和二分次数想改精度就改 exp5.m 里的 eps想换被积函数就改 myfun.m主程序不用动。原文档里的运行结果截图显示积分值收敛到小数后若干位稳定具体数值随 eps 和初始步长略有差异。核对结果时可以拿一个能解析积分的函数替换 myfun.m 试跑比如 f x.*y它在一象限矩形上的二重积分有闭式解 (b²−a²)(d²−c²)/4用这个基准去验证你的外推循环写没写对比盯着 tan 的结果瞎猜靠谱得多。3. 三次样条插值μ/α/d 三对角组装、追赶法求 M 与车门曲线 S(i0.5)3.1 从公式到三对角方程组第一型三次样条的意思是给定两个端点的一阶导数值 y₀′ 和 yₙ′其余节点只给函数值。设节点 x₀ x₁ … xₙ步长 hⱼ xⱼ₊₁ − xⱼ二阶导数弯矩Mⱼ 满足内点μⱼMⱼ₋₁ 2Mⱼ αⱼMⱼ₊₁ dⱼj 1…n−1其中 μⱼ hⱼ₋₁/(hⱼ₋₁hⱼ)αⱼ 1 − μⱼ hⱼ/(hⱼ₋₁hⱼ)dⱼ 6·f[xⱼ₋₁,xⱼ,xⱼ₊₁] 6·((yⱼ₊₁−yⱼ)/hⱼ − (yⱼ−yⱼ₋₁)/hⱼ₋₁)/(hⱼ₋₁hⱼ)两个边界条件补两行左端2M₀ M₁ (6/h₀)·((y₁−y₀)/h₀ − y₀′)右端Mₙ₋₁ 2Mₙ (6/hₙ₋₁)·(yₙ′ − (yₙ−yₙ₋₁)/hₙ₋₁)一共 n1 个未知数、n1 个方程系数矩阵是三对角且严格对角占优用追赶法解最省事。原文档里提到「公式中的 lamda 没有找到对应符号用 α 代替」其实就是这个 αⱼ数学书写里常记作 λⱼ代码里叫什么无所谓别和 μ 弄混就行。把上面的推导落成 MATLAB第一段是组装function M spline3(x, y, dy0, dyn) % 第一型三次样条返回各节点的二阶导数 M n length(x) - 1; % 区间数节点数为 n1 h diff(x); % h(j) x(j1) - x(j)长度 n % 差商f[x(j-1),x(j)] 与 f[x(j),x(j1)] f1 diff(y) ./ h; % 一阶差商长度 n % 组装三对角a 下对角、b 主对角、c 上对角、r 右端 a zeros(n1,1); b zeros(n1,1); c zeros(n1,1); r zeros(n1,1); % 左边界2*M0 M1 rhs b(1) 2; c(1) 1; r(1) 6/h(1) * (f1(1) - dy0); % 内点 j 1..n-1MATLAB 下标 2..n for j 2:n mu h(j-1) / (h(j-1) h(j)); alp h(j) / (h(j-1) h(j)); d 6 * (f1(j) - f1(j-1)) / (h(j-1) h(j)); a(j) mu; b(j) 2; c(j) alp; r(j) d; end % 右边界M(n-1) 2*Mn rhs a(n1) 1; b(n1) 2; r(n1) 6/h(n) * (dyn - f1(n)); M thomas(a, b, c, r); % 追赶法求解 end逻辑上有两处必须对齐下标h 的长度是 n区间数f1 长度也是 n而 M 的长度是 n1节点数。内点循环跑 j 2…n正好覆盖 M₁ 到 Mₙ₋₁边界两行分别占 M₀ 和 Mₙ 的位置。参数 dy0、dyn 就是题面的 y₀′0.8、y₁₀′0.2。3.2 追赶法不调 \ 也能解MATLAB 里直接写 M A\b 是最省事的但作业的意图是让你把追赶法写出来。三对角方程组的追赶法就是先消元再回代O(n) 复杂度function x thomas(a, b, c, r) % 追赶法解三对角方程组a 下对角、b 主对角、c 上对角、r 右端 n length(r); cp zeros(n,1); dp zeros(n,1); cp(1) c(1) / b(1); dp(1) r(1) / b(1); % 追前向消元 for i 2:n m b(i) - a(i) * cp(i-1); % 主元 cp(i) c(i) / m; dp(i) (r(i) - a(i) * dp(i-1)) / m; end % 赶回代 x zeros(n,1); x(n) dp(n); for i n-1:-1:1 x(i) dp(i) - cp(i) * x(i1); end end需要注意 a(1) 和 c(n) 虽然用不到但数组必须存在传零进去更安全。这段代码和直接 A\b 的结果应当在 1e-12 量级内一致写完先拿一个 3×3 的小系统手算对照通过之后再上 11 个节点的车门数据。3.3 回代求 S(x) 与打印 S(i0.5)拿到 M 之后每个小区间 [xⱼ, xⱼ₊₁] 上的三次多项式写成S(x) Mⱼ(xⱼ₊₁−x)³/(6hⱼ) Mⱼ₊₁(x−xⱼ)³/(6hⱼ) (yⱼ − Mⱼhⱼ²/6)(xⱼ₊₁−x)/hⱼ (yⱼ₊₁ − Mⱼ₊₁hⱼ²/6)(x−xⱼ)/hⱼ对应的求值函数用一个循环定位区间再套公式function s evalS(x, y, M, xq) % 在查询点 xq 处计算三次样条 S(xq) n length(x) - 1; s zeros(size(xq)); for t 1:numel(xq) xv xq(t); j find(x xv, 1, last); % 定位所在区间 if j n, j n; end % 右端点保护 h x(j1) - x(j); A (x(j1) - xv)^3 / (6*h); B (xv - x(j))^3 / (6*h); C (y(j) - M(j) *h^2/6) * (x(j1) - xv) / h; D (y(j1) - M(j1)*h^2/6) * (xv - x(j)) / h; s(t) M(j)*A M(j1)*B C D; end end主程序把题面数据填进去按 0.5 的间隔求值并打印x 0:10; y [2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80]; dy0 0.8; dyn 0.2; % 端点导数条件 M spline3(x, y, dy0, dyn); % 求弯矩 xq 0.5:1:9.5; % i0.5, i0..9 sq evalS(x, y, M, xq); for i 1:numel(xq) fprintf(S(%.1f) %.6f\n, xq(i), sq(i)); end跑完先做两个自检把 xq 换成原节点 x算出的 S(xⱼ) 应当逐个等于 yⱼ误差在 1e-10 以内把端点导数换成有限差分近似结果应当和现在相差不大。两条都过了再去看 S(i0.5) 的数值那才是题面要交的输出。4. 常微分方程初值问题RK4 起步、AB4-AM4 预测校正与四种格式的误差排序4.1 右端函数与精确解各占一个文件题面给的初值问题是 y′ −x²y²0 ≤ x ≤ 1.5y(0) 3步长 h 0.1精确解 y(x) 3/(1x³)。把右端函数和精确解拆成两个文件是为了让四种方法共用同一套输入function z fxy(x, y) % 右端函数 f(x,y) -x^2 * y^2 z -x.^2 .* y.^2; end function z fx(x) % 精确解用于误差统计 z 3 ./ (1 x.^3); end这看起来是小事但四种方法各写一遍 −x²y²改题的时候就要改四个地方漏一个就会出现「三种方法对、一种方法离谱」的假象排查半天发现是抄错了。4.2 RK4四阶单步方法经典四级四阶 Runge-Kuttayᵢ₊₁ yᵢ h/6·(k₁ 2k₂ 2k₃ k₄)k₁ f(xᵢ, yᵢ)k₂ f(xᵢ h/2, yᵢ h/2·k₁)k₃ f(xᵢ h/2, yᵢ h/2·k₂)k₄ f(xᵢ h, yᵢ h·k₃)function [x, y] rk4(f, a, b, y0, h) % 四阶 Runge-Kutta x (a:h:b); N length(x); y zeros(N,1); y(1) y0; for i 1:N-1 k1 f(x(i), y(i)); k2 f(x(i)h/2, y(i)h/2*k1); k3 f(x(i)h/2, y(i)h/2*k2); k4 f(x(i)h, y(i)h*k3); y(i1) y(i) h/6*(k1 2*k2 2*k3 k4); end end单步方法的好处是自启动只要给 y(0) 就能一路算到底不需要额外提供前几个点的值。AB4 系列恰恰相反四阶格式要用到前四步的 f 值必须先用 RK4 把前 4 个点算出来。4.3 AB4 显式与 AM4 隐式的预测校正四阶 Adams-Bashforth 显式格式yᵢ₊₁ yᵢ h/24·(55fᵢ − 59fᵢ₋₁ 37fᵢ₋₂ − 9fᵢ₋₃)四阶 Adams-Moulton 隐式格式yᵢ₊₁ yᵢ h/24·(9fᵢ₊₁ 19fᵢ − 5fᵢ₋₁ fᵢ₋₂)隐式里带 fᵢ₊₁ f(xᵢ₊₁, yᵢ₊₁)未知数在右边所以先用显式预测一个 ŷᵢ₊₁再代入隐式做一次校正这就是预测-校正法PECEfunction [x, y] ab4am4(f, a, b, y0, h, mode) % mode: ab4 纯显式 | pc 预测校正 | ipc 改进的预测校正 x (a:h:b); N length(x); y zeros(N,1); y(1) y0; % 用 RK4 提供前 4 个点 for i 1:3 k1 f(x(i), y(i)); k2 f(x(i)h/2, y(i)h/2*k1); k3 f(x(i)h/2, y(i)h/2*k2); k4 f(x(i)h, y(i)h*k3); y(i1) y(i) h/6*(k1 2*k2 2*k3 k4); end for i 4:N-1 f0 f(x(i), y(i)); f1 f(x(i-1), y(i-1)); f2 f(x(i-2), y(i-2)); f3 f(x(i-3), y(i-3)); % 预测AB4 显式 yp y(i) h/24*(55*f0 - 59*f1 37*f2 - 9*f3); switch mode case ab4 y(i1) yp; case pc % 校正AM4 隐式用 yp 代 f(i1) yc y(i) h/24*(9*f(x(i1), yp) 19*f0 - 5*f1 f2); y(i1) yc; case ipc % 改进预测值与校正值加权 yc y(i) h/24*(9*f(x(i1), yp) 19*f0 - 5*f1 f2); y(i1) 251/270*yc 19/270*yp; end end end三处细节值得拎出来说。第一h/24 这个系数四阶 Adams 公式的标准形式分母就是 24写成 h/12 再配 55 的系数会让步长整体偏大一倍程序不会报错只是结果差了量级这种坑在避坑章里还会再提。第二前 4 个点用 RK4 算出来之后循环从 i 4 开始MATLAB 从 1 开始计数对应数学上的第 4 步此时 f₀、f₁、f₂、f₃ 都拿得到。第三改进格式里的权重 251/270 和 19/270 是题面给的校正式权重约 0.93预测式约 0.07本质上是把校正值再往预测值方向拉一点抑制预测校正系统的振荡。4.4 误差统计与方法排序主程序把四种方法都跑一遍跟精确解对比用 2-范数衡量整体误差h 0.1; a 0; b 1.5; y0 3; [x, ye] odeEulerGrid(a, b, h); % 生成网格并算精确解 ye fx(x); [~, y1] rk4 ( fxy, a, b, y0, h); [~, y2] ab4am4(fxy, a, b, y0, h, ab4); [~, y3] ab4am4(fxy, a, b, y0, h, pc); [~, y4] ab4am4(fxy, a, b, y0, h, ipc); e (y) norm(y - ye); % 2-范数误差 fprintf(RK4 err%.6e\n, e(y1)); fprintf(AB4 err%.6e\n, e(y2)); fprintf(AB4-AM4 err%.6e\n, e(y3)); fprintf(改进型 err%.6e\n, e(y4));原文档给出的误差分别是 3.8e-5、4.809e-3、9.2e-4、5.31e-4排序是 RK4 改进型 AB4-AM4 AB4-AM4 AB4。这个排序很符合理论预期RK4 是四阶单步每步都重新评估斜率稳定性最好AB4 是显式多步误差累积快加一步校正能把误差压下去一个量级再对预测值和校正值加权还能再降一点。把这张表和误差数字一起看比只背「RK4 精度高」有用得多。5. 常见问题排查五个不报错但结果全错的坑5.1 坑一Adams 公式系数分母写成 12现象AB4 和 AB4-AM4 的结果都偏大跟精确解的趋势一致但数值系统性偏高RK4 却完全正常。原因四阶 Adams-Bashforth 的标准形式是 yᵢ₊₁ yᵢ h/24·(55fᵢ − 59fᵢ₋₁ 37fᵢ₋₂ − 9fᵢ₋₃)。有些推导把括号内系数整体乘 2 写成 h/12·(110fᵢ − 118fᵢ₋₁ 74fᵢ₋₂ − 18fᵢ₋₃)也等价。但只改了分母、没改括号内系数等价于把步长放大了一倍误差自然上升一个量级而且程序不会报任何错。解决写完公式先拿一个能解析求解的简单方程对一遍比如 y′ y、y(0) 1精确解是 eˣ跑十步看误差量级是否落在 O(h⁴) 附近。对不上就老老实实回教材核对每一项系数别凭记忆。5.2 坑二样条回代时 M 的下标错位现象S(i0.5) 打出来的数值在部分区间明显跳变尤其在 i 6、7 附近曲线看起来「折了一下」。原因回代公式里同时出现 Mⱼ 和 Mⱼ₊₁MATLAB 里对应 M(j) 与 M(j1)。如果定位区间时用了 find(x xv, 1, last) 却把结果直接当成 j1 用或者循环里 j 从 1 跑到 n 而没有对右端点做保护区间和弯矩就会整体错一位只有端点附近碰巧对得上。解决写完求值函数先做自检——把查询点设成原节点S(xⱼ) 必须等于 yⱼ。这一条能抓住 90% 的下标错位。再检查 xq 打在 x(n1) 上时有没有越界加一句 if j n, j n; end 兜底。5.3 坑三二分循环没有上限程序卡死现象exp5.m 跑起来之后命令行一直刷 iter几分钟都不停最后只能 CtrlC。原因外推判据用的是相邻两次近似值之差小于 eps如果 eps 设得比机器精度还小比如 1e-16或者被积函数在积分区域内有接近奇点的行为差值永远降不到阈值以下。解决循环里加 maxit 上限并且每轮把 diff 打印出来。正常收敛时 diff 应当每轮下降大约一个数量级复化梯形二分一次误差降 4 倍外推后更快如果 diff 在某个量级来回震荡说明已经到浮点精度极限直接把当前值当结果别硬等。5.4 坑四复化梯形端点权重漏乘 1/2现象积分结果和参考值差一个小量而且二分次数增加后误差不按预期下降一直卡在 1e-3 左右。原因复化梯形公式里首尾两个节点的权重是 1/2中间节点是 1。循环写 for ix 0:nx 时如果忘记判断端点全部按 1 累加等价于多加了两端各半个函数值误差量级正好是 O(h)。解决在内外两层循环里都显式写端点判断或者更稳妥地用向量化写法w ones(nx1,1); w([1 end]) 0.5; s h * sum(w .* fvals)。写完之后拿 f x.*y 这种有闭式解的函数验证一遍误差应该随 h 减小而快速下降。5.5 坑五AB4 起步的索引从 0 还是从 1 开始现象AB4 的结果从第二步开始就明显偏离而 RK4 完全正常。原因数学推导里 x₀ 是第一个点f₀ 是第一步的斜率MATLAB 数组从 1 开始f(1) 对应数学上的 f₀。如果循环写 for i 4:N-1 却把 f1 取成 f(i-3)、f2 取成 f(i-2)整个差商窗口就整体右移了一位前几步看起来还行越往后偏得越多。解决写循环之前先在纸上列一张对照表把数学下标 0、1、2、3 和 MATLAB 下标 1、2、3、4 写清楚再动手。跑完之后单独打印前五步的 y 值跟 RK4 的结果逐点比前四步必须完全一致因为就是 RK4 算的第五步开始才允许有微小差异。6. 二次开发把三份作业拼成一个可复用的小型数值算法库三份作业各自能跑但接口风格不统一样条那边传 x、y、dy0、dynODE 这边传函数句柄和步长重积分又传四个边界加两个步长。想把这套东西变成以后能反复用的小工具第一件事是统一调用方式。我的做法是保留底层的单个算法文件再加一层薄薄的调度脚本用字符串选择方法返回值固定成结构体function res numlab(task, opt) % 统一入口task 取 spline | integral | ode switch task case spline M spline3(opt.x, opt.y, opt.dy0, opt.dyn); res.M M; res.sq evalS(opt.x, opt.y, M, opt.xq); % 查询点上的函数值 case integral res.I integrate2d(opt.f, opt.ab, opt.cd, opt.eps, opt.maxit); case ode res.y odeSolve(opt.f, opt.span, opt.y0, opt.h, opt.method); otherwise error(未知任务类型: %s, task); end endodeSolve 内部再用 switch 分发到 rk4、ab4am4 的不同 mode这样以后加一个新方法只需要在一个地方加分支。第二件事是把误差统计做成表格四种方法跑完直接出一张对照表比看一屏 fprintf 直观方法2-范数误差相对 RK4 的倍数每步函数调用次数RK43.8e-51.04改进型 AB4-AM45.31e-4约 142AB4-AM49.2e-4约 242AB44.809e-3约 1271这张表比单纯背结论有用RK4 精度最高但每步要算四次 fAB4 只要一次预测校正两次。实际工程里如果 f 的计算很贵比如每次都要解一个子问题AB4-AM4 这类多步法反而更划算精度和 RK4 差一个量级代价只有一半。第三件事是给每个函数加一个自检开关默认关掉调试时打开就跑一遍基准算例样条用线性函数 y 2x1 验证三次样条对线性函数应当精确复现积分用 x.*y 验证ODE 用 y′ y 验证。基准过了再上真实数据。从那以后我每次写完数值代码都强制先跑一遍能解析求解的基准算例再去看真实问题的结果。翻车的次数没少但至少能立刻分清是公式抄错了还是数据本身有问题排查时间从半天缩到十几分钟。这份上机资料里的三个主脚本和对应子函数都在包里接手之后建议先跑通 exp5 看输出格式再换自己的函数进去试最后再回头看样条和 ODE 那两块顺序反了容易被一堆下标绕晕。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
