简介一份面向需要运用MATLAB解决高等应用数学问题的学习者与科研人员的源码包内含318个可直接运行的.m程序覆盖微积分与微分方程、线性代数与矩阵运算、数值优化、统计与概率、信号处理、图像处理、控制理论、动态系统模拟、智能优化算法及并行计算等多个应用方向。压缩包共319个文件以313个.m源程序为主另有4个.asv自动备份、1个.mdl系统模型和1个.pdf说明文档整体大小48.23MB。目前已有576人学习代码按典型问题组织便于按需查阅。通过研读这些实例读者可以掌握MATLAB在符号计算、数值求解、算法设计与程序实现方面的核心技巧理解数学理论与编程实践的对应关系并能够将相关方法迁移到自己的研究或工程项目中。1. 318个源程序包的背后高等应用数学的 MATLAB 解法地图解压高等应用数学问题的matlab求解318个源程序.zip之后最先看到的是一串按examp 章号 例题号命名的脚本examp2_1.m、examp3_27.m、examp4_38.m其中还夹杂着examp8_20.asv、examp6_20.asv这类编辑器自动保存产生的备份文件。这套资源的真正价值不在单个脚本的技巧而在于它按教材章节把高等应用数学中最常遇到的问题——微分方程、线性代数、数值优化、统计建模——逐一落到可运行的 MATLAB 代码上。对正在做课程设计、写论文或者处理工程计算任务的人来说找到对应章节的脚本改参数、换边界条件几分钟内就能看到中间结果远比从零推导算法再实现一遍高效。318 个源程序也提供了从入门到进阶的完整梯度新手看语法和调用方式熟手看一个数学问题从建模到代码的组织顺序。2. 微分方程求解的双路线dsolve 符号解与 ode45、ode15s 数值解2.1 先用 dsolve 探路符号解建立结构认知微分方程是高等应用数学中出现频率最高的板块之一源程序里examp6_20、examp7_9这类编号集中出现对应的就是常微分方程与差分方程章节。拿到一个具体的 ODE我的习惯是先试试符号求解即使最终必然要走数值路线符号解也能在一开始就告诉你解的形态衰减、振荡、发散还是趋于常数。syms y(t) % 声明符号变量 y(t) eqn diff(y, t) -2*y t; % 定义一阶常微分方程 cond y(0) 1; % 初值条件 ySol(t) dsolve(eqn, cond) % 返回符号解函数逻辑说明diff(y, t)是对符号函数 y 求关于 t 的一阶导eqn是一个方程对象在这里用于定义方程而不是比较。dsolve的输入可以只传方程也可以连同初值条件一起传输出是一个符号函数可以直接参与fplot、diff、int等后续符号操作。这组命令解决的是“拿到微分方程先判断它长什么样”的问题而不是直接跳到数值计算。参数说明二阶方程需要两个初值条件写法是cond [y(0) 1, Dy(0) 0]其中Dy表示一阶导数D2y表示二阶导数这是符号工具箱的固定记号。若返回结果里出现RootOf或隐式积分项说明解析解存在但无法用初等函数显式表达这时应当换数值路线。提示判断符号解是否可用最稳的方法是把ySol代回原方程并化简simplify(lhs(eqn) - rhs(eqn))结果应为 0。2.2 数值解主路线ode45 的自适应步长与容差控制当方程无法解析求解或者解析解形式过于复杂ode45就是默认选择。它基于 Dormand-Prince 显式 Runge-Kutta 对用 4 阶公式推进、5 阶公式估计误差步长由误差控制自动加密或放宽因此不需要用户手工指定步长。tspan [0 5]; % 积分区间 y0 1; % 初值 [tn, yn] ode45((t, y) -2*y t, tspan, y0); plot(tn, yn, o); hold on; fplot(ySol, [0 5], r--); legend(ode45 数值解, dsolve 解析解);逻辑说明ode45第一个参数是函数句柄格式固定为接受(t, y)并返回一阶导数值第二个参数tspan是积分区间也可以传0:0.1:5这种节点向量强制输出在指定节点上第三个参数y0是初值多变量系统时是列向量。返回的tn是自适应选出的节点序列yn是对应解两者长度一致。把数值解和解析解画在同一张图上是检查实现正确性的第一道工序两条线明显分离一般意味着方程写错或初值不对。参数说明用户不直接指定步长而是通过odeset调整容差间接控制。opts odeset(RelTol, 1e-6, AbsTol, 1e-8)是常用起点。RelTol控制相对误差AbsTol控制绝对误差当解接近 0 时相对误差失去意义此时绝对误差才是决定精度的因素。把RelTol从1e-6压到1e-9通常会显著增加计算量但对工程结果精度提升有限除非在做基准测试否则没必要追求极端容差。2.3 刚性系统与 ode15s 的切换判据ode45并非万能。当方程的特征值跨越多个数量级时显式方法为了保证稳定会把步长压到极小表现为求解时间异常长或者在命令窗口持续输出 “Iteration is not making good progress” 之类的警告。特征值差距大的系统被称为刚性系统这种情况下应该切换到ode15s。opts odeset(RelTol, 1e-6, AbsTol, 1e-8); f_rigid (t, y) -1000*(y - cos(t)) - sin(t); [t, y] ode15s(f_rigid, [0 10], 2, opts);逻辑说明这个方程的线性部分系数达到 1000解在短时间内经历快速暂态之后进入慢变阶段正是典型的刚性结构。ode15s使用多步变阶后向差分公式隐式迭代对刚性问题的稳定性远好于ode45。在 tspan 相同、初值相同的情况下ode45可能以万级步数逼近ode15s只需要几十步。参数说明同前容差控制是核心。ode15s在同容差下往往比ode45耗费更多单步计算但在刚性问题上总步数少几个数量级整体仍然快得多。判断是否刚性的经验法则是ode45跑不动、但把tspan缩短来看解本身很光滑这种矛盾就是刚性信号。求解器适用问题典型特征失败表现ode45非刚性显式 RK4/5 阶步数过多、长时间无结果ode23低精度非刚性显式 RK2/3 阶容差苛刻时效率低ode15s刚性/隐式多步变阶 BDF非刚性问题精度偏低ode113光滑变阶 Adams精度高、光滑问题解不光滑时失效2.4 把被积函数拆成独立文件而不是堆匿名函数318 个源程序里有不少脚本采用“主脚本 子函数”的组织方式这个习惯值得保留。实际项目里被积函数往往需要调用参数、查表或做中间计算塞在(t,y)匿名函数里可读性和可维护性都会变差。常见做法是把方程写成独立的 function 文件function dydt myODE(t, y) % t: 自变量方程与 t 无关时也保留占位 % y: 状态向量 dydt -2*y t; end主脚本里调用[tn, yn] ode45(myODE, [0 5], 1);。如果参数需要从外部传入用嵌套函数捕获或者用(t,y) myODE(t, y, k)再包一层都可以。这样后续做参数扫描、批量初值实验时只需要在外层循环里更新参数方程文件本身不用动。源程序里那些按章节编号的脚本很多就是这样组织的先读主脚本再跳进被调用的函数能很快看清一个数学问题被拆分成了哪几个计算环节。3. 线性代数核心算法LU、QR、Cholesky 分解的选型与验证3.1 矩阵求逆不是第一选项分解才是工程计算里经常看到x inv(A) * b但 MATLAB 的官方实现和绝大多数数值计算库在解线性方程组时都不会先求逆。原因是A\b底层按矩阵结构自动选择算法一般方阵走 LU 分解加回代对称正定走 Cholesky三角矩阵走前代/回代超定方程走 QR。显式求逆既浪费计算量数值稳定性也差同样的矩阵用inv和用分解得到的解残差可能差好几个数量级。源程序examp3_27、examp4_38这类编号覆盖的就是线性方程组与矩阵分解章节。我一般会在代码里把几种分解都显式写出来做交叉验证而不是只依赖\的自动调度A gallery(lehmer, 8); % 对称正定 Lehmer 矩阵条件数约 1e7 量级 b ones(8, 1); x1 A \ b; % 自动选路对称正定走 Cholesky [L, U, P] lu(A); % LU 分解P 是置换矩阵 x2 U \ (L \ (P * b)); R chol(A); % Cholesky 分解A R*R x3 R \ (R \ b); % 相对残差对比三种方式应接近机器精度 res [norm(A*x1 - b)/norm(b), ... norm(A*x2 - b)/norm(b), ... norm(A*x3 - b)/norm(b)]; disp(res);逻辑说明lu返回的P是行置换矩阵所以回代顺序是P*b不能直接写L\(U\b)否则结果错得离谱。chol要求矩阵对称正定gallery(lehmer, 8)返回的 Lehmer 矩阵满足这个条件但其特征值相差很大、条件数高非常适合用来暴露数值问题。三种方式求出的x1、x2、x3应当高度一致相对残差在1e-14附近属于正常。参数说明A\b不只是语法糖它背后是完整的算法调度。理解这一点很重要因为同样的\遇到非方阵时会自动切换语义超定方程A (mn)得到最小二乘解欠定方程得到基础解系中的一个解。所以看到别人的代码里大量使用\时不要急着改成显式分解先弄清楚矩阵的形状和需求。分解方法适用条件主要用途MATLAB函数LU一般方阵线性方程组、行列式计算luQR长方形矩阵、超定方程最小二乘、特征值迭代qrCholesky对称正定快速求解、协方差分解chol3.2 特征值问题全量 eig 与稀疏 eigs 的适用边界特征值计算是矩阵分析的另一个高频主题PCA、模态分析、图拉普拉斯谱聚类本质都是特征值分解问题。小规模矩阵直接用eig返回全部特征值当矩阵维度上万且稀疏只关心少数几个模最大或接近某个目标值的特征值时要换eigs。d eig(A); % 全部特征值8 阶矩阵直接算 opts.tol 1e-6; % 迭代法收敛容差 [V, D] eigs(A, 6, largestabs, opts); % 求模最大的 6 个特征值 res_eig norm(A*V - V*D, fro) / norm(A, fro);逻辑说明eigs返回的特征向量按列排列在V中D是对角阵A*V - V*D就是残差矩阵。这里用 Frobenius 范数做归一化得到相对残差。这个残差如果超过1e-6通常意味着容差设置过松或矩阵条件数过差需要把opts.tol再压一挡或者检查矩阵是否真的稀疏。参数说明largestabs表示求模最大的特征值对应旧版本里的lm。还可以传smallestabs求模最小、largestreal求实部最大等。对于广义特征值问题Ax λBx直接调用eigs(A, B, k, sigma)即可这在结构动力学和模态分析脚本里很常见。3.3 条件数与解的可信度解矩阵方程、算特征值之后我习惯性检查一个数条件数cond(A)。它衡量方程组解对输入扰动的放大倍数量级接近1e7时解的有效数字大约比机器精度少 7 位双精度下只剩约 9 位有效数字。c cond(A); fprintf(cond(A) %.3e\n, c);逻辑说明cond(A)计算的是最大奇异值与最小奇异值之比是判断矩阵是否病态的最直观指标。配合前面的相对残差能区分两种常见情况残差大但条件数正常说明代码有 bug残差小但条件数极大说明结果本身可信度有限应该考虑预处理、换基或改用精度更高的问题表达方式而不是继续调整求解参数。工程上常见的误用是只看残差不看条件数。残差小只代表A*x接近b并不代表x接近真实解条件数大时一个微小的扰动就能让解发生巨大变化。所以验证一套线性代数求解代码cond和残差要同时看二者缺一不可。4. 数值优化与统计建模从 fmincon 到 ga 与 kmeans4.1 把实际问题规约成标准优化形式源程序进入到examp9_25、examp10_32这类编号话题就转到数值优化。用 MATLAB 优化工具箱做优化第一步不是调函数而是把问题改写成标准形式目标函数f(x)最小化约束统一为线性不等式A*x b、线性等式Aeq*x beq、变量上下界lb/ub、非线性约束c(x)0和ceq(x)0四类。这样做的好处是同一套问题定义可以直接在fmincon、ga、patternsearch之间切换换算法时不需要改约束建模逻辑。以下面的问题为例目标函数是凸二次型约束是圆形区域加一个半平面fun (x) (x(1)-2)^2 (x(2)-1)^2; % 目标函数 x0 [0; 0]; % 初始点 A [-1 -1]; b -1; % 线性不等式x1 x2 1 lb [-10; -10]; ub [10; 10]; % 变量上下界 nonlcon (x) deal(x(1)^2 x(2)^2 - 4, []); % 非线性不等式 0 xopt fmincon(fun, x0, A, b, [], [], lb, ub, nonlcon);逻辑说明fmincon的完整签名是fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon)没有等式约束所以Aeq、beq传空数组。nonlcon返回两个输出第一个是不等式约束的值要求小于等于 0第二个是等式约束要求等于 0用deal可以一条语句完成返回这是源程序里比较常见的写法。参数说明x0是局部优化算法的起点。凸问题从哪出发结果一致非凸问题则几乎由初值决定结局。建议至少取 5 到 10 个随机起点分别求解对比目标函数值而不是只看一次结果。我的做法是用for循环配合fmincon批量跑再取最优值这个方法在源程序的基础上加几行就能实现。4.2 fmincon 的内置算法与关键公差fmincon在近几个主版本里默认算法是 interior-point内点法它对中小规模约束优化问题表现稳定求解器按约束类型自动选择子算法。常用设置是把显示选项和容差显式写出来方便排查收敛过程opts optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... OptimalityTolerance, 1e-8, ... ConstraintTolerance, 1e-8, ... MaxIterations, 1000); xopt fmincon(fun, x0, A, b, [], [], lb, ub, nonlcon, opts);参数说明OptimalityTolerance控制一阶最优性条件的容差太小会让求解器在最优解附近反复迭代ConstraintTolerance控制约束违反程度工程问题1e-6已经足够实验室精度追求才压到1e-10以下。Display设为iter能在命令窗口实时观察迭代步长和目标值变化排查不收敛问题时比只看最终结果有效得多。如果目标函数不可导内点法会报梯度相关错误这时可以改用Algorithm, sqp它对目标函数的光滑性要求更低。但更好的策略是直接换到全局搜索类的无导数算法也就是下一节的内容。4.3 全局搜索兜底ga 遗传算法的参数经验当目标存在多个局部极值fmincon的多起点策略仍然可能漏掉全局最优这时我会上ga。遗传算法不需要梯度对目标函数形态几乎无要求代价是收敛慢且结果有随机性。这个取舍在高等应用数学问题里很常见能接受计算时间换来的是更稳的全局搜索能力。opts optimoptions(ga, ... PopulationSize, 100, ... MaxGenerations, 200, ... Display, iter); [xg, fval] ga(fun, 2, A, b, [], [], lb, ub, nonlcon, opts);逻辑说明ga的第二个参数是变量个数这里fun是二维问题所以传 2。约束建模方式与fmincon完全一致这也是前面强调标准形式的原因。ga不需要初值x0它用随机种群启动所以每次运行结果会有差异。参数说明PopulationSize默认 50非线性约束强时可以升到 100 到 200MaxGenerations默认是变量数的 100 倍这里二维问题就是 200 代。Displayiter可以观察每次迭代的最优值和平均适应度如果连续多代最优值不变说明已经收敛。正式使用前用rng(42)固定随机种子方便复现和对比调参。求解器适用场景是否需梯度随机性关键参数fminunc无约束光滑问题需要无OptimalityTolerancefmincon有约束光滑问题需要无ConstraintTolerancega非光滑/离散/多极值不需要有PopulationSize, MaxGenerationspatternsearch非光滑小规模不需要无MeshSize, PollMethod4.4 回归、概率分布与 kmeans 聚类的配合优化之外源程序中也有大量统计相关内容。线性回归用fitlm一步到位返回系数、p 值和置信区间概率分布拟合用fitdist聚类问题默认kmeans。这三者经常串联使用先做回归理解变量关系再对残差做分布假设检验最后对样本聚类找分组结构。rng(42); data [randn(100, 2) [2 2]; randn(100, 2) - [2 2]]; [idx, C] kmeans(data, 2, Replicates, 5); gscatter(data(:, 1), data(:, 2), idx); pd fitdist(data(:, 1), Normal); % 单变量正态分布拟合 mdl fitlm(data(:, 1), data(:, 2)); % 一元线性回归逻辑说明kmeans是基于距离的划分迭代算法对初始中心敏感。Replicates参数让 MATLAB 从多个随机初始中心出发返回总组内距离最小的解。fitdist返回的概率分布对象可以直接调用pdf、cdf和paramci等方法适合后续做置信区间分析。fitlm的输出里包含系数估计、标准误和 F 统计量是回归分析的统一入口。参数说明如果数据量纲差异大比如一列是上万量级、另一列是 0 到 1 之间务必先标准化否则聚类和回归结果会被量纲大的特征主导。标准化用zscore即可这一步看起来简单实际效果非常关键。数据再复杂一些需要上模型时Deep Learning Toolbox 里的trainNetwork也能接手回归或分类任务但特征量纲处理的原则完全一致这套源程序里体现的建模思路可以直接迁移。5. .asv 文件恢复与运行时间回归源码包调试的实用习惯5.1 .asv 备份文件不是垃圾318 个源程序里混着examp8_20.asv、examp6_20.asv这样的文件这是 MATLAB 编辑器自动保存产生的副本。原.m文件被误删或改坏时把.asv改名为.m就能恢复到最近一次自动保存的状态。编辑器默认每隔几分钟自动保存一份.asv所以备份内容通常只比你最后一次保存晚几分钟。当目录里混了一堆.asv需要统一处理时批量重命名即可。Linux/macOS 终端里执行for f in *.asv; do mv $f ${f%.asv}.m; doneWindows 的 PowerShell 里对应写法是Get-ChildItem *.asv | Rename-Item -NewName { $_.Name -replace \.asv$, .m }。需要注意若同名.m已存在先对比文件修改时间再决定覆盖方向别拿旧备份覆盖新代码。5.2 数值结果的残差验证与 assert 断言我习惯在脚本尾部写断言用机器判断结果是否可信而不是靠肉眼盯数值。对前面的 ODE 例子yRef double(arrayfun((t) ySol(t), tn)); assert(norm(yn - yRef, inf) 1e-4, 数值解与解析解偏差过大);逻辑说明arrayfun把符号解ySol映射到数值节点tn上double完成符号到双精度转换norm(..., inf)取最大偏差。assert第一个条件为假时直接抛异常并输出提示这样批量跑源程序时万一某个参数改动让结果发散能立刻定位到具体脚本而不是一路跑到底再回来肉眼排查。5.3 用 timeit 替代 tic/toc 计时在优化算法或大规模矩阵计算里算法改进效果不能靠一次tic/toc判断因为它受系统负载影响大。timeit会多次运行函数句柄并取稳健统计值更适合做性能回归测试f1 () ode45((t, y) -2*y t, [0 5], 1); t1 timeit(f1);参数说明timeit的输入必须是无参数或参数已捕获的函数句柄运行次数由它内部决定一般会在几百到几千次之间权衡精度与耗时。对比两个算法时分别timeit并重复三次以上再取均值结论比单次tic/toc可靠得多。如果需要更大的性能收益把热点代码写成 C 再编译成 mex 文件是 MATLAB 里运行外部语言的标准做法能把循环密集的数学内核提速数倍。另外注意一个数值显示陷阱1e100在命令行直接显示为1.0000e100继续参与某些计算时容易溢出成Inf。用sym(1e100)包一层即可看到完整精确值避免在调试大数问题时被浮点显示误导。本文还有配套的精品资源点击获取 MATLAB高等应用数学实战：微分方程、矩阵分解与优化算法详解
简介一份面向需要运用MATLAB解决高等应用数学问题的学习者与科研人员的源码包内含318个可直接运行的.m程序覆盖微积分与微分方程、线性代数与矩阵运算、数值优化、统计与概率、信号处理、图像处理、控制理论、动态系统模拟、智能优化算法及并行计算等多个应用方向。压缩包共319个文件以313个.m源程序为主另有4个.asv自动备份、1个.mdl系统模型和1个.pdf说明文档整体大小48.23MB。目前已有576人学习代码按典型问题组织便于按需查阅。通过研读这些实例读者可以掌握MATLAB在符号计算、数值求解、算法设计与程序实现方面的核心技巧理解数学理论与编程实践的对应关系并能够将相关方法迁移到自己的研究或工程项目中。1. 318个源程序包的背后高等应用数学的 MATLAB 解法地图解压高等应用数学问题的matlab求解318个源程序.zip之后最先看到的是一串按examp 章号 例题号命名的脚本examp2_1.m、examp3_27.m、examp4_38.m其中还夹杂着examp8_20.asv、examp6_20.asv这类编辑器自动保存产生的备份文件。这套资源的真正价值不在单个脚本的技巧而在于它按教材章节把高等应用数学中最常遇到的问题——微分方程、线性代数、数值优化、统计建模——逐一落到可运行的 MATLAB 代码上。对正在做课程设计、写论文或者处理工程计算任务的人来说找到对应章节的脚本改参数、换边界条件几分钟内就能看到中间结果远比从零推导算法再实现一遍高效。318 个源程序也提供了从入门到进阶的完整梯度新手看语法和调用方式熟手看一个数学问题从建模到代码的组织顺序。2. 微分方程求解的双路线dsolve 符号解与 ode45、ode15s 数值解2.1 先用 dsolve 探路符号解建立结构认知微分方程是高等应用数学中出现频率最高的板块之一源程序里examp6_20、examp7_9这类编号集中出现对应的就是常微分方程与差分方程章节。拿到一个具体的 ODE我的习惯是先试试符号求解即使最终必然要走数值路线符号解也能在一开始就告诉你解的形态衰减、振荡、发散还是趋于常数。syms y(t) % 声明符号变量 y(t) eqn diff(y, t) -2*y t; % 定义一阶常微分方程 cond y(0) 1; % 初值条件 ySol(t) dsolve(eqn, cond) % 返回符号解函数逻辑说明diff(y, t)是对符号函数 y 求关于 t 的一阶导eqn是一个方程对象在这里用于定义方程而不是比较。dsolve的输入可以只传方程也可以连同初值条件一起传输出是一个符号函数可以直接参与fplot、diff、int等后续符号操作。这组命令解决的是“拿到微分方程先判断它长什么样”的问题而不是直接跳到数值计算。参数说明二阶方程需要两个初值条件写法是cond [y(0) 1, Dy(0) 0]其中Dy表示一阶导数D2y表示二阶导数这是符号工具箱的固定记号。若返回结果里出现RootOf或隐式积分项说明解析解存在但无法用初等函数显式表达这时应当换数值路线。提示判断符号解是否可用最稳的方法是把ySol代回原方程并化简simplify(lhs(eqn) - rhs(eqn))结果应为 0。2.2 数值解主路线ode45 的自适应步长与容差控制当方程无法解析求解或者解析解形式过于复杂ode45就是默认选择。它基于 Dormand-Prince 显式 Runge-Kutta 对用 4 阶公式推进、5 阶公式估计误差步长由误差控制自动加密或放宽因此不需要用户手工指定步长。tspan [0 5]; % 积分区间 y0 1; % 初值 [tn, yn] ode45((t, y) -2*y t, tspan, y0); plot(tn, yn, o); hold on; fplot(ySol, [0 5], r--); legend(ode45 数值解, dsolve 解析解);逻辑说明ode45第一个参数是函数句柄格式固定为接受(t, y)并返回一阶导数值第二个参数tspan是积分区间也可以传0:0.1:5这种节点向量强制输出在指定节点上第三个参数y0是初值多变量系统时是列向量。返回的tn是自适应选出的节点序列yn是对应解两者长度一致。把数值解和解析解画在同一张图上是检查实现正确性的第一道工序两条线明显分离一般意味着方程写错或初值不对。参数说明用户不直接指定步长而是通过odeset调整容差间接控制。opts odeset(RelTol, 1e-6, AbsTol, 1e-8)是常用起点。RelTol控制相对误差AbsTol控制绝对误差当解接近 0 时相对误差失去意义此时绝对误差才是决定精度的因素。把RelTol从1e-6压到1e-9通常会显著增加计算量但对工程结果精度提升有限除非在做基准测试否则没必要追求极端容差。2.3 刚性系统与 ode15s 的切换判据ode45并非万能。当方程的特征值跨越多个数量级时显式方法为了保证稳定会把步长压到极小表现为求解时间异常长或者在命令窗口持续输出 “Iteration is not making good progress” 之类的警告。特征值差距大的系统被称为刚性系统这种情况下应该切换到ode15s。opts odeset(RelTol, 1e-6, AbsTol, 1e-8); f_rigid (t, y) -1000*(y - cos(t)) - sin(t); [t, y] ode15s(f_rigid, [0 10], 2, opts);逻辑说明这个方程的线性部分系数达到 1000解在短时间内经历快速暂态之后进入慢变阶段正是典型的刚性结构。ode15s使用多步变阶后向差分公式隐式迭代对刚性问题的稳定性远好于ode45。在 tspan 相同、初值相同的情况下ode45可能以万级步数逼近ode15s只需要几十步。参数说明同前容差控制是核心。ode15s在同容差下往往比ode45耗费更多单步计算但在刚性问题上总步数少几个数量级整体仍然快得多。判断是否刚性的经验法则是ode45跑不动、但把tspan缩短来看解本身很光滑这种矛盾就是刚性信号。求解器适用问题典型特征失败表现ode45非刚性显式 RK4/5 阶步数过多、长时间无结果ode23低精度非刚性显式 RK2/3 阶容差苛刻时效率低ode15s刚性/隐式多步变阶 BDF非刚性问题精度偏低ode113光滑变阶 Adams精度高、光滑问题解不光滑时失效2.4 把被积函数拆成独立文件而不是堆匿名函数318 个源程序里有不少脚本采用“主脚本 子函数”的组织方式这个习惯值得保留。实际项目里被积函数往往需要调用参数、查表或做中间计算塞在(t,y)匿名函数里可读性和可维护性都会变差。常见做法是把方程写成独立的 function 文件function dydt myODE(t, y) % t: 自变量方程与 t 无关时也保留占位 % y: 状态向量 dydt -2*y t; end主脚本里调用[tn, yn] ode45(myODE, [0 5], 1);。如果参数需要从外部传入用嵌套函数捕获或者用(t,y) myODE(t, y, k)再包一层都可以。这样后续做参数扫描、批量初值实验时只需要在外层循环里更新参数方程文件本身不用动。源程序里那些按章节编号的脚本很多就是这样组织的先读主脚本再跳进被调用的函数能很快看清一个数学问题被拆分成了哪几个计算环节。3. 线性代数核心算法LU、QR、Cholesky 分解的选型与验证3.1 矩阵求逆不是第一选项分解才是工程计算里经常看到x inv(A) * b但 MATLAB 的官方实现和绝大多数数值计算库在解线性方程组时都不会先求逆。原因是A\b底层按矩阵结构自动选择算法一般方阵走 LU 分解加回代对称正定走 Cholesky三角矩阵走前代/回代超定方程走 QR。显式求逆既浪费计算量数值稳定性也差同样的矩阵用inv和用分解得到的解残差可能差好几个数量级。源程序examp3_27、examp4_38这类编号覆盖的就是线性方程组与矩阵分解章节。我一般会在代码里把几种分解都显式写出来做交叉验证而不是只依赖\的自动调度A gallery(lehmer, 8); % 对称正定 Lehmer 矩阵条件数约 1e7 量级 b ones(8, 1); x1 A \ b; % 自动选路对称正定走 Cholesky [L, U, P] lu(A); % LU 分解P 是置换矩阵 x2 U \ (L \ (P * b)); R chol(A); % Cholesky 分解A R*R x3 R \ (R \ b); % 相对残差对比三种方式应接近机器精度 res [norm(A*x1 - b)/norm(b), ... norm(A*x2 - b)/norm(b), ... norm(A*x3 - b)/norm(b)]; disp(res);逻辑说明lu返回的P是行置换矩阵所以回代顺序是P*b不能直接写L\(U\b)否则结果错得离谱。chol要求矩阵对称正定gallery(lehmer, 8)返回的 Lehmer 矩阵满足这个条件但其特征值相差很大、条件数高非常适合用来暴露数值问题。三种方式求出的x1、x2、x3应当高度一致相对残差在1e-14附近属于正常。参数说明A\b不只是语法糖它背后是完整的算法调度。理解这一点很重要因为同样的\遇到非方阵时会自动切换语义超定方程A (mn)得到最小二乘解欠定方程得到基础解系中的一个解。所以看到别人的代码里大量使用\时不要急着改成显式分解先弄清楚矩阵的形状和需求。分解方法适用条件主要用途MATLAB函数LU一般方阵线性方程组、行列式计算luQR长方形矩阵、超定方程最小二乘、特征值迭代qrCholesky对称正定快速求解、协方差分解chol3.2 特征值问题全量 eig 与稀疏 eigs 的适用边界特征值计算是矩阵分析的另一个高频主题PCA、模态分析、图拉普拉斯谱聚类本质都是特征值分解问题。小规模矩阵直接用eig返回全部特征值当矩阵维度上万且稀疏只关心少数几个模最大或接近某个目标值的特征值时要换eigs。d eig(A); % 全部特征值8 阶矩阵直接算 opts.tol 1e-6; % 迭代法收敛容差 [V, D] eigs(A, 6, largestabs, opts); % 求模最大的 6 个特征值 res_eig norm(A*V - V*D, fro) / norm(A, fro);逻辑说明eigs返回的特征向量按列排列在V中D是对角阵A*V - V*D就是残差矩阵。这里用 Frobenius 范数做归一化得到相对残差。这个残差如果超过1e-6通常意味着容差设置过松或矩阵条件数过差需要把opts.tol再压一挡或者检查矩阵是否真的稀疏。参数说明largestabs表示求模最大的特征值对应旧版本里的lm。还可以传smallestabs求模最小、largestreal求实部最大等。对于广义特征值问题Ax λBx直接调用eigs(A, B, k, sigma)即可这在结构动力学和模态分析脚本里很常见。3.3 条件数与解的可信度解矩阵方程、算特征值之后我习惯性检查一个数条件数cond(A)。它衡量方程组解对输入扰动的放大倍数量级接近1e7时解的有效数字大约比机器精度少 7 位双精度下只剩约 9 位有效数字。c cond(A); fprintf(cond(A) %.3e\n, c);逻辑说明cond(A)计算的是最大奇异值与最小奇异值之比是判断矩阵是否病态的最直观指标。配合前面的相对残差能区分两种常见情况残差大但条件数正常说明代码有 bug残差小但条件数极大说明结果本身可信度有限应该考虑预处理、换基或改用精度更高的问题表达方式而不是继续调整求解参数。工程上常见的误用是只看残差不看条件数。残差小只代表A*x接近b并不代表x接近真实解条件数大时一个微小的扰动就能让解发生巨大变化。所以验证一套线性代数求解代码cond和残差要同时看二者缺一不可。4. 数值优化与统计建模从 fmincon 到 ga 与 kmeans4.1 把实际问题规约成标准优化形式源程序进入到examp9_25、examp10_32这类编号话题就转到数值优化。用 MATLAB 优化工具箱做优化第一步不是调函数而是把问题改写成标准形式目标函数f(x)最小化约束统一为线性不等式A*x b、线性等式Aeq*x beq、变量上下界lb/ub、非线性约束c(x)0和ceq(x)0四类。这样做的好处是同一套问题定义可以直接在fmincon、ga、patternsearch之间切换换算法时不需要改约束建模逻辑。以下面的问题为例目标函数是凸二次型约束是圆形区域加一个半平面fun (x) (x(1)-2)^2 (x(2)-1)^2; % 目标函数 x0 [0; 0]; % 初始点 A [-1 -1]; b -1; % 线性不等式x1 x2 1 lb [-10; -10]; ub [10; 10]; % 变量上下界 nonlcon (x) deal(x(1)^2 x(2)^2 - 4, []); % 非线性不等式 0 xopt fmincon(fun, x0, A, b, [], [], lb, ub, nonlcon);逻辑说明fmincon的完整签名是fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon)没有等式约束所以Aeq、beq传空数组。nonlcon返回两个输出第一个是不等式约束的值要求小于等于 0第二个是等式约束要求等于 0用deal可以一条语句完成返回这是源程序里比较常见的写法。参数说明x0是局部优化算法的起点。凸问题从哪出发结果一致非凸问题则几乎由初值决定结局。建议至少取 5 到 10 个随机起点分别求解对比目标函数值而不是只看一次结果。我的做法是用for循环配合fmincon批量跑再取最优值这个方法在源程序的基础上加几行就能实现。4.2 fmincon 的内置算法与关键公差fmincon在近几个主版本里默认算法是 interior-point内点法它对中小规模约束优化问题表现稳定求解器按约束类型自动选择子算法。常用设置是把显示选项和容差显式写出来方便排查收敛过程opts optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... OptimalityTolerance, 1e-8, ... ConstraintTolerance, 1e-8, ... MaxIterations, 1000); xopt fmincon(fun, x0, A, b, [], [], lb, ub, nonlcon, opts);参数说明OptimalityTolerance控制一阶最优性条件的容差太小会让求解器在最优解附近反复迭代ConstraintTolerance控制约束违反程度工程问题1e-6已经足够实验室精度追求才压到1e-10以下。Display设为iter能在命令窗口实时观察迭代步长和目标值变化排查不收敛问题时比只看最终结果有效得多。如果目标函数不可导内点法会报梯度相关错误这时可以改用Algorithm, sqp它对目标函数的光滑性要求更低。但更好的策略是直接换到全局搜索类的无导数算法也就是下一节的内容。4.3 全局搜索兜底ga 遗传算法的参数经验当目标存在多个局部极值fmincon的多起点策略仍然可能漏掉全局最优这时我会上ga。遗传算法不需要梯度对目标函数形态几乎无要求代价是收敛慢且结果有随机性。这个取舍在高等应用数学问题里很常见能接受计算时间换来的是更稳的全局搜索能力。opts optimoptions(ga, ... PopulationSize, 100, ... MaxGenerations, 200, ... Display, iter); [xg, fval] ga(fun, 2, A, b, [], [], lb, ub, nonlcon, opts);逻辑说明ga的第二个参数是变量个数这里fun是二维问题所以传 2。约束建模方式与fmincon完全一致这也是前面强调标准形式的原因。ga不需要初值x0它用随机种群启动所以每次运行结果会有差异。参数说明PopulationSize默认 50非线性约束强时可以升到 100 到 200MaxGenerations默认是变量数的 100 倍这里二维问题就是 200 代。Displayiter可以观察每次迭代的最优值和平均适应度如果连续多代最优值不变说明已经收敛。正式使用前用rng(42)固定随机种子方便复现和对比调参。求解器适用场景是否需梯度随机性关键参数fminunc无约束光滑问题需要无OptimalityTolerancefmincon有约束光滑问题需要无ConstraintTolerancega非光滑/离散/多极值不需要有PopulationSize, MaxGenerationspatternsearch非光滑小规模不需要无MeshSize, PollMethod4.4 回归、概率分布与 kmeans 聚类的配合优化之外源程序中也有大量统计相关内容。线性回归用fitlm一步到位返回系数、p 值和置信区间概率分布拟合用fitdist聚类问题默认kmeans。这三者经常串联使用先做回归理解变量关系再对残差做分布假设检验最后对样本聚类找分组结构。rng(42); data [randn(100, 2) [2 2]; randn(100, 2) - [2 2]]; [idx, C] kmeans(data, 2, Replicates, 5); gscatter(data(:, 1), data(:, 2), idx); pd fitdist(data(:, 1), Normal); % 单变量正态分布拟合 mdl fitlm(data(:, 1), data(:, 2)); % 一元线性回归逻辑说明kmeans是基于距离的划分迭代算法对初始中心敏感。Replicates参数让 MATLAB 从多个随机初始中心出发返回总组内距离最小的解。fitdist返回的概率分布对象可以直接调用pdf、cdf和paramci等方法适合后续做置信区间分析。fitlm的输出里包含系数估计、标准误和 F 统计量是回归分析的统一入口。参数说明如果数据量纲差异大比如一列是上万量级、另一列是 0 到 1 之间务必先标准化否则聚类和回归结果会被量纲大的特征主导。标准化用zscore即可这一步看起来简单实际效果非常关键。数据再复杂一些需要上模型时Deep Learning Toolbox 里的trainNetwork也能接手回归或分类任务但特征量纲处理的原则完全一致这套源程序里体现的建模思路可以直接迁移。5. .asv 文件恢复与运行时间回归源码包调试的实用习惯5.1 .asv 备份文件不是垃圾318 个源程序里混着examp8_20.asv、examp6_20.asv这样的文件这是 MATLAB 编辑器自动保存产生的副本。原.m文件被误删或改坏时把.asv改名为.m就能恢复到最近一次自动保存的状态。编辑器默认每隔几分钟自动保存一份.asv所以备份内容通常只比你最后一次保存晚几分钟。当目录里混了一堆.asv需要统一处理时批量重命名即可。Linux/macOS 终端里执行for f in *.asv; do mv $f ${f%.asv}.m; doneWindows 的 PowerShell 里对应写法是Get-ChildItem *.asv | Rename-Item -NewName { $_.Name -replace \.asv$, .m }。需要注意若同名.m已存在先对比文件修改时间再决定覆盖方向别拿旧备份覆盖新代码。5.2 数值结果的残差验证与 assert 断言我习惯在脚本尾部写断言用机器判断结果是否可信而不是靠肉眼盯数值。对前面的 ODE 例子yRef double(arrayfun((t) ySol(t), tn)); assert(norm(yn - yRef, inf) 1e-4, 数值解与解析解偏差过大);逻辑说明arrayfun把符号解ySol映射到数值节点tn上double完成符号到双精度转换norm(..., inf)取最大偏差。assert第一个条件为假时直接抛异常并输出提示这样批量跑源程序时万一某个参数改动让结果发散能立刻定位到具体脚本而不是一路跑到底再回来肉眼排查。5.3 用 timeit 替代 tic/toc 计时在优化算法或大规模矩阵计算里算法改进效果不能靠一次tic/toc判断因为它受系统负载影响大。timeit会多次运行函数句柄并取稳健统计值更适合做性能回归测试f1 () ode45((t, y) -2*y t, [0 5], 1); t1 timeit(f1);参数说明timeit的输入必须是无参数或参数已捕获的函数句柄运行次数由它内部决定一般会在几百到几千次之间权衡精度与耗时。对比两个算法时分别timeit并重复三次以上再取均值结论比单次tic/toc可靠得多。如果需要更大的性能收益把热点代码写成 C 再编译成 mex 文件是 MATLAB 里运行外部语言的标准做法能把循环密集的数学内核提速数倍。另外注意一个数值显示陷阱1e100在命令行直接显示为1.0000e100继续参与某些计算时容易溢出成Inf。用sym(1e100)包一层即可看到完整精确值避免在调试大数问题时被浮点显示误导。本文还有配套的精品资源点击获取