数值分析上机实战:三次样条、重积分外推与ODE四法对比

📅 发布时间:2026/10/15 9:13:37
数值分析上机实战:三次样条、重积分外推与ODE四法对比
简介资源为东南大学研究生数值分析课程的上机作业文档作者用MATLAB完成三次样条插值函数、重积分计算与常微分方程初值问题数值解三个核心上机题。内容覆盖第一型三次样条插值程序编制与车门曲线数据验证、复化梯形公式结合外推思想计算重积分、RK4法、AB4法与AB4-AM4预测校正方法求解常微分方程初值问题并附计算误差对比和精度排序分析。文档包含完整算法推导、计算公式、程序实现思路及运行结果适合正在学习数值分析课程、需要参考MATLAB编程实现或准备相关上机实验的研究生和本科生。压缩包内为1个docx文档大小566KB结构紧凑便于查阅。已有904人浏览学习可作为课程作业写作、算法理解与编程实践的参考资料。文档对每个上机题都给出了问题描述、算法说明、计算结果和反思总结读者可对照公式逐步理解实现过程也可借鉴作者遇到的输入错误排查与子函数模块化设计经验。1. 数值分析上机作业三次样条、重积分与ODE四法对比一次跑通数值分析这门课的上机作业是理论最能露馅的环节课本上公式推得漂亮一落到 MATLAB 代码里下标、端点、初值全是坑。我拆的这份资源包含三道题——三次样条插值汽车门曲线、重积分复化梯形二分外推、常微分方程 RK4 与 Adams 系列四法对比正好覆盖插值、数值积分、ODE 数值解三个最常考的类型。适合正在做数值分析上机题的人、MATLAB 刚入门的研究生以及期末复习想找一套完整可跑模板的本科生。三道题的通用程序和结果直接对应作业要求照着改参数就行。2. 三次样条插值三弯矩方程的装配与追赶法求解2.1 第一型三次样条的公式选型为什么用三弯矩而不是三转角三次样条插值的目标是在每个小区间 [x_i, x_{i1}] 上构造三次多项式让整个区间上的函数值、一阶导数、二阶导数在节点处连续。教材里常见的做法有两套一套是三转角以节点处一阶导数为未知量另一套是三弯矩以节点处二阶导数 M_i 为未知量。这题端点条件是 y(0)0.8、y(10)0.2属于第一型端点条件三弯矩方程处理第一型边界最直接——只需在矩阵首尾各补一行端点方程整体仍然保持三对角追赶法一次解出 M 列向量。三转角那边边界会出现非三对角项代码处理麻烦不少。核心递推方程是教材里的黑匣子公式μ_j * M_{j-1} 2 * M_j α_j * M_{j1} d_j j 1, 2, ..., n-1其中μ_j h_{j-1} / (h_{j-1} h_j)α_j h_j / (h_{j-1} h_j) 1 - μ_jd_j 6 * f[x_{j-1}, x_j, x_{j1}]也就是节点 j 附近的三点二阶差商乘 6这题的 xi 是等距节点 0 到 10hj 全部等于 1所以 μ_j 和 α_j 都恒等于 0.5。等距条件下 d_j 还能再化简d_j 3 * (y_{j1} - 2*y_j y_{j-1})。简化之后矩阵装配过程就变成了纯粹的差分计算。2.2 三弯矩方程装配边界方程与矩阵结构的细节内部节点方程之外第一型边界条件贡献两个方程。左端点方程2*M_0 M_1 (6/h) * ( f[x_0, x_1] - y(0) )这里 f[x_0, x_1] (y_1 - y_0)/h等距 h1 时左端点方程化简为 2M_0 M_1 6(y_1 - y_0 - 0.8)。代入汽车门第一组数据 y02.51、y13.30得到 2M_0 M_1 6(3.30-2.51-0.8) -0.06。右端点方程同理M_{n-1} 2*M_n (6/h) * ( y(10) - f[x_{n-1}, x_n] )代入 y95.70、y105.80得到 M_9 2M_10 6(0.2 - 0.10) 0.6。把这些方程按 j0, 1, ..., n 排成三对角矩阵M 就是待求的弯矩向量。注意这里 j 的取值范围是 0 到 n矩阵是 (n1) 阶内部 j1 到 n-1 对应 n-1 个方程加上首尾两个边界方程总共 n1 个方程未知量 M_0 到 M_n 也是 n1 个刚好封闭。我在配矩阵的时候吃过一次亏把方程写成了 μ_j * M_{j-1} 2M_j α_j * M_{j-1} d_j左边 M 的下标全写成了 j-1结果矩阵完全退化。所以配矩阵前先确认对角线上是 2M_j两侧是 μ_jM_{j-1} 和 α_j*M_{j1}方向别反。2.3 MATLAB 实现追赶法解三对角与 S(i0.5) 输出先写一个通用追赶法求解函数三对角矩阵用三个对角向量表示避免存稀疏矩阵。这个函数可以复用一整类三对角问题function x tridiag(a, b, c, d) % 追赶法求解三对角方程组 % 输入: a 为下对角线(beta), b 为主对角线, c 为上对角线(gamma), d 为右端项 n length(b); x zeros(n,1); p zeros(n,1); q zeros(n,1); % 追: 消去下对角线 p(1) -c(1) / b(1); q(1) d(1) / b(1); for i 2:n denom b(i) a(i) * p(i-1); p(i) -c(i) / denom; q(i) (d(i) - a(i) * q(i-1)) / denom; end % 赶: 回代 x(n) q(n); for i n-1:-1:1 x(i) p(i) * x(i1) q(i); end end这个实现里 a、b、c 分别对应三对角矩阵的下、主、上对角线。追赶法的数值稳定性对严格对角占优矩阵很好这里 2M_j 的系数 2 显然占优不会出现除以接近零的数导致翻车。主程序如下% 三次样条插值主程序 clc; clear; % 汽车门曲线型值点 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]; % 端点条件(第一型) yp0 0.8; ypn 0.2; n length(x) - 1; h diff(x); % 各段步长, 等距情况下全为1 % 内部节点三弯矩方程系数 mu h(1:end-1) ./ (h(1:end-1) h(2:end)); alpha 1 - mu; % 计算二阶差商, d_j 6 * f[x_{j-1}, x_j, x_{j1}] d zeros(n-1, 1); for j 2:n f1 (y(j) - y(j-1)) / h(j-1); % 一阶差商 f[x_{j-1}, x_j] f2 (y(j1) - y(j)) / h(j); % 一阶差商 f[x_j, x_{j1}] f12 (f2 - f1) / (h(j-1) h(j)); % 二阶差商 d(j-1) 6 * f12; end % 组装三对角方程的系数向量 b 2 * ones(n1, 1); % 主对角线全为2 a zeros(n1, 1); % 下对角线 c zeros(n1, 1); % 上对角线 % 边界方程: 左端点 2M0 M1 a(2) 0; % 第1行无下对角 c(1) 1; dvec zeros(n1, 1); dvec(1) 6 * ((y(2)-y(1)) / h(1) - yp0); % 内部节点方程 for j 2:n a(j) mu(j-1); b(j) 2; c(j) alpha(j-1); dvec(j) d(j-1); end % 右端点方程 M_{n-1} 2M_n a(n1) 1; dvec(n1) 6 * (ypn - (y(n1)-y(n)) / h(n)); % 追赶法求解 M tridiag(a, b, c, dvec); % 计算 S(i0.5), i0..9 mid_vals zeros(1, n); for i 1:n % 中点处 x xi 0.5, 代三次样条公式 hi h(i); left M(i) * (0.5*hi)^3 / (6*hi) ... (y(i) - M(i)*hi^2/6) * 0.5; right M(i1) * (0.5*hi)^3 / (6*hi) ... (y(i1) - M(i1)*hi^2/6) * 0.5; mid_vals(i) left right; end disp(M ); disp(M); disp(S(i0.5) ); disp(mid_vals);代码逻辑分四层先算各段步长 h 和内部节点系数 μ、α再算三点二阶差商得到右端项 d然后组装三对角矩阵首尾补充端点条件最后追赶法求解 M。输出阶段求 S(i0.5) 时把 x_i 0.5 代进小区间三次表达式等距 h1 时可以写成中点加权形式。这个过程最容易错的就是边界方程的符号——左端点的 d 是一阶差商减去左端点导数值右端点刚好反过来写成左减右还是右减左必须和课本逐项核对。2.4 车门曲线数据的结果验证拟合出的 S(x) 在节点处必须严格等于 y_i这是硬约束同时端点处 S(x) 应分别接近 0.8 和 0.2。验证方法很简单把 i0.5、1.5、……、9.5 代进代码算出来的 mid_vals再看原题给出的目标值是否一致。节点处的函数值误差应该在 1e-14 量级如果有超过 1e-10 的偏差说明矩阵装配或差商计算有问题。另一个自检手段是把端点改成自然样条M_0M_n0看曲线端点是否变得平缓——变了说明程序能响应边界条件没变说明边界行没有真正参与求解。3. 重积分计算复化梯形、逐次二分与外推的工程化实现3.1 从复化梯形到 Simpson 的外推路径题目要求对二重积分 I(f) ∫c^d ∫a^b f(x,y) dx dy 用复化梯形公式逐次二分步长相邻两次近似值之差的绝对值不超过 ε0.5×10^-5 时停止然后利用外推思想算更精确的近似值。复化梯形公式的误差展开只含偶数阶项这意味着把步长减半后两次结果的误差比例是固定的。因此 T_n 和 T_2n 做线性组合S_n(f) (4/3) * T_2n(f) - (1/3) * T_n(f)这个组合能消掉误差展开中的 h² 项从二阶精度跳到四阶精度得到的就是复化 Simpson 公式的近似值。这是 Richardson 外推最基础的应用不需要额外增加求值点只靠两次梯形结果的加权就能白赚两阶精度。二重积分场景下外推在 x、y 两个方向同时减半执行形式上和一维完全一致前提是每次二分两个方向都同步加密。3.2 子函数拆分myfun 与 getT 的设计取舍这道题的代码组织特意拆成两个子函数getT.m 负责计算任意网格下的复化梯形值myfun.m 只做一件事——返回被积函数 f(x,y)tan(x²y²) 的值。这样拆分的好处是换题目时不碰主逻辑只改函数句柄。复化梯形的实现用矩阵化写法代替三重循环权重矩阵一次性生成function T getT(nx, ny, fun) % 二维复化梯形积分 % 输入: nx为x方向等分数, ny为y方向等分数, fun为函数句柄 % 输出: T为该网格下的复化梯形公式近似值 a 0; b pi/3; % 内层积分区间 c 0; d pi/6; % 外层积分区间 hx (b-a) / nx; hy (d-c) / ny; % 网格节点 x a : hx : b; y c : hy : d; % 梯形权重: 端点0.5, 内部1 wx ones(1, nx1); wx(1) 0.5; wx(end) 0.5; wy ones(1, ny1); wy(1) 0.5; wy(end) 0.5; % 被积函数矩阵 [X, Y] meshgrid(x, y); F fun(X, Y); % 二维梯形公式: 各维权重外积生成二维权重 W wy * wx; T hx * hy * sum(sum(W .* F)); end参数说明nx 和 ny 分别是 x、y 两个方向的等分数初始从 2 开始二分fun 是函数句柄通过 myfun 传入。权重外积 W wy * wx 的写法是关键——因为二维梯形本质是两个一维梯形的嵌套权重沿两个方向独立相乘。这里 meshgrid(x,y) 的输出矩阵 X 是 ny1 行 nx1 列Y 与之同型所以权重矩阵也要是相同维度直接用 wy 乘 wx 正好对齐。被积函数单独放一个 M 文件function z myfun(x, y) % 被积函数: f(x,y) tan(x^2 y^2) % 支持矩阵输入, 便于 meshgrid 后批量计算 z tan(x.^2 y.^2); end注意这里用了点乘 .^2 而不是矩阵乘 ^2。MATLAB 新手最容易在这翻车meshgrid 给进来的是矩阵x^2 会直接报维度错误或者在某些场景下默默算成矩阵乘法导致结果完全错误。所有涉及逐元素运算的地方必须带点。3.3 主程序二分循环与精度判断的顺序主程序 exp5.m 的循环控制逻辑决定最终精度。正确顺序是先算当前网格 T1再算加密一倍网格 T2比较差值是否超过 ε如果没超过把 T2 当作基准配合 T1 做外推得到 Simpson 值如果超过T1←T2继续加密。% exp5.m 主程序 clc; clear; epsilon 0.5e-5; % 允许误差 nx 2; ny 2; T1 getT(nx, ny, myfun); nx2 2*nx; ny2 2*ny; T2 getT(nx2, ny2, myfun); % 逐次二分, 直到相邻两次差值不超过 epsilon while abs(T2 - T1) epsilon T1 T2; nx nx2; ny ny2; nx2 2*nx; ny2 2*ny; T2 getT(nx2, ny2, myfun); end % 外推: 复化梯形组合为复化Simpson S 4/3 * T2 - 1/3 * T1; fprintf(T_%d %.10f\n, nx, T2); fprintf(S_%d %.10f\n, nx, S); fprintf(最终网格: nx%d, ny%d\n, nx2, ny2);精度判断放在二分之后、外推之前。如果把外推值 S 放进 while 循环当判断对象逻辑会乱——外推得到的已经是四阶精度拿它跟二阶的 T 系列比较没有意义。程序的效率瓶颈在 tan 函数求值每次网格加密求值点翻四倍二分五六轮之后网格密度已经足够运行时间仍然可控。如果发现计算特别慢优先检查是否忘记把被积函数写成向量化版本——逐点 foor 循环调用函数句柄是最慢的写法。4. 常微分方程初值问题RK4 与 AB 族四种方法的代码对比4.1 四套公式的差异与精度预期题目初值问题y -x²y²0≤x≤1.5y(0)3步长 h0.1精确解 y(x)3/(1x³)。要求实现 RK4、AB4Adams-Bashforth 四步显式、AB4-AM4 预测校正、改进的 AB4-AM4 预测校正共四种方法。RK4 是单步法每步需要四次函数求值局部截断误差 O(h⁵)。AB4 是多步法需要前四个点的函数值——这四个初始点由 RK4 提供之后计算成本大幅下降每步只求一次 f。AB4 显式方法的稳定性区间比 RK4 窄对步长更敏感。AB4-AM4 在显式预测后加一步隐式校正局部截断误差比单独 AB4 更小且每步需要两次函数求值预测一次、校正一次。改进的 AB4-AM4 增加了预测—校正—加权平均三步结构系数分母从 12 变为 24末了再把预测值和校正值按 19/270 与 251/270 加权。4.2 RK4 通用程序与初值生成function [x, y] rk4(f, a, b, h, y0) % 经典四阶Runge-Kutta方法 % 输入: f为右端函数句柄, a/b为区间端点, h为步长, y0为初始值 % 输出: x为节点向量, y为数值解向量 x a : h : b; n length(x); y zeros(size(x)); y(1) y0; for i 1:n-1 k1 f(x(i), y(i)); k2 f(x(i) h/2, y(i) h*k1/2); k3 f(x(i) h/2, y(i) h*k2/2); k4 f(x(i) h, y(i) h*k3); y(i1) y(i) h/6 * (k1 2*k2 2*k3 k4); end end四个 k 值的顺序不能换k1 是当前点斜率k2、k3 是中点斜率的两步不同估计k4 是终点斜率。如果以为 k2 和 k3 只差一个 k 值的替换关系而省略中间变量后续预测校正部分会出大问题。这个函数同时也是 AB4 的初值提供者——多步法在 i4 开始递推之前y(1) 到 y(4) 必须来自 RK4 保证精度用欧拉法或其他低阶方法提供初值会直接拉低整条曲线的精度。4.3 AB4、AB4-AM4 与改进型的实现细节三种多步法共用一个右端函数 ffunction dydx fxy(x, y) % 右端函数: y -x^2 * y^2 dydx -x.^2 .* y.^2; end精确解用于误差评估function yexact fxm(x) % 精确解: y 3 / (1 x^3) yexact 3 ./ (1 x.^3); end四种方法一次跑完并输出二范数误差% 四种方法对比主程序 clc; clear; a 0; b 1.5; h 0.1; y0 3; x a : h : b; n length(x); % RK4 [xr, yrk4] rk4(fxy, a, b, h, y0); % AB4: 前4点由RK4提供 yab4 zeros(1, n); yab4(1:4) yrk4(1:4); for i 4:n-1 f_i fxy(x(i), yab4(i)); f_im1 fxy(x(i-1), yab4(i-1)); f_im2 fxy(x(i-2), yab4(i-2)); f_im3 fxy(x(i-3), yab4(i-3)); yab4(i1) yab4(i) h/12 * (55*f_i - 59*f_im1 37*f_im2 - 9*f_im3); end % AB4-AM4预测校正 yam4 zeros(1, n); yam4(1:4) yrk4(1:4); for i 4:n-1 % 预测(P) f_i fxy(x(i), yam4(i)); f_im1 fxy(x(i-1), yam4(i-1)); f_im2 fxy(x(i-2), yam4(i-2)); f_im3 fxy(x(i-3), yam4(i-3)); yp yam4(i) h/12 * (55*f_i - 59*f_im1 37*f_im2 - 9*f_im3); % 校正(E) fp fxy(x(i1), yp); yam4(i1) yam4(i) h/12 * (9*fp 19*f_i - 5*f_im1 f_im2); end % 改进的AB4-AM4 (加权平均) yimp zeros(1, n); yimp(1:4) yrk4(1:4); for i 4:n-1 f_i fxy(x(i), yimp(i)); f_im1 fxy(x(i-1), yimp(i-1)); f_im2 fxy(x(i-2), yimp(i-2)); f_im3 fxy(x(i-3), yimp(i-3)); yp yimp(i) h/24 * (55*f_i - 59*f_im1 37*f_im2 - 9*f_im3); fc fxy(x(i1), yp); yc yimp(i) h/24 * (9*fc 19*f_i - 5*f_im1 f_im2); yimp(i1) (251*yc 19*yp) / 270; end % 误差比较 yexact fxm(x); err_rk4 norm(yrk4 - yexact); err_ab4 norm(yab4 - yexact); err_am4 norm(yam4 - yexact); err_imp norm(yimp - yexact); fprintf(RK4误差: %.6f\n, err_rk4); fprintf(AB4误差: %.6f\n, err_ab4); fprintf(AB4-AM4误差: %.6f\n, err_am4); fprintf(改进AB4-AM4: %.6f\n, err_imp);代码逻辑上普通 AB4-AM4 与改进型的区别集中在三个地方预测步分母是 12 还是 24校正步是否该有单独的预测值加权项最终结果是直接取校正值还是与预测值按 251/270 与 19/270 加权。改进型这样选择的原因是单独用校正值会保留部分误差结构而把预测和校正按固定比例混合能抵消更多低频误差分量——这本质上也是一种外推思想只是作用在时间轴上。多步法的迭代从 i4 开始这是 MATLAB 数组下标 1-based 下对应第 5 个点x0.4 处。如果从 i5 开始相当于跳过一个有效递推步会导致 x0.3 到 x0.4 这个区间没有推进。4.4 四法对比结论精度排序与原因原题运行结果中四种方法的二范数误差如下方法误差二范数精度名次RK40.0000381改进的 AB4-AM40.0005312AB4-AM40.0009203AB40.0048094RK4 精度最高的原因是每个步长内做了四次斜率估计等效 Richardson 外推到了四阶且单步构造没有任何历史点依赖。AB4 纯显式多步法最差是因为四个历史点引入的插值误差没有校正环节去抑制。改进型能明显优于普通预测校正关键在于最终输出的加权平均把预测误差和校正误差做了部分对消代价是每步多算一次 f 值。这个排序不是绝对的——如果步长继续减半AB4 的误差会按四阶速度下降与 RK4 的差距会缩小但在这题 h0.1 的设置下RK4 优势明显。5. 上机避坑三次样条、重积分与 ODE 的五个翻车现场5.1 三弯矩方程下标注错M_{j1} 写成 M_{j-1}现象程序跑通节点处插值结果正确但区间中点的 S(i0.5) 数值与预期曲线形状明显不对某些段甚至出现折线感。原因教材公式是 μ_j M_{j-1} 2M_j α_j M_{j-1} d_j 的原稿排版里存在笔误右侧 M 下标应当是 j1照抄后矩阵变成双下三角结构丢失了上三角信息。我看的这份作业资源里正好踩了这个坑找同学帮忙运行对比才定位到。解决装配矩阵前画一个 3×3 示意小矩阵手动填几行确定非零位置。三对角矩阵的下标规律是行 j 的三个非零列分别是 j-1、j、j1对应系数 μ_j、2、α_j。边界行第 0 行只有 2M_0 M_1 两个非零最后一行的 M_{n-1} 2M_n 同理。5.2 d_j 漏乘 6 或差商公式用反现象弯矩 M 的量级明显偏小插值曲线整体偏离数据点但节点处仍然吻合。原因是样条方程右端项 d_j 6 乘以二阶差商漏乘后相当于把弯矩系统整体缩小 6 倍节点值误差为 0 是因为插值条件强制了节点相等但节点之间的曲率全错了。原因公式记忆混乱把一阶差商和二阶差商混在一个式子里。解决验证方法很粗暴——计算一个等距三点y0, y1, y2的 d_j如果 x 距离为 1手算二阶差商 (y2 - 2y1 y0)/2d_j 应该是它的 6 倍即 3(y2 - 2*y1 y0)。先用这种三个点的小案例单测 d_j 再跑全流程。5.3 外推时机错误先外推后判断还是先判断后外推现象重积分结果精度不稳定有时候正好达标有时候需要多二分两轮才满足输出的 S 值与 T 值差距不符合四阶外推预期。原因把 S_n 4/3 T_2n - 1/3 T_n 的公式套进了 while 循环判断条件拿 Simpson 值跟上一轮的 T 比较两个不同精度的东西自然没法形成稳定收敛判据。解决判断基准只使用同阶量——T2 与 T1 都是复化梯形值只是网格密度不同。循环跳出后才做外推一次组合即可。外推输出后如果想继续加密获得更高精度需要把 S 当作新的基准再做一轮 Simpson 外推但不适用于这道题的停止条件。5.4 多步法初值来源错误用欧拉法提供 AB4 初始点现象RK4 误差 1e-4 量级AB4 误差却达到 1e-2比预想大了一个数量级且 0.4 之前段与之后段误差分布明显不同。原因AB4 的启动点如果用欧拉法或改进欧拉法前四点的局部误差就有 O(h²) 量级多步递推把这份误差以四阶方法的方式继续传播——前面差一点后面跟着差一大截。RK4 提供初值误差是 O(h⁵)完全匹配 AB4 的精度档位。解决启动代码yab4(1:4) yrk4(1:4)保证前四点是四阶精度。同理AB4-AM4 和改动型的初始点也必须统一从 RK4 复制不能各自用不同方法生成——否则四种方法的初值不一致后面的误差对比毫无意义。5.5 meshgrid 求解时权重矩阵维度忘记转置现象getT 返回的积分值远大于实际值或者随着网格加密不收敛反而震荡。检查 sum(sum(W .* F)) 时报矩阵维度不匹配或匹配但大小不对。原因meshgrid(x,y) 生成的 X、Y 矩阵维度是 (ny1)×(nx1)而 wx 是 1×(nx1) 行向量wy 是 (ny1)×1 列向量时权重矩阵才能对齐。写成 wx * wy 得到的是 (nx1)×(ny1)方向正好相反。解决确认W wy * wx其中 wy 的长度等于 ny1 对应行数wx 长度等于 nx1 对应列数。最稳妥的方式是size(W)与size(F)完全一致再继续。6. 精度验证的复用技巧用精确解的收敛阶反查代码错误三道题里ODE 部分有精确解 y(x) 3/(1x³)这是验证方法实现正确性的黄金标尺。我拿到代码后做的第一件事不是直接跑题目而是临时把步长从 0.1 改成 0.05 和 0.025看四种方法的误差如何变化。四阶方法的收敛表现应当是步长减半误差缩小约 16 倍因为 h⁴ 项主导。如果实测 RK4 从 0.000038 降到 0.0000027 左右缩小约 14~18 倍说明实现正确误差主要来自截断误差而非程序 bug。如果缩小倍数只有 4 倍左右说明实际只有二阶精度——八成是某个 k 值写错了或者多步法某一步把 f_i 和 f_{i-1} 下标错位了。这个验证思路不依赖理论分析完全靠数值实验就能定位问题。三次样条部分同理没有精确解但可以用已知的简单函数做反向测试。比如用 y sin(x) 在 0 到 10 上取节点做插值端点导数用解析值代进去跑完程序看 S 在 0.5、1.5 这些中点处与 sin 真实值的误差。如果误差量级在 1e-4 到 1e-3对应四次样条逼近的典型表现说明程序正确如果误差是 1e-1说明矩阵装错了。这种用已知函数当黑匣子探针的办法比逐行读代码排查效率高得多。从那以后我每次拿到这类数值分析上机资源都先强制自己走一遍已知解反查 步长减半看收敛阶的流程再跑正式题目。这一套验证下来公式输入错误基本十分钟内暴露。希望帮到你。本文还有配套的精品资源点击获取