在建立常微分方程后,我们常常会问的一个问题是: 这个常微分方程系统中是否存在平衡点?如果存在,这个平衡点是否稳定? 本节将展示如何通过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
x
x代表物种1的种群个体数,
x
x
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
n
n元方程组即可。
然而,在现实中, f 1 , f 2 , . . . , f n f_1,f_2,...,f_n f1,f2,...,fn可能是非常复杂的函数,导致很难得出该方程组的解析解。因而,我们常常采用数值方法对其进行求解。
ecode
包目前采用Newton迭代法法来求解常微分方程的平衡点。
Newton迭代法
Newton迭代法是一种用于求解方程组的数值方法。如果要求解一个一元方程,
f ( x ) = 0 f(x)=0 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 f(x)=0 f(x)=0的真实解,因为 x n + 1 x_{n+1} xn+1在 x n x_n xn的基础上沿着函数 f ( x ) f(x) f(x)的梯度下降。当 ∣ x n + 1 − x n ∣ |x_{n+1}-x_n| ∣xn+1−xn∣小于一个精度值 ε ε ε时,Newton迭代法停止, x n + 1 x_{n+1} xn+1为方程 f ( x ) = 0 f(x)=0 f(x)=0的近似解。
需要注意,Newton迭代法并不是总能求出方程 f ( x ) = 0 f(x)=0 f(x)=0的解。当 x 0 x_0 x0的值设定的不合适时,经过Newton迭代后, x n + 1 x_{n+1} xn+1可能陷入函数 f ( x ) f(x) f(x)的某个大于0的极小值中,或者落到 f ( x ) f(x) f(x)的定义域外。
当 f ( x ) \textbf f(\textbf x) f(x)是一个多元方程组时,需要求解的方程是
f ( x ) = 0 \textbf f(\textbf x)=0 f(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为多元函数 f ( x ) \textbf f(\textbf x) f(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
x
x,另一个维度是
y
y
y,该相点的值是
(
x
,
y
)
=
(
1
,
1
)
(x,y)=(1,1)
(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
)
(1-2)
(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
:精度值 ε ε ε。当两个迭代前后的两个相点之间的距离小于精度值 ε ε ε时,Newton迭代法停止,该函数返回计算结果。该参数的默认值为0.001。
代码运行后的结果表明,当用户指定 ( 0.5 , 0.5 ) (0.5,0.5) (0.5,0.5)为初始解时,Newton迭代法能够顺利地找出模型 ( 1 − 2 ) (1-2) (1−2)的平衡点 ( 1 / 3 , 1 / 3 ) (1/3,1/3) (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 ) (1,5) (1,5)为起始解,Newton迭代法就无法找到平衡点。屏幕中显示了Newton迭代法的迭代过程。当迭代到第6步时,相点已经在 ( − 0.0002 , 0.9996 ) (-0.0002,0.9996) (−0.0002,0.9996),这已经超出了模型变量的取值范围 0 < x , y < 1000 0<x,y<1000 0<x,y<1000,因而Newton迭代法失败。
由此可见,利用eode_get_cripoi
函数探寻一个常微分方程的平衡点时,需要仔细设置初始解,并尝试多个不同的初始解,从而得到比较满意的结果。