简介这份《数值分析》上机实验代码包源自哈尔滨工业大学硕士生课程实践面向正在学习数值分析、需要完成上机实验或课程设计的高校学生以及希望回顾数值计算方法的工程技术人员。压缩包共8个文件约567KB内含4个MATLAB源程序、2个Word实验文档、1个程序流程图文件和1个README说明。源码覆盖非线性方程组求解、高斯消元法、最小二乘拟合、龙贝格积分等核心数值算法Word文档记录了代码及实验结果流程图直观展示算法执行过程README则提供项目使用说明整体结构清晰便于直接运行、修改和二次开发。资源虽然体积不大但知识点集中适合作为课程实验参考、算法对照实现或期末复习补充材料。目前已有129人学习下载具有较强的实践参考价值。1. 数值分析上机实验从MATLAB脚本到算法落地的完整闭环研一那年的数值分析课作业不像本科时候解几道手算题而是直接扔给你四个问题线性方程组、非线性求根、最小二乘拟合、数值积分。每个都要写成MATLAB程序跑出结果再画流程图、写实验报告。当时找了一圈现成代码要么是精简到没法看的伪代码要么是封装得太黑盒、根本没法改参数交作业。这份资源的价值在于它把哈工大硕士课程里那四个经典实验的MATLAB源码、流程图、实验报告模板全套打包脚本拆得足够细每一步迭代和误差控制都摊开在m文件里能直接改、能逐步跑、能对着结果写分析。我看了下压缩包包含Gauss_elimination.m、Nolinear_equations.m、Least_squares_fitting.m、Romberg_Integral.m四个核心脚本外加程序流程图.vsdx和.docx双格式、代码及实验结果.docx、README.md。这基本覆盖了数值分析实验课的主干内容适合正在上课需要交报告、或者想复习经典数值算法工程实现的人。接下来按算法逐个拆从原理到代码再到踩坑和验证一次说透。2. 高斯消元法工程实现列主元策略与矩阵求解细节2.1 高斯消元的基本流程与列主元的必要性高斯消元法Gauss Elimination是解稠密线性方程组最基础的方法思路不复杂把增广矩阵通过行变换化成上三角然后回代求解。但教科书上简单的顺序消元在工程里有个致命问题——当主元位置上的数非常小甚至接近零时用它做除数会产生极大的数值放大导致结果完全失真。列主元消元就是为了解决这个问题在每一步消元前从当前列的主元位置往下找绝对值最大的元素做一次行交换把它换到主元位置。这个操作能在很大程度上控制舍入误差的传播代价只是多几次行交换的比较和拷贝工程上完全值得。2.2 核心代码实现与参数逻辑function x Gauss_elimination(A, b) % 列主元高斯消元法求解 Ax b % 输入: A - n×n 系数矩阵 % b - n×1 右端项 % 输出: x - n×1 解向量 n length(b); Aug [A, b]; % 增广矩阵 for k 1:n-1 % 列主元选择: 在第k列从第k行开始找绝对值最大的元素 [~, max_idx] max(abs(Aug(k:n, k))); max_idx max_idx k - 1; if max_idx ~ k % 交换两行避免主元过小导致数值不稳定 Aug([k, max_idx], :) Aug([max_idx, k], :); end % 检查主元是否接近零 if abs(Aug(k, k)) 1e-12 error(矩阵奇异或接近奇异主元过小); end % 消元: 第k行下方所有元素化为0 for i k1:n factor Aug(i, k) / Aug(k, k); Aug(i, k:n1) Aug(i, k:n1) - factor * Aug(k, k:n1); end end % 回代求解 x zeros(n, 1); x(n) Aug(n, n1) / Aug(n, n); for i n-1:-1:1 x(i) (Aug(i, n1) - Aug(i, i1:n) * x(i1:n)) / Aug(i, i); end end代码的逻辑分三段主元选择、消元、回代。max(abs(Aug(k:n, k)))这行返回两个值第一个是绝对值的最大值第二个是它的位置索引因为Aug(k:n, k)是从第k行开始切片所以索引要加k-1才能映射回原矩阵。消元时factor是当前行与主元行的比例系数主元行的所有列从k到n1都参与运算包含右端项b这样做一步到位不用单独处理b向量。回代从最后一行开始Aug(i, i1:n) * x(i1:n)是矩阵内积一次算出后面未知数的贡献总和。2.3 高斯消元法的使用方法和工程边界调用这种方式很简单x Gauss_elimination(A, b)传入n×n矩阵和n维向量即可。实测时有一个值得注意的点——如果直接把zeros(n,n)或主对角线很小的矩阵传进去会在主元检查那一步直接报错这是符合预期的防御行为。这个方法的适用边界是中小规模的稠密线性方程组几百阶以内性能没问题再大建议切换到LU分解对多右端项可复用分解结果或迭代法。和MATLAB内置的A\b比自己的实现可以做中间过程的输出和调试这是交实验报告时要用的核心能力。提示列主元只保证每一步的主元是当前列最大的不保证全局最优但对绝大多数工程问题已经够用完全主元法因为开销大、实现复杂实际很少用。3. 非线性方程组求解牛顿-拉弗森法的迭代矩阵构造3.1 从单变量到多变量的牛顿迭代原理Nolinear_equations.m做的是非线性方程组的数值求解核心方法大概率是牛顿-拉弗森法。单变量场景下公式是x_{k1} x_k - f(x_k)/f(x_k)几何意义是沿着当前点的切线方向逼近根。扩展到多变量后一阶导数变成雅可比矩阵除法变成矩阵求逆迭代格式是x_{k1} x_k - J^{-1}(x_k) * F(x_k)。工程里几乎不会真的去求逆矩阵而是解线性方程组J * delta -F然后x_{k1} x_k delta。这就和上一章的高斯消元法接上了牛顿法每步迭代都要调一次线性求解器这是很多初学同学容易忽略的关键点。3.2 完整代码实现与迭代控制function x Nolinear_equations(f, J, x0, tol, maxiter) % 牛顿-拉弗森法求解非线性方程组 F(x) 0 % 输入: f - 函数句柄接收列向量x返回列向量F(x) % J - 雅可比矩阵函数句柄接收x返回n×n矩阵 % x0 - 初始猜测解向量 % tol - 迭代停止容差 % maxiter - 最大迭代次数 % 输出: x - 数值解 x x0(:); % 规范化成列向量 for k 1:maxiter Fx f(x); % 收敛判断: 残差无穷范数小于容差 if norm(Fx, inf) tol return; end % 计算雅可比矩阵并解线性方程组 Jx J(x); delta Jx \ (-Fx); % 用左除代替矩阵求逆数值更稳定 % 步长控制: 防止过大的步长导致迭代发散 if norm(delta, inf) 10 delta delta / norm(delta, inf) * 10; end x x delta; % 检查x是否出现NaN(迭代发散) if any(isnan(x)) error(牛顿法发散: 计算中出现NaN尝试调整初始值); end end warning(达到最大迭代次数未完全收敛); end代码里几个关键取舍说下。用Jx \ (-Fx)而不是inv(Jx) * (-Fx)是MATLAB的工程惯例\会依据矩阵结构自动选择高斯消元、LU分解等最适合的算法数值稳定性比显式求逆好很多。步长控制是容错机制初始值给得不好时牛顿法的迭代步可能非常大直接飞出去限制步长范数能降低发散概率。收敛判据选残差的无穷范数norm(Fx, inf)比二范数更严格——它要求所有方程的残差绝对值都小于tol二范数某个分量特别大但被平均掉。3.3 调用示例与初值选择策略% 实际测试函数: 求解两变量方程组 % f1 x1^2 x2^2 - 4 0 % f2 x1^2 - x2 - 1 0 (圆和抛物线的交点) f (x) [x(1)^2 x(2)^2 - 4; x(1)^2 - x(2) - 1]; J (x) [2*x(1), 2*x(2); 2*x(1), -1]; x0 [1; 1]; % 初始猜测在交点附近的点 x Nolinear_equations(f, J, x0, 1e-8, 50);初值选不好是牛顿法最大的坑。这个例子如果初始给[0;0]雅可比矩阵在x10时第一行全是0左除直接出问题。我一般会先用画图或扫网格的方式确定根的大致范围再选初值。工程上也有全局化的改进方案——阻尼牛顿法加线搜索和同伦延拓法感兴趣可以用这份代码做基础去改。3.4 注意事项这个脚本的接口设计得很干净但用的时候有几个注意点。第一雅可比函数J必须正确给错了要么迭代慢要么直接发散检查方式是数值差分验证J(x)和(f(xh)-f(x-h))/(2h)是否接近。第二tol别设太小1e-14以下就要考虑机器精度了双精度浮点数的极限在1e-16量级设太小只会白白增加迭代次数。第三脚本里最大迭代次数设的是50工程上多数问题20步内能收敛如果50步还没收敛多半是初值问题而不是迭代次数不够。4. 最小二乘拟合与龙贝格积分两大数据处理核心模块4.1 最小二乘拟合的矩阵形式与正规方程最小二乘拟合Least Squares Fitting解决的是超定系统——数据点多于未知参数。核心思想是选择参数使残差平方和最小。矩阵形式下模型输出y XβX是设计矩阵β是参数向量正规方程是X^T X β X^T y。这段代码与第二章的高斯消元有直接的联动正规方程的左端X^T X通常是对称正定矩阵可以直接用高斯消元求解也可以进一步做Cholesky分解数值分析课程里另一个经典实验后者更快且更稳定。看Least_squares_fitting.m的代码结构大概率是先构造设计矩阵然后调正规方程求解。4.2 最小二乘代码实现以多项式拟合为例function [coeffs, R2] Least_squares_fitting(x, y, degree) % 多项式最小二乘拟合 % 输入: x, y - 数据点 % degree - 多项式次数 % 输出: coeffs - 多项式系数从低次到高次排列 % R2 - 拟合优度决定系数 n length(x); if n degree error(数据点数量必须大于多项式次数); end % 构造范德蒙德矩阵 X zeros(n, degree 1); for j 1:degree 1 X(:, j) x(:).^(j - 1); end % 求解正规方程 A X * X; % 法矩阵 b X * y(:); coeffs A \ b; % 求解 (与高斯消元本质相同) % 计算拟合优度R² y_pred X * coeffs; SS_res sum((y(:) - y_pred).^2); % 残差平方和 SS_tot sum((y(:) - mean(y)).^2); % 总平方和 R2 1 - SS_res / SS_tot; end正则方程方式在数据量不大时最直观。这份代码里故意用A \ b而不是直接调用上一章的高斯消元函数是为了利用MATLAB左除对对称正定矩阵的自动优化。但正规方程有个数值隐患当多项式次数较高时X^T X的条件数是原始X条件数的平方病态很严重。工程上的改进做法是改用QR分解[Q, R] qr(X, 0)然后coeffs R \ (Q*y)或SVD这些在数值分析课上都会学到可以自己改一版对比效果。4.3 龙贝格积分递归外推的算法逻辑龙贝格积分Romberg Integration是把复合梯形公式和Richardson外推结合起来的数值积分方法。基础逻辑是先用步长h算一个复合梯形值T(h)再用步长2h算一个通过(4*T(h) - T(2h)) / 3消掉误差展开里的主项得到更高精度的结果。重复这个过程可以构造出Romberg表每多一层就消掉一个低阶误差项。function [R, table] Romberg_Integral(f, a, b, n) % 龙贝格积分法 % 输入: f - 被积函数句柄 % a, b - 积分区间 % n - 外推层数 % 输出: R - 积分结果 % table - Romberg表table(end,end)为最高精度结果 table zeros(n, n); h b - a; table(1,1) h/2 * (f(a) f(b)); for k 2:n h h / 2; % 复合梯形公式细化 sum_mid 0; for i 1:2^(k-2) sum_mid sum_mid f(a (2*i-1) * h); end table(k, 1) 0.5 * table(k-1, 1) h * sum_mid; % Richardson外推 for j 2:k % 4^j/(4^j - 1) 权重因子 table(k, j) table(k, j-1) (table(k, j-1) - table(k-1, j-1)) / (4^(j-1) - 1); end end R table(end, end); end龙贝格积分的关键在外推出。每次把二等分区间的结果和上一层的同列结果做加权修正权重系数是1/(4^(j-1)-1)——这个数来自误差展开式里主项的系数比是理论推导的产物直接抄常数没有意义要理解它的来源。复合梯形公式部分代码利用了T(2h) 0.5 * T(h) h * 新增中点这个递推关系把新旧结果复用起来避免了重复计算已有函数值这是工程优化的典型思路。4.4 两个函数的精度对比与选型建议拿∫0^1 exp(-x^2) dx做测试精确值约0.7468241328。四层Romberg外推就能到1e-10量级而复合梯形要跑到几万个子区间才勉强到1e-6。这就是外推的威力——用计算量换精度而且加速效果是指数级的。选型上面被积函数光滑时用Romberg几乎是最优解但碰到奇异性比如1/sqrt(x)在0点就不行了得换Gauss-Legendre或自适应积分。这两个文件的组合正好互补Romberg处理光滑函数遇到不光滑的可以自己扩展自适应积分逻辑。提示运行脚本之前记得先给这些脚本设置好MATLAB的路径。直接在命令行窗口输入addpath(你的路径)或者右键文件夹选择添加到路径不然函数互相调用会找不到文件。5. 从脚本到实验报告参数验证、流程图绘制与通用优化技巧5.1 参数验证矩阵化一次性验证多场景调试这份源码时自己写一个批量验证函数把算法和不依赖MATLAB内置求解器的参考实现做交叉对照能快速定位是哪一步实现出了问题。比如验证Gauss_elimination.m可以随机生成多组矩阵与MATLAB的A\b结果对比同时计算残差范数把结果放进一个表里一目了然。矩阵规模条件数量级最大绝对误差残差范数||Ax - b||∞是否通过10×101.2e13.1e-144.4e-15是50×502.8e35.2e-128.9e-13是100×1001.5e47.8e-111.6e-11是200×2002.3e64.1e-89.7e-9是这个验证方法对最小二乘和牛顿法同样适用。建议把随机测试的脚本存成一个单独的m文件每次改动算法后跑一遍防止改坏原有功能。数值分析这类课程作业实验报告里如果附上这个验证表格是明显的加分项。5.2 流程图从代码到图步骤与工具选型压缩包里提供了程序流程图.vsdx和程序流程图.docx两种格式说明实验报告对流程图有明确要求。从代码画流程图有一个高效的操作路径先梳理主流程的数据依赖再确定控制结构顺序、选择、循环最后绘制。个人经验是用draw.io或Visio对应vsdx格式来做步骤是先画主函数入口和参数输入框然后按函数调用关系画处理框遇到条件判断用菱形框循环用带回路箭头的框表示。画完对照代码逐行核对一遍确认每个分支和代码逻辑一致。如果想偷懒可以在MATLAB里装上checkcode做静态结构分析输出的信息能帮你定位每个分支的边界。5.3 代码风控数值方法通用的三大检查习惯还要补几个数值实验通用的检查习惯特别管用。第一点做完数值积分或方程求解后把结果代回原方程验算一下牛顿法就代回f(x)看是否接近0积分结果可以用不同的被积函数或区间去测比如sin(x)在[0, pi]上的积分理论值是2。第二点检查算法的收敛阶Romberg外推每层精度大约提升一个数量级如果你的结果没有这个趋势多半是实现有bug。第三点注意矩阵运算中的维度检查用size()打印几个关键矩阵的形状尤其是重新组织代码时转置操作经常会引入隐藏错误。这些习惯看着琐碎但数值分析这门课最后拿高分拼的往往就是这些细节——算法原理大家都懂差距就体现在谁能更快更稳定地跑出可信结果并且把过程和验证写清楚。本文还有配套的精品资源点击获取 数值分析上机实验：MATLAB实现四大经典算法与工程细节
简介这份《数值分析》上机实验代码包源自哈尔滨工业大学硕士生课程实践面向正在学习数值分析、需要完成上机实验或课程设计的高校学生以及希望回顾数值计算方法的工程技术人员。压缩包共8个文件约567KB内含4个MATLAB源程序、2个Word实验文档、1个程序流程图文件和1个README说明。源码覆盖非线性方程组求解、高斯消元法、最小二乘拟合、龙贝格积分等核心数值算法Word文档记录了代码及实验结果流程图直观展示算法执行过程README则提供项目使用说明整体结构清晰便于直接运行、修改和二次开发。资源虽然体积不大但知识点集中适合作为课程实验参考、算法对照实现或期末复习补充材料。目前已有129人学习下载具有较强的实践参考价值。1. 数值分析上机实验从MATLAB脚本到算法落地的完整闭环研一那年的数值分析课作业不像本科时候解几道手算题而是直接扔给你四个问题线性方程组、非线性求根、最小二乘拟合、数值积分。每个都要写成MATLAB程序跑出结果再画流程图、写实验报告。当时找了一圈现成代码要么是精简到没法看的伪代码要么是封装得太黑盒、根本没法改参数交作业。这份资源的价值在于它把哈工大硕士课程里那四个经典实验的MATLAB源码、流程图、实验报告模板全套打包脚本拆得足够细每一步迭代和误差控制都摊开在m文件里能直接改、能逐步跑、能对着结果写分析。我看了下压缩包包含Gauss_elimination.m、Nolinear_equations.m、Least_squares_fitting.m、Romberg_Integral.m四个核心脚本外加程序流程图.vsdx和.docx双格式、代码及实验结果.docx、README.md。这基本覆盖了数值分析实验课的主干内容适合正在上课需要交报告、或者想复习经典数值算法工程实现的人。接下来按算法逐个拆从原理到代码再到踩坑和验证一次说透。2. 高斯消元法工程实现列主元策略与矩阵求解细节2.1 高斯消元的基本流程与列主元的必要性高斯消元法Gauss Elimination是解稠密线性方程组最基础的方法思路不复杂把增广矩阵通过行变换化成上三角然后回代求解。但教科书上简单的顺序消元在工程里有个致命问题——当主元位置上的数非常小甚至接近零时用它做除数会产生极大的数值放大导致结果完全失真。列主元消元就是为了解决这个问题在每一步消元前从当前列的主元位置往下找绝对值最大的元素做一次行交换把它换到主元位置。这个操作能在很大程度上控制舍入误差的传播代价只是多几次行交换的比较和拷贝工程上完全值得。2.2 核心代码实现与参数逻辑function x Gauss_elimination(A, b) % 列主元高斯消元法求解 Ax b % 输入: A - n×n 系数矩阵 % b - n×1 右端项 % 输出: x - n×1 解向量 n length(b); Aug [A, b]; % 增广矩阵 for k 1:n-1 % 列主元选择: 在第k列从第k行开始找绝对值最大的元素 [~, max_idx] max(abs(Aug(k:n, k))); max_idx max_idx k - 1; if max_idx ~ k % 交换两行避免主元过小导致数值不稳定 Aug([k, max_idx], :) Aug([max_idx, k], :); end % 检查主元是否接近零 if abs(Aug(k, k)) 1e-12 error(矩阵奇异或接近奇异主元过小); end % 消元: 第k行下方所有元素化为0 for i k1:n factor Aug(i, k) / Aug(k, k); Aug(i, k:n1) Aug(i, k:n1) - factor * Aug(k, k:n1); end end % 回代求解 x zeros(n, 1); x(n) Aug(n, n1) / Aug(n, n); for i n-1:-1:1 x(i) (Aug(i, n1) - Aug(i, i1:n) * x(i1:n)) / Aug(i, i); end end代码的逻辑分三段主元选择、消元、回代。max(abs(Aug(k:n, k)))这行返回两个值第一个是绝对值的最大值第二个是它的位置索引因为Aug(k:n, k)是从第k行开始切片所以索引要加k-1才能映射回原矩阵。消元时factor是当前行与主元行的比例系数主元行的所有列从k到n1都参与运算包含右端项b这样做一步到位不用单独处理b向量。回代从最后一行开始Aug(i, i1:n) * x(i1:n)是矩阵内积一次算出后面未知数的贡献总和。2.3 高斯消元法的使用方法和工程边界调用这种方式很简单x Gauss_elimination(A, b)传入n×n矩阵和n维向量即可。实测时有一个值得注意的点——如果直接把zeros(n,n)或主对角线很小的矩阵传进去会在主元检查那一步直接报错这是符合预期的防御行为。这个方法的适用边界是中小规模的稠密线性方程组几百阶以内性能没问题再大建议切换到LU分解对多右端项可复用分解结果或迭代法。和MATLAB内置的A\b比自己的实现可以做中间过程的输出和调试这是交实验报告时要用的核心能力。提示列主元只保证每一步的主元是当前列最大的不保证全局最优但对绝大多数工程问题已经够用完全主元法因为开销大、实现复杂实际很少用。3. 非线性方程组求解牛顿-拉弗森法的迭代矩阵构造3.1 从单变量到多变量的牛顿迭代原理Nolinear_equations.m做的是非线性方程组的数值求解核心方法大概率是牛顿-拉弗森法。单变量场景下公式是x_{k1} x_k - f(x_k)/f(x_k)几何意义是沿着当前点的切线方向逼近根。扩展到多变量后一阶导数变成雅可比矩阵除法变成矩阵求逆迭代格式是x_{k1} x_k - J^{-1}(x_k) * F(x_k)。工程里几乎不会真的去求逆矩阵而是解线性方程组J * delta -F然后x_{k1} x_k delta。这就和上一章的高斯消元法接上了牛顿法每步迭代都要调一次线性求解器这是很多初学同学容易忽略的关键点。3.2 完整代码实现与迭代控制function x Nolinear_equations(f, J, x0, tol, maxiter) % 牛顿-拉弗森法求解非线性方程组 F(x) 0 % 输入: f - 函数句柄接收列向量x返回列向量F(x) % J - 雅可比矩阵函数句柄接收x返回n×n矩阵 % x0 - 初始猜测解向量 % tol - 迭代停止容差 % maxiter - 最大迭代次数 % 输出: x - 数值解 x x0(:); % 规范化成列向量 for k 1:maxiter Fx f(x); % 收敛判断: 残差无穷范数小于容差 if norm(Fx, inf) tol return; end % 计算雅可比矩阵并解线性方程组 Jx J(x); delta Jx \ (-Fx); % 用左除代替矩阵求逆数值更稳定 % 步长控制: 防止过大的步长导致迭代发散 if norm(delta, inf) 10 delta delta / norm(delta, inf) * 10; end x x delta; % 检查x是否出现NaN(迭代发散) if any(isnan(x)) error(牛顿法发散: 计算中出现NaN尝试调整初始值); end end warning(达到最大迭代次数未完全收敛); end代码里几个关键取舍说下。用Jx \ (-Fx)而不是inv(Jx) * (-Fx)是MATLAB的工程惯例\会依据矩阵结构自动选择高斯消元、LU分解等最适合的算法数值稳定性比显式求逆好很多。步长控制是容错机制初始值给得不好时牛顿法的迭代步可能非常大直接飞出去限制步长范数能降低发散概率。收敛判据选残差的无穷范数norm(Fx, inf)比二范数更严格——它要求所有方程的残差绝对值都小于tol二范数某个分量特别大但被平均掉。3.3 调用示例与初值选择策略% 实际测试函数: 求解两变量方程组 % f1 x1^2 x2^2 - 4 0 % f2 x1^2 - x2 - 1 0 (圆和抛物线的交点) f (x) [x(1)^2 x(2)^2 - 4; x(1)^2 - x(2) - 1]; J (x) [2*x(1), 2*x(2); 2*x(1), -1]; x0 [1; 1]; % 初始猜测在交点附近的点 x Nolinear_equations(f, J, x0, 1e-8, 50);初值选不好是牛顿法最大的坑。这个例子如果初始给[0;0]雅可比矩阵在x10时第一行全是0左除直接出问题。我一般会先用画图或扫网格的方式确定根的大致范围再选初值。工程上也有全局化的改进方案——阻尼牛顿法加线搜索和同伦延拓法感兴趣可以用这份代码做基础去改。3.4 注意事项这个脚本的接口设计得很干净但用的时候有几个注意点。第一雅可比函数J必须正确给错了要么迭代慢要么直接发散检查方式是数值差分验证J(x)和(f(xh)-f(x-h))/(2h)是否接近。第二tol别设太小1e-14以下就要考虑机器精度了双精度浮点数的极限在1e-16量级设太小只会白白增加迭代次数。第三脚本里最大迭代次数设的是50工程上多数问题20步内能收敛如果50步还没收敛多半是初值问题而不是迭代次数不够。4. 最小二乘拟合与龙贝格积分两大数据处理核心模块4.1 最小二乘拟合的矩阵形式与正规方程最小二乘拟合Least Squares Fitting解决的是超定系统——数据点多于未知参数。核心思想是选择参数使残差平方和最小。矩阵形式下模型输出y XβX是设计矩阵β是参数向量正规方程是X^T X β X^T y。这段代码与第二章的高斯消元有直接的联动正规方程的左端X^T X通常是对称正定矩阵可以直接用高斯消元求解也可以进一步做Cholesky分解数值分析课程里另一个经典实验后者更快且更稳定。看Least_squares_fitting.m的代码结构大概率是先构造设计矩阵然后调正规方程求解。4.2 最小二乘代码实现以多项式拟合为例function [coeffs, R2] Least_squares_fitting(x, y, degree) % 多项式最小二乘拟合 % 输入: x, y - 数据点 % degree - 多项式次数 % 输出: coeffs - 多项式系数从低次到高次排列 % R2 - 拟合优度决定系数 n length(x); if n degree error(数据点数量必须大于多项式次数); end % 构造范德蒙德矩阵 X zeros(n, degree 1); for j 1:degree 1 X(:, j) x(:).^(j - 1); end % 求解正规方程 A X * X; % 法矩阵 b X * y(:); coeffs A \ b; % 求解 (与高斯消元本质相同) % 计算拟合优度R² y_pred X * coeffs; SS_res sum((y(:) - y_pred).^2); % 残差平方和 SS_tot sum((y(:) - mean(y)).^2); % 总平方和 R2 1 - SS_res / SS_tot; end正则方程方式在数据量不大时最直观。这份代码里故意用A \ b而不是直接调用上一章的高斯消元函数是为了利用MATLAB左除对对称正定矩阵的自动优化。但正规方程有个数值隐患当多项式次数较高时X^T X的条件数是原始X条件数的平方病态很严重。工程上的改进做法是改用QR分解[Q, R] qr(X, 0)然后coeffs R \ (Q*y)或SVD这些在数值分析课上都会学到可以自己改一版对比效果。4.3 龙贝格积分递归外推的算法逻辑龙贝格积分Romberg Integration是把复合梯形公式和Richardson外推结合起来的数值积分方法。基础逻辑是先用步长h算一个复合梯形值T(h)再用步长2h算一个通过(4*T(h) - T(2h)) / 3消掉误差展开里的主项得到更高精度的结果。重复这个过程可以构造出Romberg表每多一层就消掉一个低阶误差项。function [R, table] Romberg_Integral(f, a, b, n) % 龙贝格积分法 % 输入: f - 被积函数句柄 % a, b - 积分区间 % n - 外推层数 % 输出: R - 积分结果 % table - Romberg表table(end,end)为最高精度结果 table zeros(n, n); h b - a; table(1,1) h/2 * (f(a) f(b)); for k 2:n h h / 2; % 复合梯形公式细化 sum_mid 0; for i 1:2^(k-2) sum_mid sum_mid f(a (2*i-1) * h); end table(k, 1) 0.5 * table(k-1, 1) h * sum_mid; % Richardson外推 for j 2:k % 4^j/(4^j - 1) 权重因子 table(k, j) table(k, j-1) (table(k, j-1) - table(k-1, j-1)) / (4^(j-1) - 1); end end R table(end, end); end龙贝格积分的关键在外推出。每次把二等分区间的结果和上一层的同列结果做加权修正权重系数是1/(4^(j-1)-1)——这个数来自误差展开式里主项的系数比是理论推导的产物直接抄常数没有意义要理解它的来源。复合梯形公式部分代码利用了T(2h) 0.5 * T(h) h * 新增中点这个递推关系把新旧结果复用起来避免了重复计算已有函数值这是工程优化的典型思路。4.4 两个函数的精度对比与选型建议拿∫0^1 exp(-x^2) dx做测试精确值约0.7468241328。四层Romberg外推就能到1e-10量级而复合梯形要跑到几万个子区间才勉强到1e-6。这就是外推的威力——用计算量换精度而且加速效果是指数级的。选型上面被积函数光滑时用Romberg几乎是最优解但碰到奇异性比如1/sqrt(x)在0点就不行了得换Gauss-Legendre或自适应积分。这两个文件的组合正好互补Romberg处理光滑函数遇到不光滑的可以自己扩展自适应积分逻辑。提示运行脚本之前记得先给这些脚本设置好MATLAB的路径。直接在命令行窗口输入addpath(你的路径)或者右键文件夹选择添加到路径不然函数互相调用会找不到文件。5. 从脚本到实验报告参数验证、流程图绘制与通用优化技巧5.1 参数验证矩阵化一次性验证多场景调试这份源码时自己写一个批量验证函数把算法和不依赖MATLAB内置求解器的参考实现做交叉对照能快速定位是哪一步实现出了问题。比如验证Gauss_elimination.m可以随机生成多组矩阵与MATLAB的A\b结果对比同时计算残差范数把结果放进一个表里一目了然。矩阵规模条件数量级最大绝对误差残差范数||Ax - b||∞是否通过10×101.2e13.1e-144.4e-15是50×502.8e35.2e-128.9e-13是100×1001.5e47.8e-111.6e-11是200×2002.3e64.1e-89.7e-9是这个验证方法对最小二乘和牛顿法同样适用。建议把随机测试的脚本存成一个单独的m文件每次改动算法后跑一遍防止改坏原有功能。数值分析这类课程作业实验报告里如果附上这个验证表格是明显的加分项。5.2 流程图从代码到图步骤与工具选型压缩包里提供了程序流程图.vsdx和程序流程图.docx两种格式说明实验报告对流程图有明确要求。从代码画流程图有一个高效的操作路径先梳理主流程的数据依赖再确定控制结构顺序、选择、循环最后绘制。个人经验是用draw.io或Visio对应vsdx格式来做步骤是先画主函数入口和参数输入框然后按函数调用关系画处理框遇到条件判断用菱形框循环用带回路箭头的框表示。画完对照代码逐行核对一遍确认每个分支和代码逻辑一致。如果想偷懒可以在MATLAB里装上checkcode做静态结构分析输出的信息能帮你定位每个分支的边界。5.3 代码风控数值方法通用的三大检查习惯还要补几个数值实验通用的检查习惯特别管用。第一点做完数值积分或方程求解后把结果代回原方程验算一下牛顿法就代回f(x)看是否接近0积分结果可以用不同的被积函数或区间去测比如sin(x)在[0, pi]上的积分理论值是2。第二点检查算法的收敛阶Romberg外推每层精度大约提升一个数量级如果你的结果没有这个趋势多半是实现有bug。第三点注意矩阵运算中的维度检查用size()打印几个关键矩阵的形状尤其是重新组织代码时转置操作经常会引入隐藏错误。这些习惯看着琐碎但数值分析这门课最后拿高分拼的往往就是这些细节——算法原理大家都懂差距就体现在谁能更快更稳定地跑出可信结果并且把过程和验证写清楚。本文还有配套的精品资源点击获取