常微分方程建模R包ecode三——探寻平衡点

在建立常微分方程后，我们常常会问的一个问题是： 这个常微分方程系统中是否存在平衡点？如果存在，这个平衡点是否稳定？ 本节将展示如何通过ecode包来回答这个问题。

我们首先利用ecode包建立一个Lotka–Volterra竞争模型（详见第一节），
d x d t = ( r 1 − a 11 x − a 12 y ) x , ( 1 ) d y d t = ( r 2 − a 21 x − a 22 y ) y , ( 2 ) \frac{dx}{dt}=(r_1-a_{11}x-a_{12}y)x, \quad (1) \\ \frac{dy}{dt}=(r_2-a_{21}x-a_{22}y)y, \quad (2) dtdx=(r1−a11x−a12y)x,(1)dtdy=(r2−a21x−a22y)y,(2)
其中，$x$代表物种1的种群个体数，$x$代表物种2的种群个体数， r 1 , r 2 r_1, r_2 r1,r2为种群增长率， a 11 , a 12 , a 21 , a 22 a_{11},a_{12},a_{21},a_{22} a11,a12,a21,a22为两物种之间的竞争系数。

library(ecode)

eq1 <- function(x, y, r1 = 1, a11 = 1, a12 = 2) (r1 - a11 * x - a12 * y) * x
eq2 <- function(x, y, r2 = 1, a21 = 2, a22 = 1) (r2 - a21 * x - a22 * y) * y
x <- eode(dxdt = eq1, dydt = eq2)

一、数学原理

对任意一个常微分方程组，
d x 1 d t = f 1 ( x 1 , x 2 , . . . , x n ) , d x 2 d t = f 2 ( x 1 , x 2 , . . . , x n ) , . . . d x m d t = f n ( x 1 , x 2 , . . . , x n ) . \frac{dx_1}{dt}=f_1(x_1,x_2,...,x_n),\\ \frac{dx_2}{dt}=f_2(x_1,x_2,...,x_n),\\ ...\\ \frac{dx_m}{dt}=f_n(x_1,x_2,...,x_n). dtdx1=f1(x1,x2,...,xn),dtdx2=f2(x1,x2,...,xn),...dtdxm=fn(x1,x2,...,xn).
若要求出其平衡点 ( x 1 ∗ , x 2 ∗ , . . . , x n ∗ ) (x_1^*,x_2^*,...,x_n^*) (x1∗,x2∗,...,xn∗)，则需令，
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. 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. f1(x1∗,x2∗,...,xn∗)=0,f2(x1∗,x2∗,...,xn∗)=0,...fn(x1∗,x2∗,...,xn∗)=0.
然后求解这个$n$元方程组即可。

然而，在现实中， f 1 , f 2 , . . . , f n f_1,f_2,...,f_n f1,f2,...,fn可能是非常复杂的函数，导致很难得出该方程组的解析解。因而，我们常常采用数值方法对其进行求解。

ecode包目前采用Newton迭代法法来求解常微分方程的平衡点。

