群落组成的平衡点与极限环

研究一个动力学系统，例如群落动态，我们会问：群落组成是否会达到平衡？如果达到平衡，这个平衡是否稳健/稳定（robust/stability）？如果稳定，到达稳定点的时间（或问弛豫时间，relaxation time）又是多少？如果不能达到平衡点，群落动态是否又会收敛到一个极限环？如果会，这个极限环是否稳定？达到极限环的时间又是多少？如果即没有稳定点、又没有极限环，是否意味着群落最终会处于混沌状态？

本文接下来逐一解决这些问题。

一、平衡点

1.1 平衡点是否存在

在常微分方程理论中，平衡点又可称为驻点（critical point）、奇点（singular point）。假定群落动态可用自治系统表示：

$\begin{equation}\begin{cases}\dot{x_1}=f_1(x_1,x_2,...,x_n)\\\dot{x_2}=f_2(x_1,x_2,...,x_n)\\...\\\dot{x_n}=f_n(x_1,x_2,...x_n)\end{cases}\end{equation},\quad(1)$

该群落的平衡点满足：

$\begin{equation}\begin{cases}f_1(x_1^*,x_2^*,...,x_n^*)=0\\f_2(x_1^*,x_2^*,...,x_n^*)=0\\...\\f_n(x_1^*,x_2^*,...x_n^*)=0\end{cases}\end{equation}$

若 $x_1,...,x_n$ 均为某个物种的种群数量，则：

$f_i=r_i(x_1,x_2,...,x_n)x_i,\quad x_i>0$

其中， $r_i$ 为种群增长率（population growth rate，pgr）。平衡点可表示为：

$\begin{equation}\begin{cases}r_1(x_1^*,x_2^*,...,x_n^*)=0\\r_2(x_1^*,x_2^*,...,x_n^*)=0\\...\\r_n(x_1^*,x_2^*,...x_n^*)=0\end{cases}\end{equation}$

等式 $r_i(x_1^*,x_2^*,...,x_n^*)=0$ 称为物种i的零增长曲线（zero net growth isocline，ZNGI）。

1.2 平衡点是否稳定

平衡点分为稳定平衡点和不稳定平衡点。相点位于稳定平衡点时，尽管它受扰，一段时间后也会恢复到原来的位置。

将系统(1)线性化，

$\begin{equation}\begin{cases}\dot{x_1}=\frac{∂f_1}{∂x_1}x_1+\frac{∂f_1}{∂x_2}x_2+...+\frac{∂f_1}{∂x_n}x_n\\\dot{x_2}=\frac{∂f_2}{∂x_1}x_1+\frac{∂f_2}{∂x_2}x_2+...+\frac{∂f_2}{∂x_n}x_n\\...\\\dot{x_n}=\frac{∂f_n}{∂x_1}x_1+\frac{∂f_n}{∂x_n}x_2+...+\frac{∂f_n}{∂x_n}x_n\end{cases}\end{equation}$

其系统矩阵可表示为，

$\textbf A=\begin{pmatrix} \frac{∂f_1}{∂x_1}&\frac{∂f_1}{∂x_2}&...&\frac{∂f_1}{∂x_n}\\\frac{∂f_2}{∂x_1}&\frac{∂f_2}{∂x_2}&...&\frac{∂f_2}{∂x_n}\\...&...&...&...\\\frac{∂f_n}{∂x_1}&\frac{∂f_n}{∂x_2}&...&\frac{∂f_n}{∂x_n}\\\end{pmatrix}$

当 $\textbf A^*$ 的特征值 $\lambda_1^*,...,\lambda_n^*<0$ 时，平衡点 $(x_1^*,x_2^*,...,x_n^*)$ 是稳定平衡点，否则，它是不稳定平衡点。

若 $x_1,...,x_n$ 均为某个物种的种群数量，则：

$f_i=r_i(x_1,x_2,...,x_n)x_i,\quad x_i>0$

$\frac{∂f_i}{∂x_j}=\begin{equation}\begin{cases}r_i+x_i\frac{∂r_i}{∂x_j},\quad \text{when }j=i\\x_i\frac{∂r_i}{∂x_j},\quad\text{when }j\neq j\end{cases}\end{equation}$

此时系统矩阵为，

$\textbf A=\begin{pmatrix} r_1+x_1\frac{∂r_1}{∂x_1}&x_1\frac{∂r_1}{∂x_2}&...&x_1\frac{∂r_1}{∂x_n}\\x_2\frac{∂r_2}{∂x_1}&r_2+x_2\frac{∂r_2}{∂x_2}&...&x_2\frac{∂r_2}{∂x_n}\\...&...&...&...\\x_n\frac{∂r_n}{∂x_1}&x_n\frac{∂r_n}{∂x_2}&...&r_n+x_n\frac{∂r_n}{∂x_n}\\\end{pmatrix}$

