1. 项目概述从“种群竞争”到数学建模的实战跨越最近在整理往年带学生做数学建模竞赛的资料发现“种群竞争模型”几乎是一个绕不开的经典课题。无论是国赛、美赛还是各类校赛只要题目涉及到生态、资源分配、市场博弈甚至社交媒体上的信息传播这个模型的影子就可能出现。2023年的赛题中虽然没有直接以“两种生物争地盘”这样直白的表述出现但其内核——两个或多个主体在有限资源下的动态博弈——被包装在了各种新颖的场景之下。很多新手队伍看到题目描述复杂就感觉无从下手其实剥开现象看本质很可能就是在考察你对经典竞争模型的理解、改造和应用能力。简单来说种群竞争模型描述的是两个或多个物种为了共同的生存资源如食物、空间而相互制约的动态过程。它的价值远不止于解释“兔子多了草不够吃”这样的生态学现象。在当下我们可以用它来分析互联网平台上不同内容类型的流量争夺可以用来模拟有限预算下几个营销渠道的效果博弈甚至可以用来理解城市发展中不同产业对人才和政策的竞争关系。这次我就结合2023年一些赛题的思路抛开纯理论推导重点聊聊怎么用MATLAB把这个模型从纸面公式变成能跑出结果、能解释问题的实用工具。无论你是正在备赛的学生还是对动态系统建模感兴趣的爱好者这篇内容都能给你一套可直接上手操作的“工具箱”。2. 模型内核解析洛特卡-沃尔泰拉方程到底在说什么提到种群竞争就不得不提洛特卡-沃尔泰拉方程Lotka-Volterra competition equations。别被这个名字吓到它本质上就是一组考虑了“内部摩擦”和“外部压力”的加强版指数增长方程。2.1 核心方程拆解每个参数的意义我们先看最经典的两物种竞争模型dN1/dt r1 * N1 * (1 - (N1 α12 * N2) / K1) dN2/dt r2 * N2 * (1 - (N2 α21 * N1) / K2)这里每个符号都不是天书我们一个一个拆开看N1, N2: 物种1和物种2在时间t的数量。这是我们要解的核心变量。r1, r2: 物种的内禀增长率。可以理解为在理想条件下食物无限、空间无限、没有竞争对手种群增长的“最大本事”。比如细菌的r值就远大于大象。K1, K2: 环境容纳量。这是模型的关键约束代表在只有该物种自身竞争种内竞争的情况下环境能支撑的最大种群数量。它由资源总量决定。α12 和 α21:竞争系数。这是整个模型的灵魂也是最容易用错的地方。α12表示物种2对物种1的竞争效应。具体来说一个物种2的个体对物种1造成的资源压力相当于α12个物种1的个体。如果α120.5就意味着每增加2个物种2对物种1资源消耗的影响相当于增加1个物种1。α21同理表示物种1对物种2的竞争效应。注意竞争系数α不是“谁厉害”的简单评分。它是一个折算系数用于将不同物种的个体数量统一到对同一资源消耗的“标准单位”下。理解成“兑换率”更准确。2.2 四种结局的直观理解两个物种竞争长期来看无非四种结局这完全由竞争系数和环境容纳量的关系决定物种1胜出物种2灭绝当α12 K1/K2且α21 K2/K1时发生。意味着物种2对1的抑制很强α12大但物种1对2的抑制很弱α21小同时物种1自己的承载力K1相对较高。结果是物种2被彻底排挤。物种2胜出物种1灭绝当α12 K1/K2且α21 K2/K1时发生。情况与上一种相反。稳定共存当α12 K1/K2且α21 K2/K1时发生。这是最有趣的局面意味着双方对彼此的抑制力都相对较弱弱到无法将对方彻底排除。两者会达到一个非零的平衡点在这个点上两个种群的数量保持稳定。不稳定共存胜负取决于初始数量当α12 K1/K2且α21 K2/K1时发生。这意味着双方对彼此的抑制力都非常强形成了“既生瑜何生亮”的局面。最终谁能存活完全看起步时谁的数量多初始条件赢家通吃。这个结局在数学上存在一个不稳定的平衡点现实中稍有扰动就会打破。实操心得很多同学在建模时拍脑袋给α赋值比如觉得A物种强就设α2B物种弱就设α0.5这很容易导致模型行为与常识背离。一定要根据题目中隐含的“资源等价关系”来估算α。例如如果题目说“一个单位的企业甲消耗的资源是企业的两倍”那么在模拟两者竞争同一市场资源时企业乙对企业甲的竞争系数α_乙甲就可能接近2因为一个乙相当于两个甲的资源消耗力。3. MATLAB实现全流程从零搭建你的仿真系统理论懂了关键在实现。下面我们一步步在MATLAB里搭建一个可运行、可调节、可视化的竞争模型。3.1 环境准备与模型函数定义首先我们创建一个名为lotka_volterra_competition.m的函数文件。这样做的好处是参数清晰易于调试和重复调用。function dN lotka_volterra_competition(t, N, params) % 洛特卡-沃尔泰拉竞争模型微分方程 % 输入: % t: 时间 (未直接使用但ODE求解器要求此参数) % N: 当前种群数量向量 [N1; N2] % params: 结构体包含所有参数 r1, r2, K1, K2, alpha12, alpha21 % 输出: % dN: 微分方程结果 [dN1/dt; dN2/dt] % 解包参数 r1 params.r1; r2 params.r2; K1 params.K1; K2 params.K2; alpha12 params.alpha12; alpha21 params.alpha21; % 解包当前种群数量 N1 N(1); N2 N(2); % 计算微分方程 dN1_dt r1 * N1 * (1 - (N1 alpha12 * N2) / K1); dN2_dt r2 * N2 * (1 - (N2 alpha21 * N1) / K2); % 输出 dN [dN1_dt; dN2_dt]; end为什么用结构体params传递参数直接传递6个参数当然可以但在主程序中反复调整参数做测试时结构体让代码更整洁。你只需要修改params里的字段而不需要改动函数调用接口避免了参数顺序错乱的风险。3.2 主程序编写求解与可视化接下来我们写一个主脚本main_competition.m来调用这个函数完成求解和画图。% 清除环境 clear; close all; clc; % 1. 设置模型参数这里以“不稳定共存”为例制造悬念 params.r1 0.5; % 物种1增长率 params.r2 0.5; % 物种2增长率 params.K1 1000; % 物种1环境承载力 params.K2 800; % 物种2环境承载力 params.alpha12 1.2; % 物种2对1的竞争系数 params.alpha21 1.1; % 物种1对2的竞争系数 % 根据公式计算临界条件 cond1 params.alpha12 params.K1 / params.K2; % 应大于 1.25 cond2 params.alpha21 params.K2 / params.K1; % 应大于 0.8 fprintf(alpha12 K1/K2 ? %d (%.2f %.2f)\n, cond1, params.alpha12, params.K1/params.K2); fprintf(alpha21 K2/K1 ? %d (%.2f %.2f)\n, cond2, params.alpha21, params.K2/params.K1); if cond1 cond2 disp(- 参数设定为“不稳定共存”类型结局取决于初始数量。); end % 2. 设置初始条件和时间范围 N0 [100; 100]; % 初始数量 [N1; N2]试试改成[150; 50]看结果如何 tspan [0 50]; % 模拟时间范围 0到50个单位时间 % 3. 使用ODE45求解微分方程 % ‘(t,N)’创建了一个匿名函数将固定的params传递给我们的模型函数 [t, N] ode45((t,N) lotka_volterra_competition(t, N, params), tspan, N0); % 4. 绘制种群数量随时间变化图 figure(Position, [100, 100, 1200, 400]); % 设置图形窗口大小 subplot(1, 2, 1); plot(t, N(:, 1), b-, LineWidth, 2); hold on; plot(t, N(:, 2), r--, LineWidth, 2); grid on; xlabel(时间); ylabel(种群数量); title(种群数量动态变化); legend(物种1 (N1), 物种2 (N2), Location, best); hold off; % 5. 绘制相平面图相位肖像 subplot(1, 2, 2); plot(N(:,1), N(:,2), k-, LineWidth, 1.5); hold on; scatter(N0(1), N0(2), 100, go, filled); % 标记起点 scatter(N(end,1), N(end,2), 100, rs, filled); % 标记终点 xlabel(物种1数量 N1); ylabel(物种2数量 N2); title(相平面图 (N1-N2 Phase Portrait)); legend(演化轨迹, 起点, 终点, Location, best); grid on; hold off; % 6. 计算并显示平衡点令微分方程为0求解 % 平衡点方程组 % N1 alpha12*N2 K1 % N2 alpha21*N1 K2 A [1, params.alpha12; params.alpha21, 1]; B [params.K1; params.K2]; equilibrium_point A \ B; % 解线性方程组 fprintf(\n理论平衡点不一定是稳定点: N1* %.2f, N2* %.2f\n, equilibrium_point(1), equilibrium_point(2));实操要点ODE45的选择对于这种非刚性的常微分方程组ode45是首选它结合了四阶和五阶龙格-库塔法在精度和速度上平衡得很好。如果模型变得非常复杂或出现“刚性”问题某些变量变化极快才需要考虑ode15s。相平面图的价值时间序列图告诉我们“怎么变”相平面图则揭示了两个变量之间的内在关系。轨迹线直观展示了系统演化的路径起点和终点的标记让你一眼看清趋势。在论文中放入这样一张图能极大提升分析的专业性。平衡点计算代码中求解的平衡点是数学上的静态解。但它的稳定性需要进一步通过雅可比矩阵特征值来判断。如果特征值实部都小于0则是稳定平衡点吸引子如果有大于0的则是不稳定点鞍点或排斥子。3.3 参数敏感性分析与情景模拟真正的建模高手不会只满足于跑通一个案例。我们要测试不同参数下系统的行为这就是参数敏感性分析也是论文中“模型分析”部分的干货。% 参数敏感性分析示例改变竞争系数alpha12观察结局 figure(Position, [100, 100, 1200, 800]); alpha12_values [0.5, 1.0, 1.2, 1.5]; % 测试四个不同的alpha12值 N0_fixed [100; 100]; for i 1:length(alpha12_values) params_test params; % 复制基础参数 params_test.alpha12 alpha12_values(i); % 求解 [t_test, N_test] ode45((t,N) lotka_volterra_competition(t, N, params_test), tspan, N0_fixed); % 绘图 subplot(2, 2, i); plot(t_test, N_test(:,1), b-, LineWidth, 1.5); hold on; plot(t_test, N_test(:,2), r--, LineWidth, 1.5); grid on; xlabel(时间); ylabel(数量); title(sprintf(\\alpha_{12} %.1f, alpha12_values(i))); legend(N1, N2, Location, best); hold off; % 判断结局 if N_test(end, 2) 1e-2 fprintf(当alpha12%.1f时物种2灭绝。\n, alpha12_values(i)); elseif N_test(end, 1) 1e-2 fprintf(当alpha12%.1f时物种1灭绝。\n, alpha12_values(i)); else fprintf(当alpha12%.1f时两物种共存。\n, alpha12_values(i)); end end sgtitle(不同竞争系数α12下的种群动态); % 总标题通过这个循环你可以清晰地看到仅仅改变一个竞争系数系统的长期命运就可能从共存变为一方灭绝。在论文中你可以用类似的代码批量运行制作出展示参数变化如何影响结局的图表这比干巴巴的文字论述有力得多。4. 从经典模型到赛题应用2023年思路延伸经典模型是骨架赛题应用需要血肉。我们看看如何把“种群竞争”的思维应用到更广泛的场景。4.1 场景拓展不止于生物学市场竞争模型两个公司N1, N2争夺有限的市场总需求K。增长率r可以理解为公司的市场扩张能力或营销效率。竞争系数α可以理解为产品的替代性强度。如果产品高度同质化如两家卖同样矿泉水α值会接近1甚至大于1容易导致价格战和一方出局不稳定共存。如果产品有差异化如一家卖咖啡一家卖茶饮α值会小于1可能实现稳定共存细分市场。社交媒体信息传播两种不同类型的信息如新闻vs.谣言争夺用户的有限注意力K。增长率r代表信息的“爆点”潜力或传播速率。竞争系数α代表一种信息出现对另一种信息关注度的压制程度。模型可以模拟在热点事件中真实信息和虚假信息的博弈过程。城市土地规划商业用地N1与绿化用地N2争夺城市有限的空间资源K。增长率r可能代表不同利益集团的推动力。竞争系数α则反映了政策导向例如是否允许商业用地侵占绿化指标。通过调节α可以模拟不同政策下的土地分配结局。4.2 模型改良让模型更“聪明”经典LV模型假设是线性的、确定性的。实际赛题中我们需要让它更贴合现实。加入随机项随机微分方程环境不是一成不变的。可以在微分方程中加入随机噪声项模拟环境波动、随机事件的影响。这能让你分析结果的稳健性。% 简化的Euler-Maruyama方法思路示意 dN1 r1*N1*(1-(N1alpha12*N2)/K1)*dt sigma1*N1*randn*sqrt(dt);考虑时变参数增长率r或承载力K可能随时间周期性变化如季节性资源或受外部政策影响。可以将r1、K1定义为时间t的函数r1(t)。引入第三种群或捕食者构建一个更复杂的食物网或竞争网络。例如物种1和2竞争同时都被物种3捕食。这会形成更丰富的动力学行为如周期性震荡。空间显式模型经典模型是“一锅粥”式的假设个体均匀混合。可以引入元胞自动机或反应扩散方程考虑种群在空间上的分布和扩散研究竞争与空间格局的相互影响。实操心得在竞赛论文中不要一上来就摆出最复杂的改良模型。经典模型 - 参数拟合/敏感性分析 - 指出经典模型的不足 - 提出你的改良模型 - 对比分析改良效果这样的行文逻辑更清晰也更能体现你的思考深度。5. 常见问题与调试技巧实录在实际编程和备赛过程中你肯定会遇到各种问题。下面是我和学生们踩过的坑以及解决办法。5.1 模型求解与数值问题问题现象可能原因排查与解决技巧种群数量出现负值1. 时间步长太大2. 竞争过于激烈数值求解器“冲过头”。1. 使用ode45的options参数限制最大步长options odeset(MaxStep, 0.1);然后传入求解器。2. 在模型函数中加入判断if N1 0; dN1_dt 0; end强制数量不为负。但这会改变模型性质需在论文中说明。解算速度非常慢1. 参数设置导致方程“刚性”部分变量变化极快2. 模拟时间tspan过长。1. 尝试换用刚性求解器ode15s或ode23s。2. 检查参数数量级是否差异巨大如r0.01, K10000尝试归一化处理将N除以K变成无量纲的相对数量。平衡点计算为NaN或Inf竞争系数导致系数矩阵A奇异不可逆即1 - alpha12*alpha21 0。检查参数合理性。在生物意义上alpha12*alpha211是临界情况通常应避免。调整参数或直接分析此临界状态的意义。5.2 参数设定与结果解释问题跑出来的结果两个物种都快速增长然后维持在高位几乎没有竞争迹象。排查检查竞争系数α是否设得太小比如都设为0.1或者环境承载力K设得过大。这相当于资源极度丰富竞争效应不明显。合理设定K和α的比例关系是关键。技巧先做量纲分析。假设N1和N2代表个体数K11000那么初始值N0设为10或100是合理的。如果N01K11000增长空间太大图形前期会有一段很长的平缓期。根据r值如0.5/时间单位估算种群翻倍时间大约是1.4个单位时间以此调整tspan让演化过程完整展示。5.3 论文图表美化与表达图表除了最基本的时间序列图和相图可以绘制参数空间图以alpha12和alpha21为坐标轴划分出四个结局区域并用散点标出你模拟的参数点落在哪个区域。平衡点稳定性图计算雅可比矩阵特征值用颜色图表示平衡点稳定性随参数的变化。动态演示使用for循环和drawnow命令制作种群数量动态变化的动画在答辩或视频摘要中非常出彩。表达在论文中描述模型时避免直接贴大段代码。应该用公式、流程图和文字描述算法思路。将核心代码如模型函数、关键求解步骤放在附录。解释参数时务必说明其实际意义和取值依据是参考了文献还是根据题目数据估算的。最后我个人最深的体会是种群竞争模型就像一个“数学积木”。经典方程是基础块真正的挑战在于你如何根据具体问题挑选、修改、拼接这些积木甚至自己创造新的积木。在2023年的很多赛题里胜出的队伍并非使用了多么高深的算法而是因为他们把“竞争”这个概念理解透了并用一个恰当的、自洽的模型将其清晰地表达了出来。MATLAB是实现这个想法的强大工具但比工具更重要的是你对问题本质的洞察力和将现实抽象为数学模型的能力。多练、多改、多思考不同的参数组合会产生什么结果你就能逐渐培养出这种“建模直觉”。下次再遇到“竞争”、“博弈”、“此消彼长”这类关键词时希望你能第一时间想到这个模型并知道如何让它为你所用。