Newton迭代法
Newton迭代法是一种用于求解方程组的数值方法。如果要求解一个一元方程，
$$f(x)=0$$
则首先随机猜测这个方程组的解是 x = x 0 x=x_0 x=x0，随后不断进行如下运算，
x n 1 = x n − f ( x n ) f ′ ( x n ) , n ∈ N x_{n 1}=x_n-\frac{f(x_n)}{f'(x_n)},\quad n∈\textbf N xn 1=xn−f′(xn)f(xn),n∈N
此时 x n 1 x_{n 1} xn 1相较于 x n x_n xn而言更接近于方程$f(x)=0$的真实解，因为 x n 1 x_{n 1} xn 1在 x n x_n xn的基础上沿着函数$f(x)$的梯度下降。当 ∣ x n 1 − x n ∣ |x_{n 1}-x_n| ∣xn 1−xn∣小于一个精度值$\epsilon$时，Newton迭代法停止， x n 1 x_{n 1} xn 1为方程$f(x)=0$的近似解。
需要注意，Newton迭代法并不是总能求出方程$f(x)=0$的解。当 x 0 x_0 x0的值设定的不合适时，经过Newton迭代后， x n 1 x_{n 1} xn 1可能陷入函数$f(x)$的某个大于0的极小值中，或者落到$f(x)$的定义域外。
当$\text{f}(\text{x})$是一个多元方程组时，需要求解的方程是
$$\text{f}(\text{x})=0$$
或者表达为，
f 1 ( x 1 , x 2 , . . . , x k ) = 0 f 2 ( x 1 , x 2 , . . . , x k ) = 0 . . . f k ( x 1 , x 2 , . . . , x k ) = 0 f_1(x_1,x_2,...,x_k)=0\\ f_2(x_1,x_2,...,x_k)=0\\...\\ f_k(x_1,x_2,...,x_k)=0 f1(x1,x2,...,xk)=0f2(x1,x2,...,xk)=0...fk(x1,x2,...,xk)=0
假设设定的初始解为 x 0 \textbf x_0 x0，那么随后就要不断进行如何运算，
x n 1 = x n − J − 1 ( x n ) f ( x n ) \textbf x_{n 1}=\textbf x_n - \textbf J^{-1}(\textbf x_n)\textbf f(\textbf x_n) xn 1=xn−J−1(xn)f(xn)
其中， J − 1 \textbf J^{-1} J−1为多元函数$\text{f}(\text{x})$的Jacobian矩阵，
[ ∂ f 1 ∂ x 1 ∂ f 1 ∂ x 2 . . . ∂ f 1 ∂ x k ∂ f 2 ∂ x 1 ∂ f 2 ∂ x 2 . . . ∂ f 2 ∂ x k . . . ∂ f k ∂ x 1 ∂ f k ∂ x 2 . . . ∂ f k ∂ x k ] \begin{bmatrix} \frac{∂f_1}{∂x_1} & \frac{∂f_1}{∂x_2} & ... & \frac{∂f_1}{∂x_k}\\ \frac{∂f_2}{∂x_1} & \frac{∂f_2}{∂x_2} & ... & \frac{∂f_2}{∂x_k}\\ ...\\ \frac{∂f_k}{∂x_1} & \frac{∂f_k}{∂x_2} & ... & \frac{∂f_k}{∂x_k} \end{bmatrix} ∂x1∂f1∂x1∂f2...∂x1∂fk∂x2∂f1∂x2∂f2∂x2∂fk.........∂xk∂f1∂xk∂f2∂xk∂fk

二、pp对象

上一节中提到，**相点（phase point）**是相空间中的任意一点，用于描述一个常微分方程系统在某一时刻的状态。

在ecode包中，有专门的pp对象来描述一个相点。若要创建一个pp对象，就需要调用其构造函数pp()

point <- pp(list(x = 1, y = 1))
point
## phase point:
## x = 1 
## y = 1 

该函数创造了一个相点对象point。该相点位于二维空间中，一个维度为$x$，另一个维度是$y$，该相点的值是$(x,y)=(1,1)$。

ecode包中提供相点的四则运算函数，例如

pp(list(x = 1, y = 1))   pp(list(x = 2, y = 3))
## phase point:
## x = 3 
## y = 4

对两个位于同一空间内的相点施行加法运算，则运算后的结果也是同一空间内的相点，相点的坐标是原来两个相点坐标的和。同样，减法、乘法和除法的运算规则类似。

三、eode_get_cripoi函数

ecode包中提供eode_get_cripoi函数来获取一个常微分方程系统的平衡点。我们现在可以尝试获取常微分方程系统$(1-2)$的平衡点，

eq_point <- eode_get_cripoi(x, init_value = pp(list(x = 0.5, y = 0.5)), 
                            eps = 0.001)
eq_point
## phase point:
## x = 0.3333333 
## y = 0.3333333 

函数eode_get_cripoi中各个参数的含义是，

• x：常微分方程系统。一个eode对象。
• init_value：初始解 x 0 \textbf x_0 x0。一个pp对象， 其中pp对象的维度要与常微分方程系统中模型变量的名字对应。如前所述，用户必须指定一个初始解，这样Newton迭代法才能开始运行。
• eps：精度值$\epsilon$。当两个迭代前后的两个相点之间的距离小于精度值$\epsilon$时，Newton迭代法停止，该函数返回计算结果。该参数的默认值为0.001。

代码运行后的结果表明，当用户指定$(0.5,0.5)$为初始解时，Newton迭代法能够顺利地找出模型$(1-2)$的平衡点$(1/3,1/3)$。

然而，如果我们换一个初始解，

eq_point <- eode_get_cripoi(x, init_value = pp(list(x = 1, y = 5)), 
                            eps = 0.001)
## Fail to find an equilibrium point. The function will return iteration history
##               x         y
## 1             1         5
## 2     0.5049506  2.722774
## 3     0.2525784  1.610832
## 4     0.1155298  1.113498
## 5    0.03020005 0.9806953
## 6 -0.0002282437 0.9996615
## Error in eode_get_cripoi(x, init_value = pp(list(x = 1, y = 5)), eps = 0.001) : 
##   Fail to find an equilibrium point. Please try other initial values

这表明，如果以$(1,5)$为起始解，Newton迭代法就无法找到平衡点。屏幕中显示了Newton迭代法的迭代过程。当迭代到第6步时，相点已经在$(-0.0002,0.9996)$，这已经超出了模型变量的取值范围$0<x,y<1000$，因而Newton迭代法失败。

由此可见，利用eode_get_cripoi函数探寻一个常微分方程的平衡点时，需要仔细设置初始解，并尝试多个不同的初始解，从而得到比较满意的结果。

这篇好文章是转载于：学新通技术网