因而，平衡点 $(x_1^*,x_2^*,...,x_n^*)$ 是稳定点的条件为：方程

$|\lambda\textbf E-\textbf A^*|=\left|\begin{array}{cccc} \lambda-(r_1^*+x_1^*(\frac{∂r_1}{∂x_1})^*)&-x_1^*(\frac{∂r_1}{∂x_2})^*&...&-x_1^*(\frac{∂r_1}{∂x_n})^*\\-x_2^*(\frac{∂r_2}{∂x_1})^*&\lambda-(r_2^*+x_2^*(\frac{∂r_2}{∂x_2})^*)&...&-x_2^*(\frac{∂r_2}{∂x_n})^*\\...&...&...&...\\-x_n^*(\frac{∂r_n}{∂x_1})^*&-x_n^*(\frac{∂r_n}{∂x_2})^*&...&\lambda-(r_n^*+x_n^*(\frac{∂r_n}{∂x_n})^*)\\\end{array}\right| =0$

的解必须都为负数。

对于双物种系统，平衡点 $(x_1^*,x_2^*,...,x_n^*)$ 的稳定条件可继续简化为：

$\begin{equation}\begin{cases} (\frac{∂f_1}{∂x_1})^*+(\frac{∂f_2}{∂x_2})^*<0 \\(\frac{∂f_1}{∂x_1})^*(\frac{∂f_2}{∂x_2})^*-(\frac{∂f_1}{∂x_2})^*(\frac{∂f_2}{∂x_1})^*>0 \end{cases}\end{equation}$

即，

$\begin{equation}\begin{cases} (\frac{∂r_1}{∂x_1})^*x_1^*+(\frac{∂r_2}{∂x_2})^*x_2^*<0 \\(\frac{∂r_1}{∂x_1})^*(\frac{∂r_2}{∂N_2})^*-(\frac{∂r_1}{∂x_2})^*(\frac{∂r_2}{∂x_1})^*>0 \end{cases}\end{equation}$

1.3 到达稳定点的时间

理论上，相点的轨迹无限接近于平衡点，而不会真正到达，但相点的移动速度却有快慢之分。就像一个小球要滚到山谷，沿着陡坡滚要比沿着缓坡滚要来得快。而弛豫时间就是用来度量山坡是否陡峭的一个指标，

$\frac{\varphi(t_{rel})-\varphi(∞)}{\varphi(0)-\varphi(∞)}=\frac1e$

其中， $\varphi(t)$ 是系统的积分曲线， $\varphi(t_{rel}),\varphi(0),\varphi(∞)$ 分别是时间等于弛豫时间、初始时刻、位于平衡点时刻时相点的位置。

1.4 案例：Lotka–Volterra模型

$\begin{equation}\begin{cases}\dot{x}=kx-axy\\\dot{y}=-ly+bxy \end{cases}\end{equation}$

其中 $x,y,k,a,l,b>0$ 。ZNGI为，

$\begin{equation}\begin{cases}kx^*-ax^*y^*=0\\-ly^*+bx^*y^*=0 \end{cases}\end{equation}$

平衡点为 $x^*=l/b,\ y^*=k/a$ 。系统矩阵为

$\textbf A=\begin{pmatrix}k-ay&-ax\\by&-l+bx\\\end{pmatrix}$

其在平衡点处的值为，

$\textbf A^*=\begin{pmatrix}0&-\frac{al}{b}\\\frac{kb}{a}&0\\\end{pmatrix}$

特征值 $\lambda_1=\sqrt{kl}\ i,\lambda_2=-\sqrt{kl}\ i$ ，实部为0，虚部共轭。

所以平衡点 $(l/b,k/a)$ 为不稳定的中心点。

1.5 案例：Rosenzweig–MacArthur模型

$\frac{dN}{dt}=f(N,P)=r_0(1-\frac NK)N-ρ_{max}\frac{N}{k+N}P,$

$\frac{dP}{dt}=g(N,P)=(\gammaρ_{max}\frac{N}{k+N}-u)P.$

其ZNGI为，

$(k+N^*)r_0(1-\frac {N^*}{k})=ρ_{max}P^*,$

$\gammaρ_{max}\frac{N^*}{k+N^*}=u.$

平衡点为，

