尧图网站建设 
1. 项目缘起为什么需要一份自己的代码记录又到了数学建模竞赛的备赛季或者说是期末复习周。看着电脑里散落各处的.m文件是不是感觉头都大了这个文件夹里是去年美赛的线性规划求解那个压缩包里是上次课程大作业的非线性拟合还有一个叫test_final.m的文件打开一看里面空空如也只有一句disp(hello world)。我相信这是很多同学包括曾经的我都经历过的窘境。数学建模的核心一半在模型另一半就在实现。模型想得再天花乱坠代码跑不通一切都是零。而代码跑通之后如果没有系统性地记录下来下次遇到类似问题你大概率还是会从头开始搜索、调试、踩坑。这就是我整理这份“数学建模复习——代码记录”的初衷。它不是一个面面俱到的 MATLAB 百科全书也不是一本高深的算法理论书。它更像是我个人在多次竞赛和项目实战后沉淀下来的一个“工具箱”和“错题本”。里面记录的都是那些最常用、最容易出错、但又最核心的代码片段和思路。当你在深夜被一个规划问题卡住或者不确定该用哪种统计检验方法时希望这份记录能像一位随时待命的队友给你最直接、最可靠的参考。2. 规划问题从线性到0-1的实战代码框架规划问题是数学建模的基石无论是优化资源分配还是寻找最佳路径都离不开它。很多人觉得调用linprog或fmincon函数很简单但真正导致结果错误或求解失败的往往是模型转化和选项设置的细节。2.1 线性规划标准形式与linprog的陷阱线性规划LP的标准形式是求最小值min f*x约束条件为A*x b,Aeq*x beq,lb x ub。第一个坑就在这里你的模型是否转化成了标准形式比如如果你的目标是最大化c*x那么目标函数系数f应该是-c。% 示例生产计划问题 % 目标最大化利润 z 3*x1 5*x2 % 约束x1 4, 2*x2 12, 3*x1 2*x2 18, x1, x2 0 f [-3; -5]; % 最大化转为最小化 A [1, 0; 0, 2; 3, 2]; b [4; 12; 18]; Aeq []; beq []; lb [0; 0]; ub []; options optimoptions(linprog, Display, iter); % 显示迭代过程调试时非常有用 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub, options); if exitflag 0 fprintf(最优解: x1 %.2f, x2 %.2f\n, x(1), x(2)); fprintf(最大利润: %.2f\n, -fval); % 注意fval是转换后目标函数的最小值 else fprintf(求解失败。退出标志: %d\n, exitflag); fprintf(输出信息: %s\n, output.message); end注意exitflag的值至关重要。1表示收敛到解0表示达到最大迭代次数-2表示无可行解-3表示问题无界。务必在代码中判断exitflag而不是直接相信x的结果。我曾因为忽略了这个标志把一个无解的问题当成了有解导致后续分析全部错误。2.2 非线性规划初值、梯度与算法选择非线性规划NLP的复杂性呈指数级上升。fmincon是主力函数但它的表现极度依赖于你给的“起点”——初始值x0。% 示例求解 Rosenbrock 函数最小值一个经典测试函数 % f(x) 100*(x2 - x1^2)^2 (1 - x1)^2 % 初始点设为 [-1, 2] fun (x) 100*(x(2) - x(1)^2)^2 (1 - x(1))^2; x0 [-1, 2]; A []; b []; Aeq []; beq []; lb []; ub []; % 尝试不同的算法 options_fmincon optimoptions(fmincon, Display, final, Algorithm, interior-point); [x1, fval1] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, [], options_fmincon); options_sqp optimoptions(fmincon, Display, final, Algorithm, sqp); [x2, fval2] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, [], options_sqp); fprintf(内点法结果: x [%.6f, %.6f], fval %.6e\n, x1(1), x1(2), fval1); fprintf(SQP法结果: x [%.6f, %.6f], fval %.6e\n, x2(1), x2(2), fval2);这里的关键心得是永远不要只从一个初始点开始求解。对于复杂问题最好能从一个粗略的网格中选取多个初始点分别求解然后比较结果。这能有效避免陷入局部最优解。另外如果目标函数或约束条件是可导的提供解析梯度或雅可比矩阵给fmincon通过‘SpecifyObjectiveGradient’和‘SpecifyConstraintGradient’选项能大幅提高求解速度和稳定性。自己手推梯度很麻烦但用符号计算工具箱 (syms) 可以帮你自动求导并生成函数句柄。2.3 多目标规划帕累托前沿的生成与理解多目标规划没有唯一的“最优解”只有一系列“非劣解”帕累托最优解。MATLAB 中常用gamultiobj基于遗传算法来求解。% 示例经典ZDT1测试问题 % 目标最小化 f1(x) x1, f2(x) g(x) * (1 - sqrt(x1/g(x))) % 其中 g(x) 1 9/(n-1) * sum_{i2}^{n} x_i, n30, x_i in [0,1] nvars 30; % 变量个数 A []; b []; Aeq []; beq []; lb zeros(1, nvars); ub ones(1, nvars); % 定义多目标函数 multiFun (x) [x(1), ... (1 9*sum(x(2:end))/(nvars-1)) * (1 - sqrt(x(1)/(1 9*sum(x(2:end))/(nvars-1))))]; options optimoptions(gamultiobj, PopulationSize, 100, ParetoFraction, 0.35, ... PlotFcn, gaplotpareto); % 绘制帕累托前沿 [x, fval, exitflag, output] gamultiobj(multiFun, nvars, A, b, Aeq, beq, lb, ub, options);运行后fval是一个矩阵每一行是一个解对应的两个目标函数值。可视化这些点就能得到帕累托前沿——一条或一片权衡曲线。最重要的理解是在这条曲线上任何一个目标的改进必然以另一个目标的恶化为代价。你的任务不是找到“最好”的点而是根据实际问题的偏好例如更看重成本还是效率从曲线上选择一个合适的折中点。遗传算法是启发式搜索每次运行结果可能略有不同可以适当增大‘PopulationSize’和‘Generations’来获得更稳定的前沿。2.4 0-1规划整数约束的转化与求解0-1规划是整数规划的特例变量只能取0或1。MATLAB 使用intlinprog求解混合整数线性规划MILP。关键是指定哪些变量是整数变量。% 示例背包问题简化版 % 有5件物品重量w[2,3,4,5,9]价值v[3,4,5,8,10]背包容量C15。 % 决策变量 xi 0 或 1表示是否选择第i件物品。 % 目标最大化总价值。 f -[3; 4; 5; 8; 10]; % 最大化转为最小化 A [2, 3, 4, 5, 9]; % 重量约束 b 15; lb zeros(5,1); ub ones(5,1); intcon 1:5; % 关键指定所有5个变量都是整数变量此处为0-1 [x, fval, exitflag] intlinprog(f, intcon, A, b, [], [], lb, ub); if exitflag 0 fprintf(选择的物品编号: %s\n, mat2str(find(x 0.5))); fprintf(最大总价值: %d\n, -fval); end对于非线性0-1规划fmincon不再适用。一种常见做法是使用遗传算法 (ga)并将变量的上下界设为 [0,1]然后通过自定义的整数约束处理函数或者在适应度函数中惩罚非整数解来逼近。但更严谨的做法是使用专门的混合整数非线性规划MINLP求解器如 OPTI Toolbox 搭配的 SCIP 或 BONMIN但这需要额外安装配置。在时间紧张的竞赛中如果问题规模不大用ga并设置‘PopulationType’为‘bitstring’或‘custom’来直接操作二进制变量是一个快速的备选方案。3. 数据处理与统计分析MATLAB 作为计算引擎建模离不开数据。清洗数据、分析规律、验证假设这一套流程在 MATLAB 里可以非常流畅地完成。3.1 假设检验t-test 的正确打开方式搜索热词里提到了ttest和ttest2的区别这确实是个高频困惑点。ttest(单样本或配对t检验)检验一组数据的均值是否等于某个给定值或者检验两组配对数据的差值均值是否为零。% 单样本t检验检验一组学生的成绩均值是否为75分 scores [78, 82, 65, 90, 71, 88, 76, 85]; [h, p, ci, stats] ttest(scores, 75); % h1 拒绝原假设均值不等于75p值很小ci是均值的置信区间 % 配对样本t检验检验某种训练方法前后学生的成绩是否有差异 before [70, 65, 80, 75, 60]; after [75, 68, 85, 80, 70]; [h, p] ttest(before, after); % 默认是双侧检验检验均值差是否为0ttest2(双样本t检验)检验两个独立样本组的均值是否相等。它假设两组数据来自方差相等的正态总体同方差如果方差不齐需要使用‘Vartype’参数。% 双样本t检验比较两种教学方法下学生的成绩 groupA [78, 82, 65, 90, 71]; groupB [88, 76, 85, 92, 80, 83]; [h, p, ci, stats] ttest2(groupA, groupB); % 如果怀疑方差不齐使用 [h2, p2] ttest2(groupA, groupB, Vartype, unequal);核心经验在使用前一定要用vartest2或levene检验一下两组数据的方差齐性。直接使用默认的同方差ttest2如果数据方差异质性很大可能导致错误的结论。另外p值小于显著性水平如0.05只能说明“有统计学差异”不能说明“差异很大”。效应量如 Cohen‘s d才是衡量差异大小的指标MATLAB 需要手动计算。3.2 图像处理从基础操作到特征提取数学建模中图像常作为数据源。除了基本的读 (imread)、写 (imwrite)、显示 (imshow)更重要的是预处理和特征提取。% 示例读取图像转为灰度图进行边缘检测和阈值分割 img imread(sample.jpg); if size(img, 3) 3 grayImg rgb2gray(img); else grayImg img; end % 边缘检测Canny算子效果较好 edges edge(grayImg, canny); figure; subplot(2,2,1); imshow(img); title(原图); subplot(2,2,2); imshow(grayImg); title(灰度图); subplot(2,2,3); imshow(edges); title(边缘检测); % 阈值分割Otsu方法自动确定阈值 level graythresh(grayImg); bwImg imbinarize(grayImg, level); subplot(2,2,4); imshow(bwImg); title(二值化分割); % 计算一些简单的区域属性用于后续建模 stats regionprops(bwImg, Area, Centroid, BoundingBox); disp([检测到 , num2str(numel(stats)), 个连通区域]);踩坑记录imread读入的图像其数据类型可能是uint8(0-255)。在进行数学运算如减法、滤波前最好用im2double或double()将其转换为双精度浮点数否则可能导致溢出或精度问题。另外edge函数非常敏感其中的阈值参数需要根据图像对比度仔细调整默认参数往往不是最优的。3.3 信息熵计算量化数据的不确定性信息熵是衡量系统混乱度或信息量的重要指标在评估数据分布、特征选择中很有用。% 计算一维离散数据的信息熵 % 假设 data 是一个向量包含离散的类别或符号 data [1, 2, 1, 3, 2, 1, 1, 2, 3, 3, 3, 1]; % 方法计算每个唯一值出现的概率然后套用熵公式 H -sum(p.*log2(p)) uniqueVals unique(data); counts histcounts(data, [uniqueVals, max(uniqueVals)1]); % 统计频数 probs counts / sum(counts); probs probs(probs 0); % 去掉概率为0的项因为0*log2(0)无定义 entropy_val -sum(probs .* log2(probs)); fprintf(数据的信息熵为: %.4f bits\n, entropy_val); % 封装成函数 function H calc_entropy(data) [counts, ~] histcounts(data, BinMethod, integers); probs counts / sum(counts); probs probs(probs 0); H -sum(probs .* log2(probs)); end这个计算的核心在于概率估计。对于连续数据需要先进行离散化分箱。分箱的数目‘BinMethod’会显著影响熵值。 Sturges‘ 准则或 Freedman-Diaconis 准则是常用的自动分箱方法。在特征选择中我们常计算每个特征与目标变量的互信息这需要计算联合概率分布思路类似但维度更高。4. 高级应用与仿真Simulink 与算法实现数学建模不止于静态模型动态系统仿真和复杂算法实现是更高阶的要求。4.1 Simulink 基础与联合仿真Simulink 是进行系统动力学、控制逻辑仿真的利器。对于数学建模即使不涉及复杂的控制系统用它来搭建一个微分方程模型也直观得多。搭建一个简单的微分方程模型比如人口增长的 Logistic 模型dx/dt r*x*(1 - x/K)。在 Simulink 库中拖入一个Integrator模块其输入是导数dx/dt输出是x。用Gain(增益)、Sum(求和)、Product(乘积) 模块和Constant(常数) 模块搭建右侧的表达式r*x*(1 - x/K)。将表达式的结果连接到Integrator的输入口。设置仿真时间、初始条件在Integrator模块中设置然后运行。与 MATLAB 工作空间交互这是最实用的技巧。在 Simulink 模型中使用From Workspace模块读入 MATLAB 变量作为输入使用To Workspace模块将仿真结果输出到 MATLAB 工作空间。这样你就可以用 MATLAB 脚本预处理数据、批量运行仿真、并后处理结果。联合仿真搜索热词中提到了 Adams机械动力学软件与 MATLAB 的联合仿真。这通常通过S-Function(系统函数) 实现。Adams 可以生成一个描述机械系统的M文件 S-Function将其导入 Simulink 作为一个模块。Simulink 负责控制算法部分两者通过定义好的输入输出接口进行数据交换实现机电一体化系统的协同仿真。这个过程配置较为复杂需要确保两个软件的版本兼容并正确设置求解器参数。4.2 随机过程模拟醉汉游走模型随机游走是模拟许多自然和社会现象的经典模型如股票价格、粒子扩散等。% 一维简单随机游走醉汉游走 num_steps 1000; % 总步数 num_walkers 5; % 模拟多个醉汉 paths zeros(num_walkers, num_steps1); % 记录路径 for w 1:num_walkers position 0; paths(w, 1) position; for step 1:num_steps % 每一步以相等概率向左或向右移动一个单位 move randi([0, 1]) * 2 - 1; % 生成 -1 或 1 position position move; paths(w, step1) position; end end % 可视化 figure; hold on; for w 1:num_walkers plot(0:num_steps, paths(w, :)); end xlabel(步数); ylabel(位置); title(一维随机游走模拟); grid on; hold off; % 分析计算最终位置的均值和方差 final_positions paths(:, end); fprintf(最终位置均值: %.2f\n, mean(final_positions)); fprintf(最终位置方差: %.2f\n, var(final_positions)); % 理论上均值应为0方差约等于步数。这个模型可以轻松扩展变成二维游走每次随机选择东、南、西、北加入偏向性例如向右走的概率更大或者让步长服从某种分布如正态分布。在建模中它常用于作为更复杂模型的基准或组成部分例如在模拟搜索行为、市场波动时。4.3 自定义算法实现隐式 QR 方法求特征值当需要深入理解算法原理或者内置函数eig无法满足特殊需求时自己实现算法就很有必要。隐式 QR 方法是计算中小型矩阵全部特征值的稳定方法。function eigenvalues my_implicit_qr(A, max_iter, tol) % 简化的隐式QR算法实现带实用主义简化非工业级 % 输入方阵A最大迭代次数容忍度 % 输出特征值向量 n size(A, 1); H hess(A); % 先上Hessenberg化大幅减少计算量 eigenvalues zeros(n, 1); iter 0; while n 1 iter max_iter % 实用技巧检查次对角线元素是否可忽略 for i n:-1:2 if abs(H(i, i-1)) tol * (abs(H(i-1, i-1)) abs(H(i, i))) H(i, i-1) 0; % 发现可分离的子问题 if i n % 右下角1x1块得到一个特征值 eigenvalues(n) H(n, n); n n - 1; H H(1:n, 1:n); else % 分离为两个更小的子问题这里简化处理递归调用实际应迭代处理 % 此处为演示我们跳出循环用内置eig处理剩余部分 remaining_vals eig(H(1:n, 1:n)); eigenvalues(1:length(remaining_vals)) remaining_vals; return; end break; end end % 单步隐式QR位移使用Wilkinson位移以加速收敛 mu wilkinson_shift(H(n-1:n, n-1:n)); [Q, R] qr(H(1:n, 1:n) - mu * eye(n)); H(1:n, 1:n) R * Q mu * eye(n); iter iter 1; end if n 1 eigenvalues(1) H(1,1); end % 清理接近零的虚部 eigenvalues real(eigenvalues) 1i * (abs(imag(eigenvalues)) tol) .* imag(eigenvalues); end function mu wilkinson_shift(B) % 计算2x2矩阵的Wilkinson位移 a B(1,1); b B(1,2); d B(2,2); delta (a - d) / 2; mu d - sign(delta) * b^2 / (abs(delta) sqrt(delta^2 b^2)); end重要提醒自己实现数值线性代数算法主要是为了教学和理解。在真正的建模竞赛或工程中永远优先使用 MATLAB 内置的、经过高度优化的函数如eig,svd,chol等。它们更快、更稳定、更准确。自己写的代码可能在某些边界条件下失效比如几乎奇异矩阵。这里的实现是一个高度简化的教学版本重点展示了算法框架、Hessenberg 化降维、Wilkinson 位移加速和 deflation收缩技巧这些核心思想。5. 效率优化与调试让代码跑得更快更稳在建模中时间就是生命。一个需要跑一小时的脚本和另一个功能相同但只需一分钟的脚本有本质区别。5.1 向量化操作告别缓慢的 for 循环这是 MATLAB 性能提升的第一法则。尽可能使用矩阵和向量运算代替循环。% 低效做法计算两个向量对应元素的某种函数值例如距离计算 n 1e6; x randn(n, 1); y randn(n, 1); result_slow zeros(n, 1); tic; for i 1:n result_slow(i) sqrt(x(i)^2 y(i)^2); end time_slow toc; % 高效做法向量化 tic; result_fast sqrt(x.^2 y.^2); % 点乘 .^ 和 .* 是关键 time_fast toc; fprintf(循环耗时: %.4f 秒\n, time_slow); fprintf(向量化耗时: %.4f 秒\n, time_fast); fprintf(加速比: %.2f 倍\n, time_slow / time_fast);对于更复杂的多重循环嵌套可以考虑使用meshgrid或ndgrid生成网格然后一次性进行向量运算。但要注意这会瞬间创建巨大的矩阵可能耗尽内存。需要在向量化和内存消耗之间权衡。5.2 预分配数组避免动态增长在循环中不断向一个数组追加元素例如result [result; new_value]是性能杀手因为 MATLAB 每次都需要寻找新的连续内存并复制整个数组。% 糟糕的做法 result []; for k 1:10000 result [result; some_calculation(k)]; % 每次循环都重新分配内存 end % 正确的做法 n 10000; result zeros(n, 1); % 预先分配好内存 for k 1:n result(k) some_calculation(k); % 直接赋值 end使用zeros,ones,NaN,inf等函数进行预分配。如果不知道最终大小可以高估一个上限计算完后再截断。5.3 函数句柄与匿名函数提升代码灵活性函数句柄允许你将函数作为参数传递这在优化、积分和绘图时非常有用。% 定义目标函数用于 fminsearch rosenbrock (x) 100*(x(2)-x(1)^2)^2 (1-x(1))^2; x0 [-1.2, 1]; [x_opt, fval] fminsearch(rosenbrock, x0); % 在数组函数中使用如 arrayfun, cellfun但需谨慎有时不如直接向量化快 x_cell num2cell(1:10); squares cellfun((a) a^2, x_cell); % 对元胞数组每个元素平方 % 创建带参数的函数句柄 a 2; b 3; linear_func (x) a*x b; % 捕获了当前工作空间的 a 和 b y linear_func(5); % y 2*5313匿名函数非常简洁但过于复杂的逻辑不适合写在里面会影响可读性。对于复杂的、需要复用的功能还是应该写在独立的.m函数文件中。5.4 调试与错误排查从报错信息中快速定位MATLAB 的错误信息尤其是新版已经相当友好了。关键是要学会看。“函数或变量无法识别”比如热词中的‘deltalin’。首先检查拼写。然后检查该函数所在的文件夹是否在 MATLAB 的搜索路径中。使用which deltalin命令查看 MATLAB 是否能找到它。如果是一个自定义函数确保它的文件名和函数名一致并且当前目录或路径正确。矩阵维度不匹配这是最常见的错误之一。出错时MATLAB 会告诉你操作涉及矩阵的具体维度。立刻使用whos或在工作区查看相关变量的size。确保乘法 (*和.*)、加法、连接等操作的维度是兼容的。使用断点和调试模式在代码行号左侧点击设置断点红色圆点。运行代码时执行到断点处会暂停。此时你可以将鼠标悬停在变量上查看其当前值也可以在命令窗口直接输入变量名查看。使用F10单步执行F11步入函数F5继续运行。这是定位逻辑错误最强大的工具。try-catch块在可能出错的代码段外围加上try-catch可以捕获错误并执行备用方案避免程序完全崩溃。try data load(unreliable_data_file.mat); result complex_analysis(data); catch ME % ME 是一个包含错误信息的结构体 warning(加载或分析数据失败: %s, ME.message); result []; % 返回空值或默认值 % 可以记录错误日志fprintf(log_file, %s: %s\n, datetime, ME.message); endtic/toc和性能分析器用tic和toc包围你觉得慢的代码段来计时。对于更全面的性能分析使用profile工具。在命令窗口输入profile on运行你的代码然后输入profile viewer。它会生成一个详细的报告告诉你每行代码的执行时间和调用次数精准定位性能瓶颈。6. 可视化与成果输出让图表和论文说话好的可视化能让你和评委快速理解模型结果。而规范的输出则是最终成果的保障。6.1 高级绘图技巧超越 plot双Y轴图使用yyaxis。x linspace(0, 10); y1 sin(x); y2 exp(x/5); figure; yyaxis left; plot(x, y1, b-); ylabel(sin(x)); yyaxis right; plot(x, y2, r--); ylabel(exp(x/5)); xlabel(x); title(双Y轴示例);子图灵活排列使用tiledlayout替代旧的subplot它能更好地控制间距和标题。figure; t tiledlayout(2, 2); % 2行2列 nexttile; plot(rand(10,1)); title(图1); nexttile; scatter(rand(10,1), rand(10,1)); title(图2); nexttile(3, [1, 2]); % 跨一列占两列宽度 surf(peaks); title(跨列子图); xlabel(t, 共同的X轴标签); ylabel(t, 共同的Y轴标签); title(t, 整体标题);颜色映射与数据映射使用colormap来根据数据值着色这在显示矩阵、三维曲面时非常直观。data randn(20, 20); imagesc(data); colorbar; % 显示颜色条 colormap(jet); % 或 hot, cool, parula (默认) axis image; % 使坐标轴比例相等6.2 导出高质量图片用于论文的 EPS/PDF搜索热词里提到了导出 EPS。对于学术出版矢量格式EPS, PDF是首选因为它们无限缩放不失真。figure; plot(1:10, rand(1,10), o-, LineWidth, 2, MarkerSize, 8); xlabel(Index); ylabel(Value); title(示例图); grid on; % 方法1直接使用 exportgraphics (R2020a及以上推荐) exportgraphics(gcf, my_plot.pdf, ContentType, vector); % 导出为矢量PDF % exportgraphics(gcf, my_plot.eps, ContentType, vector); % 导出为EPS % 方法2使用 print 函数旧版本通用 set(gcf, Renderer, painters); % 确保使用矢量渲染器 print(-depsc, -tiff, -r300, my_plot.eps); % 导出EPS嵌入TIFF预览分辨率300dpi % print(-dpdf, -bestfit, my_plot.pdf); % 导出PDF自动适应页面 % 方法3调整画布大小再导出 fig gcf; fig.PaperPositionMode auto; % 使用屏幕尺寸 fig_pos fig.PaperPosition; fig.PaperSize [fig_pos(3) fig_pos(4)]; % 设置纸张大小与图形一致 print(fig, my_plot.pdf, -dpdf, -bestfit);注意如果图形中包含大量数据点如十万个散点导出为矢量文件可能会非常大。此时可以导出为高分辨率位图如PNG,‘-r600’或者使用exportgraphics并指定‘Resolution’参数。在导出前务必检查图形中所有文字标签、图例、标题的字体和大小是否合适因为一旦导出为图片就很难再修改。6.3 生成报告与文档发布你的代码MATLAB 的“发布”Publish功能可以将脚本.m文件转换为格式良好的报告HTML, Word, PDF 等包含代码、运行结果和图表。在编辑器里写好你的脚本用%%创建分节单元格。在每个单元格中先写注释Markdown格式再写代码。点击“发布”选项卡上的按钮选择输出格式。MATLAB 会按顺序执行每个单元格并将输出包括文本、变量值、图形嵌入到生成的文档中。这对于撰写建模论文的附录、制作实验报告或分享可复现的研究结果极其有用。它能确保文档中的结果与代码完全对应。 数学建模实战：MATLAB代码工具箱与调试优化指南
1. 项目缘起为什么需要一份自己的代码记录又到了数学建模竞赛的备赛季或者说是期末复习周。看着电脑里散落各处的.m文件是不是感觉头都大了这个文件夹里是去年美赛的线性规划求解那个压缩包里是上次课程大作业的非线性拟合还有一个叫test_final.m的文件打开一看里面空空如也只有一句disp(hello world)。我相信这是很多同学包括曾经的我都经历过的窘境。数学建模的核心一半在模型另一半就在实现。模型想得再天花乱坠代码跑不通一切都是零。而代码跑通之后如果没有系统性地记录下来下次遇到类似问题你大概率还是会从头开始搜索、调试、踩坑。这就是我整理这份“数学建模复习——代码记录”的初衷。它不是一个面面俱到的 MATLAB 百科全书也不是一本高深的算法理论书。它更像是我个人在多次竞赛和项目实战后沉淀下来的一个“工具箱”和“错题本”。里面记录的都是那些最常用、最容易出错、但又最核心的代码片段和思路。当你在深夜被一个规划问题卡住或者不确定该用哪种统计检验方法时希望这份记录能像一位随时待命的队友给你最直接、最可靠的参考。2. 规划问题从线性到0-1的实战代码框架规划问题是数学建模的基石无论是优化资源分配还是寻找最佳路径都离不开它。很多人觉得调用linprog或fmincon函数很简单但真正导致结果错误或求解失败的往往是模型转化和选项设置的细节。2.1 线性规划标准形式与linprog的陷阱线性规划LP的标准形式是求最小值min f*x约束条件为A*x b,Aeq*x beq,lb x ub。第一个坑就在这里你的模型是否转化成了标准形式比如如果你的目标是最大化c*x那么目标函数系数f应该是-c。% 示例生产计划问题 % 目标最大化利润 z 3*x1 5*x2 % 约束x1 4, 2*x2 12, 3*x1 2*x2 18, x1, x2 0 f [-3; -5]; % 最大化转为最小化 A [1, 0; 0, 2; 3, 2]; b [4; 12; 18]; Aeq []; beq []; lb [0; 0]; ub []; options optimoptions(linprog, Display, iter); % 显示迭代过程调试时非常有用 [x, fval, exitflag, output] linprog(f, A, b, Aeq, beq, lb, ub, options); if exitflag 0 fprintf(最优解: x1 %.2f, x2 %.2f\n, x(1), x(2)); fprintf(最大利润: %.2f\n, -fval); % 注意fval是转换后目标函数的最小值 else fprintf(求解失败。退出标志: %d\n, exitflag); fprintf(输出信息: %s\n, output.message); end注意exitflag的值至关重要。1表示收敛到解0表示达到最大迭代次数-2表示无可行解-3表示问题无界。务必在代码中判断exitflag而不是直接相信x的结果。我曾因为忽略了这个标志把一个无解的问题当成了有解导致后续分析全部错误。2.2 非线性规划初值、梯度与算法选择非线性规划NLP的复杂性呈指数级上升。fmincon是主力函数但它的表现极度依赖于你给的“起点”——初始值x0。% 示例求解 Rosenbrock 函数最小值一个经典测试函数 % f(x) 100*(x2 - x1^2)^2 (1 - x1)^2 % 初始点设为 [-1, 2] fun (x) 100*(x(2) - x(1)^2)^2 (1 - x(1))^2; x0 [-1, 2]; A []; b []; Aeq []; beq []; lb []; ub []; % 尝试不同的算法 options_fmincon optimoptions(fmincon, Display, final, Algorithm, interior-point); [x1, fval1] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, [], options_fmincon); options_sqp optimoptions(fmincon, Display, final, Algorithm, sqp); [x2, fval2] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, [], options_sqp); fprintf(内点法结果: x [%.6f, %.6f], fval %.6e\n, x1(1), x1(2), fval1); fprintf(SQP法结果: x [%.6f, %.6f], fval %.6e\n, x2(1), x2(2), fval2);这里的关键心得是永远不要只从一个初始点开始求解。对于复杂问题最好能从一个粗略的网格中选取多个初始点分别求解然后比较结果。这能有效避免陷入局部最优解。另外如果目标函数或约束条件是可导的提供解析梯度或雅可比矩阵给fmincon通过‘SpecifyObjectiveGradient’和‘SpecifyConstraintGradient’选项能大幅提高求解速度和稳定性。自己手推梯度很麻烦但用符号计算工具箱 (syms) 可以帮你自动求导并生成函数句柄。2.3 多目标规划帕累托前沿的生成与理解多目标规划没有唯一的“最优解”只有一系列“非劣解”帕累托最优解。MATLAB 中常用gamultiobj基于遗传算法来求解。% 示例经典ZDT1测试问题 % 目标最小化 f1(x) x1, f2(x) g(x) * (1 - sqrt(x1/g(x))) % 其中 g(x) 1 9/(n-1) * sum_{i2}^{n} x_i, n30, x_i in [0,1] nvars 30; % 变量个数 A []; b []; Aeq []; beq []; lb zeros(1, nvars); ub ones(1, nvars); % 定义多目标函数 multiFun (x) [x(1), ... (1 9*sum(x(2:end))/(nvars-1)) * (1 - sqrt(x(1)/(1 9*sum(x(2:end))/(nvars-1))))]; options optimoptions(gamultiobj, PopulationSize, 100, ParetoFraction, 0.35, ... PlotFcn, gaplotpareto); % 绘制帕累托前沿 [x, fval, exitflag, output] gamultiobj(multiFun, nvars, A, b, Aeq, beq, lb, ub, options);运行后fval是一个矩阵每一行是一个解对应的两个目标函数值。可视化这些点就能得到帕累托前沿——一条或一片权衡曲线。最重要的理解是在这条曲线上任何一个目标的改进必然以另一个目标的恶化为代价。你的任务不是找到“最好”的点而是根据实际问题的偏好例如更看重成本还是效率从曲线上选择一个合适的折中点。遗传算法是启发式搜索每次运行结果可能略有不同可以适当增大‘PopulationSize’和‘Generations’来获得更稳定的前沿。2.4 0-1规划整数约束的转化与求解0-1规划是整数规划的特例变量只能取0或1。MATLAB 使用intlinprog求解混合整数线性规划MILP。关键是指定哪些变量是整数变量。% 示例背包问题简化版 % 有5件物品重量w[2,3,4,5,9]价值v[3,4,5,8,10]背包容量C15。 % 决策变量 xi 0 或 1表示是否选择第i件物品。 % 目标最大化总价值。 f -[3; 4; 5; 8; 10]; % 最大化转为最小化 A [2, 3, 4, 5, 9]; % 重量约束 b 15; lb zeros(5,1); ub ones(5,1); intcon 1:5; % 关键指定所有5个变量都是整数变量此处为0-1 [x, fval, exitflag] intlinprog(f, intcon, A, b, [], [], lb, ub); if exitflag 0 fprintf(选择的物品编号: %s\n, mat2str(find(x 0.5))); fprintf(最大总价值: %d\n, -fval); end对于非线性0-1规划fmincon不再适用。一种常见做法是使用遗传算法 (ga)并将变量的上下界设为 [0,1]然后通过自定义的整数约束处理函数或者在适应度函数中惩罚非整数解来逼近。但更严谨的做法是使用专门的混合整数非线性规划MINLP求解器如 OPTI Toolbox 搭配的 SCIP 或 BONMIN但这需要额外安装配置。在时间紧张的竞赛中如果问题规模不大用ga并设置‘PopulationType’为‘bitstring’或‘custom’来直接操作二进制变量是一个快速的备选方案。3. 数据处理与统计分析MATLAB 作为计算引擎建模离不开数据。清洗数据、分析规律、验证假设这一套流程在 MATLAB 里可以非常流畅地完成。3.1 假设检验t-test 的正确打开方式搜索热词里提到了ttest和ttest2的区别这确实是个高频困惑点。ttest(单样本或配对t检验)检验一组数据的均值是否等于某个给定值或者检验两组配对数据的差值均值是否为零。% 单样本t检验检验一组学生的成绩均值是否为75分 scores [78, 82, 65, 90, 71, 88, 76, 85]; [h, p, ci, stats] ttest(scores, 75); % h1 拒绝原假设均值不等于75p值很小ci是均值的置信区间 % 配对样本t检验检验某种训练方法前后学生的成绩是否有差异 before [70, 65, 80, 75, 60]; after [75, 68, 85, 80, 70]; [h, p] ttest(before, after); % 默认是双侧检验检验均值差是否为0ttest2(双样本t检验)检验两个独立样本组的均值是否相等。它假设两组数据来自方差相等的正态总体同方差如果方差不齐需要使用‘Vartype’参数。% 双样本t检验比较两种教学方法下学生的成绩 groupA [78, 82, 65, 90, 71]; groupB [88, 76, 85, 92, 80, 83]; [h, p, ci, stats] ttest2(groupA, groupB); % 如果怀疑方差不齐使用 [h2, p2] ttest2(groupA, groupB, Vartype, unequal);核心经验在使用前一定要用vartest2或levene检验一下两组数据的方差齐性。直接使用默认的同方差ttest2如果数据方差异质性很大可能导致错误的结论。另外p值小于显著性水平如0.05只能说明“有统计学差异”不能说明“差异很大”。效应量如 Cohen‘s d才是衡量差异大小的指标MATLAB 需要手动计算。3.2 图像处理从基础操作到特征提取数学建模中图像常作为数据源。除了基本的读 (imread)、写 (imwrite)、显示 (imshow)更重要的是预处理和特征提取。% 示例读取图像转为灰度图进行边缘检测和阈值分割 img imread(sample.jpg); if size(img, 3) 3 grayImg rgb2gray(img); else grayImg img; end % 边缘检测Canny算子效果较好 edges edge(grayImg, canny); figure; subplot(2,2,1); imshow(img); title(原图); subplot(2,2,2); imshow(grayImg); title(灰度图); subplot(2,2,3); imshow(edges); title(边缘检测); % 阈值分割Otsu方法自动确定阈值 level graythresh(grayImg); bwImg imbinarize(grayImg, level); subplot(2,2,4); imshow(bwImg); title(二值化分割); % 计算一些简单的区域属性用于后续建模 stats regionprops(bwImg, Area, Centroid, BoundingBox); disp([检测到 , num2str(numel(stats)), 个连通区域]);踩坑记录imread读入的图像其数据类型可能是uint8(0-255)。在进行数学运算如减法、滤波前最好用im2double或double()将其转换为双精度浮点数否则可能导致溢出或精度问题。另外edge函数非常敏感其中的阈值参数需要根据图像对比度仔细调整默认参数往往不是最优的。3.3 信息熵计算量化数据的不确定性信息熵是衡量系统混乱度或信息量的重要指标在评估数据分布、特征选择中很有用。% 计算一维离散数据的信息熵 % 假设 data 是一个向量包含离散的类别或符号 data [1, 2, 1, 3, 2, 1, 1, 2, 3, 3, 3, 1]; % 方法计算每个唯一值出现的概率然后套用熵公式 H -sum(p.*log2(p)) uniqueVals unique(data); counts histcounts(data, [uniqueVals, max(uniqueVals)1]); % 统计频数 probs counts / sum(counts); probs probs(probs 0); % 去掉概率为0的项因为0*log2(0)无定义 entropy_val -sum(probs .* log2(probs)); fprintf(数据的信息熵为: %.4f bits\n, entropy_val); % 封装成函数 function H calc_entropy(data) [counts, ~] histcounts(data, BinMethod, integers); probs counts / sum(counts); probs probs(probs 0); H -sum(probs .* log2(probs)); end这个计算的核心在于概率估计。对于连续数据需要先进行离散化分箱。分箱的数目‘BinMethod’会显著影响熵值。 Sturges‘ 准则或 Freedman-Diaconis 准则是常用的自动分箱方法。在特征选择中我们常计算每个特征与目标变量的互信息这需要计算联合概率分布思路类似但维度更高。4. 高级应用与仿真Simulink 与算法实现数学建模不止于静态模型动态系统仿真和复杂算法实现是更高阶的要求。4.1 Simulink 基础与联合仿真Simulink 是进行系统动力学、控制逻辑仿真的利器。对于数学建模即使不涉及复杂的控制系统用它来搭建一个微分方程模型也直观得多。搭建一个简单的微分方程模型比如人口增长的 Logistic 模型dx/dt r*x*(1 - x/K)。在 Simulink 库中拖入一个Integrator模块其输入是导数dx/dt输出是x。用Gain(增益)、Sum(求和)、Product(乘积) 模块和Constant(常数) 模块搭建右侧的表达式r*x*(1 - x/K)。将表达式的结果连接到Integrator的输入口。设置仿真时间、初始条件在Integrator模块中设置然后运行。与 MATLAB 工作空间交互这是最实用的技巧。在 Simulink 模型中使用From Workspace模块读入 MATLAB 变量作为输入使用To Workspace模块将仿真结果输出到 MATLAB 工作空间。这样你就可以用 MATLAB 脚本预处理数据、批量运行仿真、并后处理结果。联合仿真搜索热词中提到了 Adams机械动力学软件与 MATLAB 的联合仿真。这通常通过S-Function(系统函数) 实现。Adams 可以生成一个描述机械系统的M文件 S-Function将其导入 Simulink 作为一个模块。Simulink 负责控制算法部分两者通过定义好的输入输出接口进行数据交换实现机电一体化系统的协同仿真。这个过程配置较为复杂需要确保两个软件的版本兼容并正确设置求解器参数。4.2 随机过程模拟醉汉游走模型随机游走是模拟许多自然和社会现象的经典模型如股票价格、粒子扩散等。% 一维简单随机游走醉汉游走 num_steps 1000; % 总步数 num_walkers 5; % 模拟多个醉汉 paths zeros(num_walkers, num_steps1); % 记录路径 for w 1:num_walkers position 0; paths(w, 1) position; for step 1:num_steps % 每一步以相等概率向左或向右移动一个单位 move randi([0, 1]) * 2 - 1; % 生成 -1 或 1 position position move; paths(w, step1) position; end end % 可视化 figure; hold on; for w 1:num_walkers plot(0:num_steps, paths(w, :)); end xlabel(步数); ylabel(位置); title(一维随机游走模拟); grid on; hold off; % 分析计算最终位置的均值和方差 final_positions paths(:, end); fprintf(最终位置均值: %.2f\n, mean(final_positions)); fprintf(最终位置方差: %.2f\n, var(final_positions)); % 理论上均值应为0方差约等于步数。这个模型可以轻松扩展变成二维游走每次随机选择东、南、西、北加入偏向性例如向右走的概率更大或者让步长服从某种分布如正态分布。在建模中它常用于作为更复杂模型的基准或组成部分例如在模拟搜索行为、市场波动时。4.3 自定义算法实现隐式 QR 方法求特征值当需要深入理解算法原理或者内置函数eig无法满足特殊需求时自己实现算法就很有必要。隐式 QR 方法是计算中小型矩阵全部特征值的稳定方法。function eigenvalues my_implicit_qr(A, max_iter, tol) % 简化的隐式QR算法实现带实用主义简化非工业级 % 输入方阵A最大迭代次数容忍度 % 输出特征值向量 n size(A, 1); H hess(A); % 先上Hessenberg化大幅减少计算量 eigenvalues zeros(n, 1); iter 0; while n 1 iter max_iter % 实用技巧检查次对角线元素是否可忽略 for i n:-1:2 if abs(H(i, i-1)) tol * (abs(H(i-1, i-1)) abs(H(i, i))) H(i, i-1) 0; % 发现可分离的子问题 if i n % 右下角1x1块得到一个特征值 eigenvalues(n) H(n, n); n n - 1; H H(1:n, 1:n); else % 分离为两个更小的子问题这里简化处理递归调用实际应迭代处理 % 此处为演示我们跳出循环用内置eig处理剩余部分 remaining_vals eig(H(1:n, 1:n)); eigenvalues(1:length(remaining_vals)) remaining_vals; return; end break; end end % 单步隐式QR位移使用Wilkinson位移以加速收敛 mu wilkinson_shift(H(n-1:n, n-1:n)); [Q, R] qr(H(1:n, 1:n) - mu * eye(n)); H(1:n, 1:n) R * Q mu * eye(n); iter iter 1; end if n 1 eigenvalues(1) H(1,1); end % 清理接近零的虚部 eigenvalues real(eigenvalues) 1i * (abs(imag(eigenvalues)) tol) .* imag(eigenvalues); end function mu wilkinson_shift(B) % 计算2x2矩阵的Wilkinson位移 a B(1,1); b B(1,2); d B(2,2); delta (a - d) / 2; mu d - sign(delta) * b^2 / (abs(delta) sqrt(delta^2 b^2)); end重要提醒自己实现数值线性代数算法主要是为了教学和理解。在真正的建模竞赛或工程中永远优先使用 MATLAB 内置的、经过高度优化的函数如eig,svd,chol等。它们更快、更稳定、更准确。自己写的代码可能在某些边界条件下失效比如几乎奇异矩阵。这里的实现是一个高度简化的教学版本重点展示了算法框架、Hessenberg 化降维、Wilkinson 位移加速和 deflation收缩技巧这些核心思想。5. 效率优化与调试让代码跑得更快更稳在建模中时间就是生命。一个需要跑一小时的脚本和另一个功能相同但只需一分钟的脚本有本质区别。5.1 向量化操作告别缓慢的 for 循环这是 MATLAB 性能提升的第一法则。尽可能使用矩阵和向量运算代替循环。% 低效做法计算两个向量对应元素的某种函数值例如距离计算 n 1e6; x randn(n, 1); y randn(n, 1); result_slow zeros(n, 1); tic; for i 1:n result_slow(i) sqrt(x(i)^2 y(i)^2); end time_slow toc; % 高效做法向量化 tic; result_fast sqrt(x.^2 y.^2); % 点乘 .^ 和 .* 是关键 time_fast toc; fprintf(循环耗时: %.4f 秒\n, time_slow); fprintf(向量化耗时: %.4f 秒\n, time_fast); fprintf(加速比: %.2f 倍\n, time_slow / time_fast);对于更复杂的多重循环嵌套可以考虑使用meshgrid或ndgrid生成网格然后一次性进行向量运算。但要注意这会瞬间创建巨大的矩阵可能耗尽内存。需要在向量化和内存消耗之间权衡。5.2 预分配数组避免动态增长在循环中不断向一个数组追加元素例如result [result; new_value]是性能杀手因为 MATLAB 每次都需要寻找新的连续内存并复制整个数组。% 糟糕的做法 result []; for k 1:10000 result [result; some_calculation(k)]; % 每次循环都重新分配内存 end % 正确的做法 n 10000; result zeros(n, 1); % 预先分配好内存 for k 1:n result(k) some_calculation(k); % 直接赋值 end使用zeros,ones,NaN,inf等函数进行预分配。如果不知道最终大小可以高估一个上限计算完后再截断。5.3 函数句柄与匿名函数提升代码灵活性函数句柄允许你将函数作为参数传递这在优化、积分和绘图时非常有用。% 定义目标函数用于 fminsearch rosenbrock (x) 100*(x(2)-x(1)^2)^2 (1-x(1))^2; x0 [-1.2, 1]; [x_opt, fval] fminsearch(rosenbrock, x0); % 在数组函数中使用如 arrayfun, cellfun但需谨慎有时不如直接向量化快 x_cell num2cell(1:10); squares cellfun((a) a^2, x_cell); % 对元胞数组每个元素平方 % 创建带参数的函数句柄 a 2; b 3; linear_func (x) a*x b; % 捕获了当前工作空间的 a 和 b y linear_func(5); % y 2*5313匿名函数非常简洁但过于复杂的逻辑不适合写在里面会影响可读性。对于复杂的、需要复用的功能还是应该写在独立的.m函数文件中。5.4 调试与错误排查从报错信息中快速定位MATLAB 的错误信息尤其是新版已经相当友好了。关键是要学会看。“函数或变量无法识别”比如热词中的‘deltalin’。首先检查拼写。然后检查该函数所在的文件夹是否在 MATLAB 的搜索路径中。使用which deltalin命令查看 MATLAB 是否能找到它。如果是一个自定义函数确保它的文件名和函数名一致并且当前目录或路径正确。矩阵维度不匹配这是最常见的错误之一。出错时MATLAB 会告诉你操作涉及矩阵的具体维度。立刻使用whos或在工作区查看相关变量的size。确保乘法 (*和.*)、加法、连接等操作的维度是兼容的。使用断点和调试模式在代码行号左侧点击设置断点红色圆点。运行代码时执行到断点处会暂停。此时你可以将鼠标悬停在变量上查看其当前值也可以在命令窗口直接输入变量名查看。使用F10单步执行F11步入函数F5继续运行。这是定位逻辑错误最强大的工具。try-catch块在可能出错的代码段外围加上try-catch可以捕获错误并执行备用方案避免程序完全崩溃。try data load(unreliable_data_file.mat); result complex_analysis(data); catch ME % ME 是一个包含错误信息的结构体 warning(加载或分析数据失败: %s, ME.message); result []; % 返回空值或默认值 % 可以记录错误日志fprintf(log_file, %s: %s\n, datetime, ME.message); endtic/toc和性能分析器用tic和toc包围你觉得慢的代码段来计时。对于更全面的性能分析使用profile工具。在命令窗口输入profile on运行你的代码然后输入profile viewer。它会生成一个详细的报告告诉你每行代码的执行时间和调用次数精准定位性能瓶颈。6. 可视化与成果输出让图表和论文说话好的可视化能让你和评委快速理解模型结果。而规范的输出则是最终成果的保障。6.1 高级绘图技巧超越 plot双Y轴图使用yyaxis。x linspace(0, 10); y1 sin(x); y2 exp(x/5); figure; yyaxis left; plot(x, y1, b-); ylabel(sin(x)); yyaxis right; plot(x, y2, r--); ylabel(exp(x/5)); xlabel(x); title(双Y轴示例);子图灵活排列使用tiledlayout替代旧的subplot它能更好地控制间距和标题。figure; t tiledlayout(2, 2); % 2行2列 nexttile; plot(rand(10,1)); title(图1); nexttile; scatter(rand(10,1), rand(10,1)); title(图2); nexttile(3, [1, 2]); % 跨一列占两列宽度 surf(peaks); title(跨列子图); xlabel(t, 共同的X轴标签); ylabel(t, 共同的Y轴标签); title(t, 整体标题);颜色映射与数据映射使用colormap来根据数据值着色这在显示矩阵、三维曲面时非常直观。data randn(20, 20); imagesc(data); colorbar; % 显示颜色条 colormap(jet); % 或 hot, cool, parula (默认) axis image; % 使坐标轴比例相等6.2 导出高质量图片用于论文的 EPS/PDF搜索热词里提到了导出 EPS。对于学术出版矢量格式EPS, PDF是首选因为它们无限缩放不失真。figure; plot(1:10, rand(1,10), o-, LineWidth, 2, MarkerSize, 8); xlabel(Index); ylabel(Value); title(示例图); grid on; % 方法1直接使用 exportgraphics (R2020a及以上推荐) exportgraphics(gcf, my_plot.pdf, ContentType, vector); % 导出为矢量PDF % exportgraphics(gcf, my_plot.eps, ContentType, vector); % 导出为EPS % 方法2使用 print 函数旧版本通用 set(gcf, Renderer, painters); % 确保使用矢量渲染器 print(-depsc, -tiff, -r300, my_plot.eps); % 导出EPS嵌入TIFF预览分辨率300dpi % print(-dpdf, -bestfit, my_plot.pdf); % 导出PDF自动适应页面 % 方法3调整画布大小再导出 fig gcf; fig.PaperPositionMode auto; % 使用屏幕尺寸 fig_pos fig.PaperPosition; fig.PaperSize [fig_pos(3) fig_pos(4)]; % 设置纸张大小与图形一致 print(fig, my_plot.pdf, -dpdf, -bestfit);注意如果图形中包含大量数据点如十万个散点导出为矢量文件可能会非常大。此时可以导出为高分辨率位图如PNG,‘-r600’或者使用exportgraphics并指定‘Resolution’参数。在导出前务必检查图形中所有文字标签、图例、标题的字体和大小是否合适因为一旦导出为图片就很难再修改。6.3 生成报告与文档发布你的代码MATLAB 的“发布”Publish功能可以将脚本.m文件转换为格式良好的报告HTML, Word, PDF 等包含代码、运行结果和图表。在编辑器里写好你的脚本用%%创建分节单元格。在每个单元格中先写注释Markdown格式再写代码。点击“发布”选项卡上的按钮选择输出格式。MATLAB 会按顺序执行每个单元格并将输出包括文本、变量值、图形嵌入到生成的文档中。这对于撰写建模论文的附录、制作实验报告或分享可复现的研究结果极其有用。它能确保文档中的结果与代码完全对应。