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

我们首先利用ecode包建立一个Lotka–Volterra竞争模型（详见第一节），
$$\begin{gathered} \frac{dx}{dt}=(r_1-a_{11}x-a_{12}y)x,(1) \\ \frac{dy}{dt}=(r_2-a_{21}x-a_{22}y)y,(2) \end{gathered}$$
其中，$x$代表物种1的种群个体数，$x$代表物种2的种群个体数，$r_1,r_2$为种群增长率，$a_{11},a_{12},a_{21},a_{22}$为两物种之间的竞争系数。

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)

一、数学原理

对任意一个常微分方程组，
$$\begin{gathered} \frac{d{x}_1}{dt}=f_1(x_1,x_2,...,x_n), \\ \frac{d{x}_2}{dt}=f_2(x_1,x_2,...,x_n), \\ ... \\ \frac{d{x}_m}{dt}=f_n(x_1,x_2,...,x_n). \end{gathered}$$
若要求出其平衡点$(x_1^{\ast },x_2^{\ast },...,x_n^{\ast })$，则需令，
$$\begin{gathered} f_1(x_1^{\ast },x_2^{\ast },...,x_n^{\ast })=0, \\ {f}_2(x_1^{\ast },x_2^{\ast },...,x_n^{\ast })=0, \\ ... \\ {f}_n(x_1^{\ast },x_2^{\ast },...,x_n^{\ast })=0. \end{gathered}$$
然后求解这个$n$元方程组即可。

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

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

Newton迭代法
Newton迭代法是一种用于求解方程组的数值方法。如果要求解一个一元方程，
$$f(x)=0$$
则首先随机猜测这个方程组的解是$x=x_0$，随后不断进行如下运算，
$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)},n\in \text{N}$$
此时$x_{n+1}$相较于$x_n$而言更接近于方程$f(x)=0$的真实解，因为$x_{n+1}$在$x_n$的基础上沿着函数$f(x)$的梯度下降。当$|x_{n+1}-x_n|$小于一个精度值$\epsilon$时，Newton迭代法停止，$x_{n+1}$为方程$f(x)=0$的近似解。
需要注意，Newton迭代法并不是总能求出方程$f(x)=0$的解。当$x_0$的值设定的不合适时，经过Newton迭代后，$x_{n+1}$可能陷入函数$f(x)$的某个大于0的极小值中，或者落到$f(x)$的定义域外。
当$\text{f}(\text{x})$是一个多元方程组时，需要求解的方程是
$$\text{f}(\text{x})=0$$
或者表达为，
$$\begin{gathered} 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 \end{gathered}$$
假设设定的初始解为${\text{x}}_0$，那么随后就要不断进行如何运算，
$${\text{x}}_{n+1}={\text{x}}_n-{\text{J}}^{-1}({\text{x}}_n)\text{f}({\text{x}}_n)$$
其中，${\text{J}}^{-1}$为多元函数$\text{f}(\text{x})$的Jacobian矩阵，
$$\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& ...& \frac{\partial f_1}{\partial x_k}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& ...& \frac{\partial f_2}{\partial x_k}\\ ...& & & \\ \frac{\partial f_k}{\partial x_1}& \frac{\partial f_k}{\partial x_2}& ...& \frac{\partial f_k}{\partial x_k}\end{array}\right]$$

二、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：初始解${\text{x}}_0$。一个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函数探寻一个常微分方程的平衡点时，需要仔细设置初始解，并尝试多个不同的初始解，从而得到比较满意的结果。