★ FEATURED ARTICLE 简介这份资源是面向数值分析、科学计算与工程仿真学习者的Matlab算法实现合集覆盖常微分方程求解、线性方程组、特征值计算与曲线拟合等典型问题适合正在修读计算方法课程或需要快速复现经典算法的本科生、研究生及工程技术人员。压缩包共12个文件全部为m脚本整体约5KB按算法模块独立成文便于逐个调试与对照学习。内容涵盖欧拉法、高斯消元、豪斯霍尔德变换、线性最小二乘法、牛顿法、幂迭代法以及Runge-Kutta家族算法并配有returnA等辅助函数可帮助读者理解迭代格式、矩阵变换与误差控制的具体实现。目前已有687人学习下载代码结构简洁适合作为课程实验、算法验证与数值方法入门的参考脚本也可在此基础上扩展为更复杂的工程计算程序。1. 计算方法 Matlab全套代码从伪代码到能跑通的工程实现很多人第一次接触数值计算是在课本上看到一堆公式和伪代码考试能过但真到用 Matlab 落地时发现连一个能跑通的脚本都写不出来。这个标题指向的就是把这层窗户纸捅破——把计算方法里那些经典算法用 Matlab 写成可直接运行、可修改、可验证的代码。它解决的不是“懂不懂原理”的问题而是“能不能在半小时内把算法变成可调试的工程代码”的问题。适合两类人正在上计算方法课、需要交作业或做课程设计的学生以及工作中需要快速验证数值方案、不想从零推导的工程师。我见过太多人卡在“知道公式但写不出循环”这一步这套代码的价值就在于把中间那段空白填上。2. 先选对算法再写代码计算方法里最容易选错的三个场景2.1 方程求根二分法、牛顿法、割线法到底用哪个方程求根是计算方法里最基础也最容易被低估的部分。很多人拿到一个非线性方程第一反应是直接上牛顿法因为收敛快。但牛顿法需要求导如果导数表达式复杂或者根本写不出来就得换割线法。二分法虽然慢但它有一个牛顿法和割线法都没有的优点只要区间端点函数值异号就一定收敛不存在发散的风险。我一般会按这个顺序判断先看能不能解析求导能就用牛顿法不能求导但能拿到两个点用割线法如果连初始点都不好选先用二分法把区间缩到足够小再切牛顿法。这个组合策略在实际工程里比单用一种方法稳得多。下面是一个牛顿法的 Matlab 实现带迭代次数上限和收敛容差function [x, iter] newton_solve(f, df, x0, tol, max_iter) % f: 原函数句柄 % df: 导函数句柄 % x0: 初始猜测值 % tol: 收敛容差 % max_iter: 最大迭代次数 x x0; for iter 1:max_iter fx f(x); dfx df(x); if abs(dfx) 1e-14 error(导数接近零牛顿法失效请换初始点或改用其他方法); end x_new x - fx / dfx; if abs(x_new - x) tol x x_new; return; end x x_new; end warning(达到最大迭代次数可能未收敛); end这段代码的关键在abs(dfx) 1e-14这个判断。牛顿法的迭代公式是 ( x_{k1} x_k - f(x_k)/f(x_k) )当导数接近零时分母极小会导致步长爆炸直接跑飞。加这个保护能避免大部分“程序跑着跑着变成 NaN”的情况。tol一般取1e-8到1e-12取决于你对精度的要求。max_iter设 100 到 500 都行正常收敛的方程一般 10 次以内就到位了。2.2 线性方程组直接法还是迭代法看矩阵规模说话线性方程组 ( Ax b ) 的求解是计算方法的核心内容之一。Matlab 里最直接的是x A \ b但很多人不知道这个反斜杠背后会根据矩阵结构自动选择算法。对于稠密小矩阵它用 LU 分解对于稀疏矩阵它会尝试各种稀疏分解。问题是当矩阵规模到几千甚至上万时直接法的内存和计算量会急剧上升这时候就得考虑迭代法。判断标准很简单矩阵维度小于 1000 且稠密直接用A \ b矩阵维度大于 5000 或者稀疏优先考虑共轭梯度法或 GMRES。但迭代法有个前提矩阵最好是对称正定的或者至少对角占优否则收敛性没保证。下面是一个预处理共轭梯度法的实现框架function [x, flag, relres, iter] pcg_solve(A, b, M, tol, max_iter) % A: 系数矩阵对称正定 % b: 右端向量 % M: 预条件子可传空 % tol: 相对残差容差 % max_iter: 最大迭代次数 n length(b); x zeros(n, 1); r b - A * x; if isempty(M) z r; else z M \ r; end p z; rz r * z; bnorm norm(b); for iter 1:max_iter Ap A * p; alpha rz / (p * Ap); x x alpha * p; r r - alpha * Ap; relres norm(r) / bnorm; if relres tol flag 0; return; end if isempty(M) z r; else z M \ r; end rz_new r * z; beta rz_new / rz; p z beta * p; rz rz_new; end flag 1; end预条件子M的选择是 PCG 的灵魂。最简单的取法是M diag(diag(A))也就是对角预条件能改善条件数但不增加太多计算量。如果矩阵来自泊松方程离散化可以考虑不完全 Cholesky 分解。tol一般取1e-6到1e-10max_iter可以设为n或者2n。注意flag返回 1 表示未收敛这时候不要直接用结果得检查矩阵性质或者换预条件子。2.3 数值积分梯形法、辛普森法、高斯求积的精度与代价数值积分看起来简单但选错方法会导致精度差好几个数量级。梯形法是一阶精度辛普森法是三阶精度高斯求积在同样节点数下精度最高。但高斯求积的节点和权重需要查表或者用 Legendre 多项式生成实现复杂度高一些。我一般这样选如果被积函数光滑且能方便地增加节点用辛普森法就够了如果被积函数有端点奇异性考虑高斯求积或者变量替换如果只是快速估算梯形法加足够多的节点也能用。Matlab 内置的integral函数用的是自适应辛普森法大多数情况下直接调它就行但如果你需要自己控制节点分布就得手写。下面是一个自适应辛普森法的递归实现function I adaptive_simpson(f, a, b, tol, depth) % f: 被积函数句柄 % a, b: 积分区间 % tol: 容差 % depth: 当前递归深度 if nargin 5 depth 0; end max_depth 20; c (a b) / 2; h b - a; fa f(a); fb f(b); fc f(c); S1 h / 6 * (fa 4 * fc fb); d (a c) / 2; e (c b) / 2; fd f(d); fe f(e); S2 h / 12 * (fa 4 * fd 2 * fc 4 * fe fb); if depth max_depth I S2; return; end if abs(S2 - S1) 15 * tol I S2 (S2 - S1) / 15; else I adaptive_simpson(f, a, c, tol/2, depth1) ... adaptive_simpson(f, c, b, tol/2, depth1); end end这里的关键是abs(S2 - S1) 15 * tol这个判断。辛普森法的误差估计是 ( |S_2 - S_1|/15 )所以当这个值小于容差时用 Richardson 外推得到更精确的结果。max_depth设 20 是防止被积函数有剧烈振荡时递归过深导致栈溢出。tol一般取1e-6到1e-10取决于被积函数的尺度。3. 把公式变成代码三个核心模块的 Matlab 实现细节3.1 矩阵分解LU、QR、Cholesky 的代码结构与调用方式矩阵分解是线性方程组、最小二乘、特征值问题的底层支撑。Matlab 内置了lu、qr、chol等函数但理解它们的输出结构和适用条件比直接调用更重要。LU 分解把 ( A ) 拆成 ( PA LU )其中 ( P ) 是置换矩阵( L ) 是下三角( U ) 是上三角。Matlab 的[L, U, P] lu(A)返回的就是这三个矩阵。注意P是置换矩阵满足 ( A P * L * U )。如果你只需要解方程直接用A \ b就行Matlab 会自动做 LU 分解并前代后回代。QR 分解把 ( A ) 拆成 ( A QR )其中 ( Q ) 是正交矩阵( R ) 是上三角。Matlab 的[Q, R] qr(A)返回的Q是正交矩阵R是上三角。QR 分解在最小二乘问题里特别有用因为正交变换不改变向量范数所以 ( |Ax - b| |QRx - b| |Rx - Qb| )这样最小二乘问题就变成了一个上三角方程组的求解。Cholesky 分解只适用于对称正定矩阵把 ( A ) 拆成 ( A LL )其中 ( L ) 是下三角。Matlab 的L chol(A, lower)返回下三角的L。Cholesky 分解的计算量大约是 LU 分解的一半而且不需要选主元数值稳定性好。但前提是矩阵必须对称正定如果矩阵只是对称但不正定chol会报错。下面是一个用 Cholesky 分解解线性方程组的完整示例% 构造一个对称正定矩阵 n 100; A gallery(poisson, sqrt(n)); % 泊松矩阵对称正定 b ones(size(A, 1), 1); % Cholesky 分解 L chol(A, lower); % 前代解 Ly b y L \ b; % 回代解 Lx y x L \ y; % 验证残差 residual norm(A * x - b) / norm(b); fprintf(相对残差: %e\n, residual);这段代码里gallery(poisson, sqrt(n))生成的是二维泊松方程离散化后的矩阵天然对称正定适合演示 Cholesky 分解。L \ b是前代L \ y是回代Matlab 会自动识别三角矩阵并选择对应的求解算法。residual用来验证解的精度正常情况下应该在1e-14量级。3.2 特征值与奇异值幂法、反幂法与 Matlab 内置函数的配合特征值问题在振动分析、主成分分析、稳定性判断里到处都是。Matlab 的eig函数能算全部特征值但对于大规模稀疏矩阵算全部特征值既慢又占内存。这时候幂法和反幂法就派上用场了。幂法用来求模最大的特征值及其对应的特征向量。基本思想是反复用矩阵乘向量向量方向会逐渐对齐到主特征向量。反幂法用来求模最小的特征值做法是对 ( A - \sigma I ) 做幂法其中 ( \sigma ) 是位移量。如果 ( \sigma ) 取在某个特征值附近反幂法会快速收敛到那个特征值对应的特征向量。下面是一个带位移的反幂法实现function [lambda, v, iter] inverse_power(A, sigma, tol, max_iter) % A: 矩阵 % sigma: 位移量 % tol: 收敛容差 % max_iter: 最大迭代次数 n size(A, 1); v rand(n, 1); v v / norm(v); for iter 1:max_iter % 解 (A - sigma*I) * w v w (A - sigma * eye(n)) \ v; lambda v * w; % Rayleigh 商 v_new w / norm(w); if norm(v_new - v) tol v v_new; lambda sigma 1 / lambda; return; end v v_new; end warning(达到最大迭代次数); end这里lambda v * w是 Rayleigh 商用来估计特征值。lambda sigma 1 / lambda是位移反幂法的特征值修正公式。sigma的选择很关键如果知道特征值的大致范围取靠近目标特征值的值能加速收敛如果完全不知道可以先跑一次eig看看特征值分布再选sigma。3.3 常微分方程初值问题欧拉法、RK4 与 ode45 的适用边界常微分方程初值问题 ( y f(t, y), y(t_0) y_0 ) 的数值解法分两大阵营单步法和多步法。欧拉法是最简单的单步法一阶精度基本只用来教学。RK4 是四阶精度的单步法稳定性和精度都不错是工程计算里的主力。Matlab 的ode45用的是自适应步长的 Runge-Kutta-Fehlberg 方法精度更高但每一步的计算量也更大。选择标准如果方程不 stiff 且对精度要求一般RK4 手写就够了如果方程 stiff 或者精度要求高用ode45或ode15s如果只是快速验证欧拉法加小步长也能凑合但不推荐。下面是一个 RK4 的实现function [t, y] rk4_solve(f, tspan, y0, h) % f: 右端函数句柄f(t, y) % tspan: [t0, tf] % y0: 初始值 % h: 步长 t0 tspan(1); tf tspan(2); t t0:h:tf; n length(t); y zeros(length(y0), n); y(:, 1) y0; for i 1:n-1 k1 f(t(i), y(:, i)); k2 f(t(i) h/2, y(:, i) h/2 * k1); k3 f(t(i) h/2, y(:, i) h/2 * k2); k4 f(t(i) h, y(:, i) h * k3); y(:, i1) y(:, i) h/6 * (k1 2*k2 2*k3 k4); end endh的选择直接影响精度和稳定性。对于 RK4步长太大精度不够步长太小计算量浪费。一般先用ode45跑一遍看看它自动选的步长范围再决定手写 RK4 时h取多少。如果方程有快速变化的区域h要相应减小。4. 避坑与排查计算方法 Matlab 代码里最常见的五个翻车点4.1 矩阵接近奇异但没报错结果全是垃圾现象A \ b能跑通但解出来的x数量级离谱或者残差norm(A*x - b)很大。原因矩阵条件数太大接近奇异。Matlab 的\不会因为矩阵接近奇异就报错它只是给出一个数值解但这个解可能完全没有意义。解决先算条件数cond(A)如果大于1e12说明矩阵病态。这时候要么用正则化比如 Tikhonov 正则化要么检查矩阵构造过程是不是有问题。另一个办法是用pinv(A) * b求伪逆解但伪逆的计算量更大只适合小矩阵。4.2 迭代法不收敛残差震荡或发散现象PCG 或 GMRES 跑了很多步残差不但不降反而上升。原因矩阵不是对称正定或者预条件子选得不好。PCG 要求矩阵对称正定如果矩阵只是对称但不正定PCG 可能不收敛。GMRES 对矩阵性质要求低一些但内存开销随迭代步数线性增长。解决先检查矩阵是否对称正定用issymmetric(A)和all(eig(A) 0)验证。如果不是换 GMRES 或者 BiCGSTAB。预条件子可以试试M diag(diag(A))或者不完全 LU 分解[L, U] ilu(A)然后用M L * U。4.3 数值积分遇到端点奇异性结果偏差巨大现象被积函数在端点附近趋于无穷梯形法或辛普森法算出来的积分值明显偏小。原因端点奇异性导致被积函数在端点附近变化剧烈固定节点的方法无法准确捕捉。解决做变量替换把奇异性消掉。比如 ( \int_0^1 x^{-1/2} f(x) dx ) 可以令 ( x t^2 )变成 ( \int_0^1 2 f(t^2) dt )奇异性就没了。或者用高斯求积因为高斯节点在端点附近分布更密。4.4 RK4 步长选太大结果完全失真现象RK4 跑出来的解和ode45差很远或者解在某个区域突然爆炸。原因步长h超过了稳定性区域。RK4 的稳定性区域是复平面上的一个有限区域如果方程的特征值落在区域外数值解会发散。解决先用ode45跑一遍看它自动选的步长范围。如果ode45的步长在1e-4量级手写 RK4 时h不要超过1e-3。另一个办法是自适应步长根据局部误差估计调整h但实现复杂度高一些。4.5 函数句柄传参错误导致维度不匹配现象调用f(t, y)时报错说矩阵维度不一致。原因Matlab 的函数句柄对参数个数和维度很敏感。如果f定义时是(t, y) ...但y传进去是行向量而函数内部按列向量处理就会出错。解决在函数内部统一用列向量或者在调用前用y(:)强制转成列向量。另一个常见错误是f定义时用了额外的参数但没通过匿名函数捕获比如f (t, y) my_func(t, y, a, b)这里的a和b必须在定义f之前就有值。5. 从能跑到好用三个让代码更稳的进阶习惯5.1 用残差和条件数做后验验证而不是只看输出很多人跑完代码看到x有值就以为对了这是最大的坑。我一般会在每个求解函数里加一段后验验证算残差norm(A*x - b) / norm(b)算条件数cond(A)如果残差大于1e-6或者条件数大于1e12就打印警告。这个习惯能帮你提前发现大部分数值问题。function x solve_and_verify(A, b) x A \ b; res norm(A * x - b) / norm(b); c cond(A); if res 1e-6 warning(相对残差偏大: %e解可能不可靠, res); end if c 1e12 warning(矩阵条件数过大: %e结果对扰动敏感, c); end end5.2 把算法参数做成可配置结构体而不是硬编码硬编码参数是代码复用的最大障碍。我习惯把容差、最大迭代次数、步长这些参数打包成一个结构体调用时传进去。这样换一个问题只需要改结构体不用动函数体。opts.tol 1e-10; opts.max_iter 500; opts.h 1e-3; [x, iter] newton_solve(f, df, x0, opts.tol, opts.max_iter);5.3 用单元测试验证边界情况而不是只跑一个例子数值代码的边界情况特别多零矩阵、单位矩阵、对角矩阵、接近奇异的矩阵。我一般会写几个小测试覆盖这些情况确保代码不会在边界上翻车。% 测试零矩阵 A zeros(3); b ones(3, 1); try x A \ b; warning(零矩阵应该报错或给出警告); catch disp(零矩阵正确报错); end % 测试单位矩阵 A eye(3); b [1; 2; 3]; x A \ b; assert(norm(x - b) 1e-14, 单位矩阵求解错误);这些习惯看起来简单但能帮你省下大量调试时间。我自己的经验是数值代码的 bug 往往不在算法本身而在参数传递、维度匹配和边界处理上。把这三块守住代码的可靠性会提升一个档次。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