$N^*=\frac{uk}{\gamma \rho_{max}-u},$

$P^*=\frac{r_0}{\rho_{max}}(1-\frac{N^*}{K})(k+N^*)=\frac{\gamma r_0k}{\gamma\rho_{max}-u}(1-\frac{uk}{(\gamma\rho_{max}-u)K}).$

因而，

$\frac{∂f}{∂N}=(1-\frac{2N}{K})r_0-\rho_{max}P\frac{k}{(k+N)^2},$

$\frac{∂f}{∂P}=-\rho_{max}\frac{N}{k+N},$

$\frac{∂g}{∂N}=\frac{P\gamma\rho_{max}k}{(k+N)^2},$

$\frac{∂g}{∂P}=\gamma\rho_{max}\frac{N}{k+N}-u.$

代入平衡点 $(N^*,P^*)$ 的值，可得，

$(\frac{∂f}{∂N})^*=\frac{-\gamma\rho_{max}k+\gamma\rho_{max}K-uk-uK}{\gamma\rho_{max}K}r_0u,$

$(\frac{∂g}{∂N})^*=\frac{\gamma\rho_{max}K-uK-uk}{K}·\frac{r_0}{\rho_{max}},$

$(\frac{∂f}{∂P})^*=-\frac{u}{\gamma},$

$(\frac{∂g}{∂P})^*=0.$

若 $(N^*,P^*)$ 为稳定点，则，

$\begin{equation}\begin{cases} (\frac{∂f}{∂N_1})^*+(\frac{∂g}{∂N_2})^*<0 \\(\frac{∂f}{∂N_1})^*(\frac{∂g}{∂N_2})^*-(\frac{∂f}{∂N_2})^*(\frac{∂g}{∂N_1})^*>0 \end{cases}\end{equation}$

即，

$\begin{equation}\begin{cases}\gamma\rho_{max}K-uK-uk>0\\\gamma\rho_{max}K-\gamma\rho_{max}k-uk-uK进一步可简化为，<img class=$

注意到 $N^*=\frac{K-k}{2}$ 刚好是二次函数 $P^*=\frac{r_0}{\rho_{max}}(1-\frac{N^*}{K})(k+N^*)$ 取最大值时自变量的值。

此为Rosenzweig–MacArthur模型具有稳定平衡点的条件。但该模型的相位曲线比较难求解，因而弛豫时间比较难计算。其实，当该模型不具有稳定平衡点时，模型会产生极限环，但极限环也同样很难求解析解。

二、极限环

研究极限环的难点在于它很难求出来，不像平衡点那样，直接通过联立几个ZNGI的方法求出。一旦求出极限环后，判断其稳定性反而相对较容易。

2.1 极限环是否存在

对于自治系统(1)，首先求出其平衡点。只有在平衡点的邻域才有可能出现极限环。假定某个平衡点 $\textbf x^*=(x_1^*,..,x_n^*)$ 附近存在极限环：

$g(\textbf x)=g(x_1,x_2,...,x_n)=0$

那么就意味着有：

$\text dg(\textbf x)/\text dt=0,\quad \text d^2g(\textbf x)/\text dt^2=0$

令 $V(\textbf x)=[g(\textbf x)-g(\textbf x^*)]^n$ ，可以证明，

$\text dV(\textbf x)/\text dt=0,\quad \text d^2V(\textbf x)/\text dt^2=0$

这样，找到极限环的关键就在于找到一个函数 $V(\textbf x)\geq 0$ ，使得其对t的一阶导数和二阶导数均为0。即：

「定理」 系统(1)的稳定点 $\textbf x^*=(x_1^*,..,x_n^*)$ 附近存在极限环 ⇔ 存在 $V(\textbf x)=V(x_1,...,x_n)\geq 0$ ，使得 $\text dV(\textbf x)/\text dt=0,\quad \text d^2V(\textbf x)/\text dt^2=0$ 。

若 $g(\textbf x)-g(\textbf x^*)\geq 0$ ，则n可以取1，此时的 $V(\textbf x)$ 可以形象地理解成是相点所在极限环的“半径”函数。而 $\text dV(\textbf x)/\text dt=0$ 可以理解为“极限环的半径不变”。

2.2 极限环是否稳定

假设自治系统(1)存在极限环 $g(\textbf x)=0$ ，该极限环包含平衡点 $\textbf x^*=(x_1^*,..,x_n^*)$ ，极限环对应“半径”函数为 $V(\textbf x)$ ，则：

$g(\textbf x)=0$ 外稳定 ⇔ 在极限环外侧 $\text dV(\textbf x)/\text dt<0,\quad \text d^2V(\textbf x)/\text dt^2>0$

