我第一次认真研究Lotka-Volterra模型是因为看到一组北美山猫和雪兔的皮毛收购数据两种动物每隔9到10年就同步振荡一次雪兔多了山猫跟着多雪兔崩了山猫也跟着崩。这种周期波动非常规律以至于你很难相信背后没有一个简单数学机制。后来把方程写出来才发现两个物种的此消彼长本质上就是一组很干净的常微分方程而MATLAB恰好是把它跑起来、画出来的最快工具。这篇内容适合所有刚接触种群生态模型的人不管你是生态学学生、数学建模选手还是想把LV模型作为练习ODE求解的MATLAB进阶玩家。1. 模型背后的生态直觉为什么捕食者与猎物会此消彼长1.1 从真实观测数据到数学抽象LV模型最早是Lotka在1925年研究化学反应时提出第二年Volterra用它解释亚得里亚海捕食鱼类占比的变化。这两个人从不同学科出发最后落在同一个方程上本身就说明数学结构并不挑应用场景。模型的核心假设是在没有捕食者时猎物以固定速率增长没有猎物时捕食者以固定速率死亡。当两者相遇捕食者吃掉猎物把一部分生物量转化为自己的繁殖能量。于是写出这样一组方程dx/dt alpha*x - beta*x*y dy/dt delta*x*y - gamma*y这里的x代表猎物数量y代表捕食者数量。alpha是猎物的自然增长率beta是捕食效率delta是捕食者将猎物转化为后代的效率gamma是捕食者的自然死亡率。四个参数全是正数模型不考虑环境承载、不考虑年龄结构、不考虑空间差异把所有复杂生态关系压缩成了一条相遇率假设——捕食者与猎物的接触次数正比于两者的乘积x*y。这个乘积是整组方程的灵魂。它假设捕食者和猎物在空间里均匀混合、随机相遇类似理想气体分子碰撞。野外当然不是这样但作为第一步近似它抓住了周期振荡的根源猎物多导致捕食者吃得饱、繁殖快捕食者增加反过来压制猎物猎物减少后捕食者饿死捕食者减少又给了猎物恢复空间。负反馈回路天然成立周期解就自然涌现。1.2 参数变化会带来什么后果我刚学这个模型时习惯性忽略了参数取值的重要性结果第一次仿真猎物数量冲到天文数字捕食者又低到几乎灭绝完全不像真实观测。后来才意识到参数不仅决定振幅还决定振荡周期和平衡点位置。基础LV模型的非零平衡点很容易算令dx/dt0且dy/dt0得到x*gamma/deltay*alpha/beta。这个结果很有力——平衡猎物密度只由捕食者死亡率和转化效率决定平衡捕食者密度只由猎物增长率和捕食效率决定。想让猎物平衡值降低可以提高捕食者转化效率或者想办法降低捕食者死亡率这从管理角度看非常有指导意义。而振荡周期约为2pi/sqrt(alphagamma)和初始种群大小没有关系这也是模型一个重要特征只要四个参数不变无论从什么初始值出发系统都会以相同周期绕圈。理解这一点后面做MATLAB仿真时你就知道该把注意力放在哪。很多初学者只顾着调初始值看曲线变化其实初始值只是决定绕圈路径半径参数才是决定整个系统行为的关键。2. MATLAB中从零到一用ODE45把基础LV模型跑起来2.1 写出符合MATLAB语法的一阶ODE方程组MATLAB求解这类初值问题思路是把高阶方程或方程组写成y f(t, y)的标准形式。LV模型本来就一阶常微分方程组直接用匿名函数或者单独的函数文件喂给ode45即可。% lv_basic.m alpha 1.1; % 猎物自然增长率 beta 0.4; % 捕食效率 delta 0.1; % 捕食者转化效率 gamma 0.4; % 捕食者死亡率 f (t, z) [alpha*z(1) - beta*z(1)*z(2); delta*z(1)*z(2) - gamma*z(2)]; tspan [0 80]; z0 [20; 8]; [t, z] ode45(f, tspan, z0); figure; plot(t, z(:,1), b-, LineWidth, 1.8); hold on; plot(t, z(:,2), r--, LineWidth, 1.8); legend(猎物 N(t), 捕食者 P(t)); xlabel(时间 t); ylabel(种群数量); grid on; figure; plot(z(:,1), z(:,2), k-, LineWidth, 1.5); xlabel(猎物 N); ylabel(捕食者 P); axis equal; grid on;这里向量z的第一位是猎物第二位是捕食者。匿名函数虽然方便但参数一多建议改成单独函数文件因为后面加Logistic项、Holling功能响应时表达式会膨胀得很厉害匿名函数嵌套匿名函数会非常难维护。我通常写一个lv_fun.m把参数放到调用处用deal结构体传递这样扫描参数时不用反复修改函数体。2.2 运行与验证能量守恒量的检验基础LV模型有个隐藏性质它存在一个守恒量。沿着任意一条相轨线下面这个组合保持恒定V delta*x - gamma*ln(x) beta*y - alpha*ln(y)这有点像物理学里的机械能守恒。你可以先用初始值计算V0再在ode45输出后计算每个时间点的V检查最大偏差。数值方法毕竟有截断误差V会有微小漂移但偏差应保持在很小的量级。这个检验是我强烈建议保留的它是判断你是否把方程写错的最快手段——如果V明显漂移说明要么符号写错要么求解精度太低。C0 delta*z0(1) - gamma*log(z0(1)) beta*z0(2) - alpha*log(z0(2)); C delta*z(:,1) - gamma*log(z(:,1)) beta*z(:,2) - alpha*log(z(:,2)); maxDev max(abs(C - C0)); disp(maxDev);我第一次算出来最大偏差只有1e-10量级当时还很惊讶后来把RelTol降低到默认值偏差就升到1e-5左右。这个指标可以作为数值精度的晴雨表。2.3 画好时间序列图和相图的几个细节时间序列图就是把两个物种数量分别画在纵轴上、时间画在横轴上直观看相位差。相图则以猎物为横轴、捕食者为纵轴每个初始条件对应一圈闭合曲线。同一参数下不同初始条件对应不同大小的环互不相交。画相图有个细节常被忽略plot(z(:,1), z(:,2))默认会连出闭合曲线但如果求解器在早期步长太大曲线会出现尖角甚至自交的假象。建议对基础模型把RelTol设到1e-8再用axis equal让横纵坐标比例一致否则圆会被拉伸成椭圆视觉上误以为轨线形状发生了变化。要叠加多条相轨线可以用循环对不同初始值分别调用ode45再用hold on叠加。这样生成的一簇闭合曲线能够直观展示模型在所有初始条件下的运动模式也是论文里展示LV模型最标准的配图。3. 让模型更接近野外Logistic自限与Holling功能响应3.1 猎物无限增长的修正基础模型里即使没有捕食者猎物也按指数增长这在野外几乎不可能。资源有限、空间有限、疾病传播等都会让种群受到密度制约。标准修正是把猎物增长率从常数alpha改成Logistic形式dx/dt r*x*(1 - x/K) - beta*x*y dy/dt delta*x*y - gamma*yr是猎物内禀增长率K是环境容纳量。这样猎物单独存在时趋于K不会爆炸捕食者加入后猎物会在K以下振荡而不是围绕一个固定平衡点做等幅循环。实际计算后你会发现系统通常收敛到稳定平衡点或衰减振荡经典的中性周期消失了。这一改动看似简单却把模型从保守系统变成了耗散系统。这意味着初始条件的影响会随时间消失长期行为主要由参数决定而不是由起点决定。处理真实观测数据时耗散系统通常更合理因为野外种群很少呈现永恒等幅振荡。3.2 Holling II型功能响应的处理基础模型假设捕食者吃猎物吃的数量与猎物密度成正比——这在猎物极稀疏时合理但在猎物密度高时捕食者总有个处理食物的时间不可能无限吃下去。Holling给出了三类功能响应函数其中II型最常用捕食率 a*x / (1 a*h*x)a是攻击率h是处理单个猎物所需时间。当x很小时该项近似a*x当x很大时趋于1/h即捕食率被处理时间饱和。这个函数图像像一把逐渐变平的弯刀本质上和化学反应里的Michaelis-Menten方程是同一套数学。加入Holling II型后完整模型变成dx/dt r*x*(1 - x/K) - a*x*y/(1 a*h*x) dy/dt e*a*x*y/(1 a*h*x) - m*y其中e是捕食者把被捕食猎物转化为新捕食者的效率m是捕食者死亡率。这里我把delta拆成了e*a因为转化率和攻击率在机制上应该是独立的。r 1.2; K 100; a 0.6; h 0.2; e 0.4; m 0.3; f (t, z) [r*z(1)*(1 - z(1)/K) - a*z(1)*z(2)/(1 a*h*z(1)); e*a*z(1)*z(2)/(1 a*h*z(1)) - m*z(2)]; [t, z] ode45(f, [0 200], [40; 4]); figure; plot(t, z(:,1), b-, LineWidth, 1.8); hold on; plot(t, z(:,2), r--, LineWidth, 1.8); legend(猎物, 捕食者); xlabel(时间 t); ylabel(种群数量); grid on;3.3 加上这些项之后平衡点发生了什么Holling II型模型有多个平衡点。平凡平衡点(0,0)不稳定无捕食者平衡点(K,0)可能稳定也可能不稳定非零平衡点需要同时满足x* m / (a*(e - m*h)) y* (r/a)*(1 - x*/K)*(1 a*h*x*)这个公式透露的生态含义非常重要。捕食者平衡密度x与猎物模型参数r和K完全无关只取决于捕食者自身的a、e、m、h但y的大小则同时受到猎物承载力的影响。只有满足e mhx才是正的否则捕食者效率太低根本无法维持自身种群。这个不等式是物种共存的基本门槛。仿真前用这个条件判断参数组合是否合理能省掉大量无效调试。4. 进一步延伸比例依赖、时滞与季节扰动4.1 Arditi-Ginzburg比例依赖模型基础模型和Holling型功能响应都有一个共同弱点它们假设捕食者搜索到的猎物数量只和猎物绝对密度成正比而实际上捕食者之间会互相干扰捕食率往往更依赖于猎物与捕食者的比例。Arditi和Ginzburg在1989年提出比例依赖模型把功能响应改为捕食率 alpha*x / (y c*x)这里的c可以理解为捕食者之间干扰程度的倒数。当y远小于x时该项接近alpha/c捕食率接近常数描述的是猎物充足而捕食者互相干扰的场景当y远大于x时该项近似alpha*x/y描述的是捕食者太多、每个个体分到的猎物很少的场景形式上是只与猎物捕食者比例相关。MATLAB中实现无非是把分母写进去dg (t, z) [r*z(1)*(1 - z(1)/K) - alpha*z(1)*z(2)/(z(2) c*z(1)); beta*z(1)*z(2)/(z(2) c*z(1)) - m*z(2)];这种模型能产生更丰富的动力学行为包括Hopf分岔导致的极限环。扫几个初始值后你会发现相图不一定收敛到平衡点也可能收敛到同一个闭合环这就是自激振荡。大多数真实捕食系统的振荡不是中性周期而是极限环这也是比例依赖模型受到重视的原因。4.2 用dde23处理繁殖时滞真实捕食者吃掉猎物后不可能立刻转化为新个体。从进食到孕育、出生还有时间延迟这个延迟可能长达一个季节。连续时间模型中加入时滞方程就变成延迟微分方程(DDE)MATLAB用dde23求解。一个带时滞的改进模型可以写成dx/dt r*x(t)*(1 - x(t)/K) - a*x(t)*y(t)/(1 a*h*x(t)) dy/dt e*a*x(t-tau)*y(t-tau)/(1 a*h*x(t-tau)) - m*y(t)这里捕食者的繁殖项使用了tau时刻之前的猎物与捕食者数量代表捕食效率经过延迟tau后才反映到新个体上。tau 1.5; lags tau; history (t) [50; 5]; sol dde23(ddefun, lags, history, [0 200]); figure; plot(sol.x, sol.y(1,:), b-, LineWidth, 1.8); hold on; plot(sol.x, sol.y(2,:), r--, LineWidth, 1.8); xlabel(时间 t); ylabel(种群数量); grid on; function dz ddefun(t, z, Z) r 1.2; K 100; a 0.6; h 0.2; e 0.4; m 0.3; x_delay Z(1, 1); y_delay Z(2, 1); dz [r*z(1)*(1 - z(1)/K) - a*z(1)*z(2)/(1 a*h*z(1)); e*a*x_delay*y_delay/(1 a*h*x_delay) - m*z(2)]; endZ(:,1)表示过去时刻t-tau的状态变量矩阵。实测下来时滞会让系统更倾向于振荡tau较小时系统趋于稳定平衡点tau超过某个临界值后开始等幅振荡这就是典型的Hopf分岔。你想在论文里展示分岔现象可以扫几个tau绘制时间序列对比比直接讨论公式直观得多。4.3 把季节变化放入参数温带生态系统的很多参数并不是常数。猎物繁殖率随季节变化、捕食者死亡率在冬天上升这些都可以通过把参数改成时间t的周期函数来近似。最简单的做法是令r(t) r0*(1 epsilonsin(2pi*t/T))其中T是一年epsilon是季节波动幅度系数。在MATLAB里不需要修改求解流程把匿名函数里的alpha换成时间段相关表达式即可r_base 1.2; epsilon 0.3; T_year 12; % 以月为单位为例 r_fun (t) r_base * (1 epsilon * sin(2*pi*t/T_year)); f (t,z) [r_fun(t)*z(1)*(1-z(1)/K) - a*z(1)*z(2)/(1a*h*z(1)); e*a*z(1)*z(2)/(1a*h*z(1)) - m*z(2)];加入季节扰动后可能出现一个现象外部周期驱动与系统自身振荡频率发生共振导致振幅随着季节周期性变化甚至进入混沌。这种“驱动系统”的建模思路比单纯加项有意思得多很多生态教科书里的复杂模式都能用这个框架解释。实际操作中记得把时间跨度设足够长让系统先越过初始过渡段再截取稳态部分分析。5. 用数值实验做稳定性分析Jacobian矩阵与参数扫描5.1 平衡点与Jacobian矩阵前面对基础LV模型的分析有一个严格数学基础局部稳定性由平衡点处的Jacobian矩阵特征值决定。对于二维系统J [dF1/dx, dF1/dy; dF2/dx, dF2/dy]如果你不想手动求偏导MATLAB的符号工具箱很省力。先定义符号变量和表达式再用jacobian函数求矩阵然后代入平衡点数值eig求特征值。这套流程在改进模型中尤其有用因为表达式越来越长手动求导很容易漏项。对基础LV模型非零平衡点(x*,y*)(gamma/delta, alpha/beta)处的Jacobian特征值是纯虚数±isqrt(alphagamma)所以平衡点是中心型表现为中性周期振荡对应相图上一圈圈闭合曲线。Logistic修正后特征值实部变成负值平衡点变为稳定焦点种群振荡会逐渐衰减。这一个实部正负的变化就决定了系统是永久振荡还是趋向稳定是整个建模分析的分水岭。5.2 参数扫描的两种方式做灵敏度分析时最常见问题是单参数扫描结果画出来只有一条线看不出参数间的交互作用。我习惯做二维扫描把参数网格化对每个组合求解记录某个指标如均衡时猎物最小值、捕食者最大值再用imagesc或pcolor画热图。r_list linspace(0.6, 1.8, 40); m_list linspace(0.1, 0.6, 40); minPrey zeros(40, 40); for i 1:40 for j 1:40 r r_list(i); m m_list(j); f (t, z) [r*z(1)*(1-z(1)/K) - a*z(1)*z(2)/(1a*h*z(1)); e*a*z(1)*z(2)/(1a*h*z(1)) - m*z(2)]; [~, z] ode45(f, [0 300], [30; 5]); minPrey(i,j) min(z(:,1)); % 去掉过渡段再做min更准 end end figure; imagesc(m_list, r_list, minPrey); xlabel(捕食者死亡率 m); ylabel(猎物增长率 r); colorbar; title(猎物最小密度随参数变化);如果你有并行计算工具箱把外层for改成parfor一次扫描的时间能降到原来的几分之一。扫描结束后可以再加一个mask判断系统是否灭绝minPrey小于某个阈值时记为0这样能直观看到参数空间中物种存活的区域比只盯着单条时间序列判断强得多。6. 我在实操中踩过的坑与解决思路6.1 求解器不是万能的ode45适合大多数非刚性问题但改进模型加了Logistic项和Holling型功能响应后方程可能变刚。一个典型信号是ode45报错说无法满足积分容差或者计算速度明显变慢。这时不要盲目调低容差先用ode15s试试通常一切顺了。实测基础LV模型用ode45没问题但加了时间延迟饵料供给之类的强反馈后建议直接用ode15s。没有这个意识的话新手常被“ode45可以解所有ODE”这句话带偏实际它在刚性问题上的效率惨不忍睹。6.2 初值、残差与数值振荡初值设置不能拍脑袋。对基础LV模型如果你的初始捕食者数量为零捕食者永远不会出现系统退化为猎物指数增长解就会很大。这类退化解不是数值错误而是模型本身的动力学限制。写代码之前先花一分钟检查初始值是否都在正象限内。还有一次我遇到时间序列尾部出现微小跳动一开始以为模型错了后来发现是步长太大导致捕食者密度接近零时数值计算出负值被Logistic项或对数项放大。解决办法是把NonNegative选项打开opts odeset(RelTol,1e-8, NonNegative,[1 2]); [t, z] ode45(f, tspan, z0, opts);这个选项让解始终非负避免为负种群数量带来的荒谬结果。代价是极少数情况下求解器可能变慢但对种群模型来说非负约束几乎总是值得的。6.3 代码结构上的建议写长模型时最忌把所有参数堆在脚本开头然后在下文反复使用。一旦改参数忘记同步结果就会悄悄出错。我的习惯是把参数打包成结构体p.r 1.2; p.K 100; p.a 0.6; p.h 0.2; p.e 0.4; p.m 0.3;然后在函数文件里用p作为第二个参数传入。这样拷贝到其他脚本做参数扫描时非常清晰不存在闭包捕获旧参数值的坑。另外一个细节是给文件起名尽量用有意义的英文不要用中文文件名也不要起LV_model_MATLAB_final这种容易覆盖的版本名。我见过太多人在文件名上翻车最后分不清哪个脚本是更新后的版本。配合Git使用是更好的习惯至少用数字版本号。最后再分享一个我的个人体会Lotka-Volterra模型最迷人的地方不在于公式本身而在于它是一个可以不断往上加复杂度的框架。基础模型让人理解振荡的根源Logistic项让模型贴近现实Holling型功能响应带来饱和效应时滞和季节驱动则展现了分岔与混沌的可能性。每一步改动都能在MATLAB里几分钟内看到结果这种即时反馈是其他分析工具很难替代的。我的建议是先从最简模型跑通闭环验证守恒量画好相图再一项一项往上加复杂度每加一项都做一次稳定性分析。这样你既不会迷失在参数里也能清晰知道每种生物学机制到底对系统行为贡献了什么。