$g(\textbf x)=0$ 外不稳定 ⇔ 在极限环外侧 $\text dV(\textbf x)/\text dt>0,\quad \text d^2V(\textbf x)/\text dt^2>0$

$g(\textbf x)=0$ 内稳定 ⇔ 在极限环内侧 $\text dV(\textbf x)/\text dt>0,\quad \text d^2V(\textbf x)/\text dt^2<img class=$ 内不稳定 ⇔ 在极限环内侧 $\text dV(\textbf x)/\text dt<0,\quad \text d^2V(\textbf x)/\text dt^2<0$

若 $g(\textbf x)=0$ 同时外稳定、内稳定，则它为稳定极限环（stable limit cycle）；若 $g(\textbf x)=0$ 内外均不稳定，则它为不稳定极限环（unstable limit cycle）；若 $g(\textbf x)=0$ 外稳定、内不稳定，或者内稳定、外不稳定，则它为半稳定极限环（semi-stable limit cycle）。

$g(\textbf x)=0$ 为中性稳定中心（neutrality-stable center） ⇔ 在极限环内外两侧都有 $\text dV(\textbf x)/\text dt=0,\quad \text d^2V(\textbf x)/\text dt^2=0$ 。

2.3 案例：教科书案例

$\dot{x_1}=x_2+k_1x_1(x_1^2+x_2^2-\beta^2)$

$\dot{x_2}=-x_1+k_2x_2(x_1^2+x_2^2-\beta^2)$

其中 $x_1,x_2∈\textbf R$ 。则 $(0,0)$ 是该系统的平衡点。对应于该点的极限环为

$g(x_1,x_2)=x_1^2+x_2^2-\beta^2=0$ （问题是，对于形式繁杂的系统而言，就很难猜出其极限环在哪里，因此这步很难）

因而， $g(x_1,x_2)-g(0,0)=x_1^2+x_2^2$ ，它刚好必然大于0，故取 $n=1$ ， $V(x_1,x_2)=x_1^2+x_2^2$ ，可以发现，

$\frac {dV(x_1,x_2)}{dt}=2(x_1^2+x_2^2-\beta^2)(k_1x_1^2+k_2x_2^2)$

$\frac{d^2V(x_1,x_2)}{dt^2}=4(x_1^2+x_2^2-\beta^2)[2k_1x_1^4+2k_2x_2^4+(k_1+k_2)^2x_1^2x_2^2-\beta^2(k_1^2x_1^2+k_2^2x_2^2)]$

可以判断，若取 $k_1=k_2=1$ ， $g(x_1,x_2)=0$ 不稳定；若取 $k_1=k_2=-1$ ， $g(x_1,x_2)=0$ 稳定。

2.4 案例：Lotka–Volterra模型

$\begin{equation}\begin{cases}\dot{x}=kx-axy\\\dot{y}=-ly+bxy \end{cases}\end{equation}$

其中， $x,y,k,a,l,b>0$ ，平衡点为 $x^*=l/b,\ y^*=k/a$ 。

用分离变量法求其积分曲线，得：

$g(x,y)=bx+ay-l\ln x-k\ln y+C=0$

则，

$g(x,y)-g(x^*,y^*)=bx+ay-l\ln x-k\ln y-(l+k-l\ln\frac lb-k\ln\frac ka)$

它恒为正数。因而取 $V(x,y)=g(x,y)-g(x^*,y^*)$ ，计算可得

$\frac {dV(x,y)}{dt}=0,\quad \frac{d^2V(x,y)}{dt^2}=0$

因此，即使在极限环内外两侧的邻域， $V(x,y)$ 也均为0，故 $g(x,y)=0$ 是中性稳定中心。

有关极限环如何求解这一问题依然所知甚少，欢迎交流。

参考文献

Ghaffari, A., Tomizuka, M. & Soltan, R.A. The stability of limit cycles in nonlinear systems. Nonlinear Dyn 56, 269–275 (2009). https://doi.org/10.1007/s11071-008-9398-3.

最后编辑于：2023.02.19 04:24:42

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 216,402评论 6赞 499
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,377评论 3赞 392
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 162,483评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,165评论 1赞 292
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,176评论 6赞 388
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,146评论 1赞 297
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,032评论 3赞 417
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,896评论 0赞 274
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,311评论 1赞 310
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,536评论 2赞 332
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,696评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,413评论 5赞 343
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,008评论 3赞 325
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,659评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,815评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,698评论 2赞 368
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,592评论 2赞 353