本节中我们考虑一个更为复杂的常微分方程模型，
d
X
C
d
t
=
ν
(
X
A
+
Y
A
)
−
β
⋅
X
C
⋅
(
Y
C
+
Y
A
)
−
(
μ
+
g
)
⋅
X
C
,
(
1
)
d
Y
C
d
t
=
β
⋅
X
C
⋅
(
Y
C
+
Y
A
)
−
(
μ
+
g
+
ρ
)
⋅
Y
C
,
(
2
)
d
X
A
d
t
=
g
⋅
X
C
−
β
⋅
X
A
⋅
(
Y
C
+
Y
A
)
,
(
3
)
d
Y
A
d
t
=
β
⋅
X
A
∗
(
Y
C
+
Y
A
)
+
g
∗
Y
C
−
ρ
∗
Y
A
(
4
)
\frac{dX_C}{dt} = \nu (X_A + Y_A) - \beta · X_C · (Y_C + Y_A) - (\mu + g) · X_C, \quad(1) \\ \frac{dY_C}{dt} = \beta · X_C · (Y_C + Y_A) - (\mu + g + \rho) · Y_C, \quad(2)\\ \frac{dX_A}{dt} = g · X_C - \beta · X_A · (Y_C + Y_A), \quad (3) \\ \frac{dY_A}{dt} = \beta · X_A * (Y_C + Y_A) + g * Y_C - \rho * Y_A \quad(4)
dtdXC=ν(XA+YA)−β⋅XC⋅(YC+YA)−(μ+g)⋅XC,(1)dtdYC=β⋅XC⋅(YC+YA)−(μ+g+ρ)⋅YC,(2)dtdXA=g⋅XC−β⋅XA⋅(YC+YA),(3)dtdYA=β⋅XA∗(YC+YA)+g∗YC−ρ∗YA(4)
该常微分方程系统用于模拟一种树木传染病动态，其中
X
C
X_C
XC代表易感树苗（susceptible sapling）的个体数，
Y
C
Y_C
YC代表感病树苗（infected sapling）的个体数，
X
A
X_A
XA代表易感成年树木的个体数（susceptible adult），
Y
A
Y_A
YA代表感病成年树木（infected adult）的个体数。显然
X
C
,
Y
C
,
X
A
,
Y
A
≥
0
X_C,Y_C,X_A,Y_A≥0
XC,YC,XA,YA≥0。
( 1 ) (1) (1)式中， ν ( X A + Y A ) \nu (X_A + Y_A) ν(XA+YA)代表繁殖产生新树苗的速率，其中 ν \nu ν是繁殖速率。 − β ⋅ X C ⋅ ( Y C + Y A ) -\beta · X_C · (Y_C + Y_A) −β⋅XC⋅(YC+YA)指的是由于疾病传染，导致易感树苗转化为感病树苗的速率，其中 β \beta β为传染率。 − ( μ + g ) ⋅ X C , - (\mu + g) · X_C, −(μ+g)⋅XC,是指由于自然死亡及树木成长，导致易感树苗被移除，或是进入到易感成年树木的速率，其中 μ \mu μ为自然死亡率， g g g为树木生长率。
( 2 ) (2) (2)式中， β ⋅ X C ⋅ ( Y C + Y A ) \beta · X_C · (Y_C + Y_A) β⋅XC⋅(YC+YA)代表由易感树苗被传染而转化为感病树苗的速率， − ( μ + g + ρ ) ⋅ Y C - (\mu + g + \rho) · Y_C −(μ+g+ρ)⋅YC则包含自然死亡、树木成长、因疾病感染而死亡这三个过程，其中 ρ \rho ρ代表由于传染病而导致的死亡率。
( 3 ) (3) (3)式中， g ⋅ X C g · X_C g⋅XC代表由于树木生长而使易感树苗转换为易感成年树木的速率， − β ⋅ X A ⋅ ( Y C + Y A ) - \beta · X_A · (Y_C + Y_A) −β⋅XA⋅(YC+YA)代表由于疾病传染，使易感成年大树转换为感病成年大树的速率。
( 4 ) (4) (4)式中， β ⋅ X A ∗ ( Y C + Y A ) \beta · X_A * (Y_C + Y_A) β⋅XA∗(YC+YA)对应于疾病传染过程， g ∗ Y C g * Y_C g∗YC对应于树木生长过程， − ρ ∗ Y A -\rho * Y_A −ρ∗YA对应于疾病导致的死亡过程。
研究一个常微分方程系统，最为直接的方法是研究其相速矢量场（phase velocity vector filed）。下面我们回顾一下与相速矢量场相关的几个重要概念，
常微分方程中的几个重要概念
相空间（phase space）：是指所有模型变量的所有可能取值的组合构成的空间。在本节案例中，相空间为 { ( X C , Y C , X A , Y A ) ∣ X C , Y C , X A , Y A ≥ 0 } \{(X_C,Y_C,X_A,Y_A)| X_C,Y_C,X_A,Y_A≥0\} {(XC,YC,XA,YA)∣XC,YC,XA,YA≥0}。
相点（phase point）：相空间中的任意一个点称为相点。相点用于描述系统的状态。在本节案例中，相点 ( X C , Y C , X A , Y A ) = ( 10 , 60 , 20 , 100 ) (X_C,Y_C,X_A,Y_A)=(10,60,20,100) (XC,YC,XA,YA)=(10,60,20,100)代表系统中有10棵易感树苗、60棵感病树苗、20棵易感成树、100棵感病成树。
相速矢量（phase velocity vector）：系统位于某一相点时，其速度大小和方向构成的矢量，叫做该相点所对应的相速矢量。在本节案例中，针对相点 ( X C , Y C , X A , Y A ) = ( 10 , 60 , 20 , 100 ) (X_C,Y_C,X_A,Y_A)=(10,60,20,100) (XC,YC,XA,YA)=(10,60,20,100)，将 X C , Y C , X A , Y A X_C,Y_C,X_A,Y_A XC,YC,XA,YA的值代入式 ( 1 − 4 ) (1-4) (1−4)，求出 ( d X C d t , d Y C d t , d X A d t , d Y A d t ) (\frac{dX_C}{dt}, \frac{dY_C}{dt} , \frac{dX_A}{dt} , \frac{dY_A}{dt}) (dtdXC,dtdYC,dtdXA,dtdYA)，其值便是该相点所对应的相速矢量。相速矢量描述了系统位于某一相点时的运动方向和快慢。
相速矢量场（phase velocity vector field）：相空间中所有相速矢量组成的集合。
相位曲线（phase curve）：相点沿相速矢量场移动所形成的运动轨迹。
在ecode包中，当函数plot作用于eode类的对象时，plot函数会自动绘制出某个常微分方程系统的相速矢量场，或相速矢量。在上一节中，我们介绍了当常微分方程系统中含有两个模型变量时，plot函数的用法。
本节所关注的模型含有四个模型变量
X
C
,
Y
C
,
X
A
,
Y
A
X_C,Y_C,X_A,Y_A
XC,YC,XA,YA，因而将介绍含有多个模型变量时，plot函数的行为。
一、绘制相速矢量场
首先我们在ecode包中构建上述模型（式
(
1
−
4
)
(1-4)
(1−4)）：
library(ecode)
dX_Cdt <- function(X_C, Y_C, X_A, Y_A, nu = 0.15, beta = 0.1, mu = 0.15, g = 0.04)
nu * (X_A + Y_A) - beta * X_C * (Y_C + Y_A) - (mu + g) * X_C
dY_Cdt <- function(X_C, Y_C, Y_A, beta = 0.1, mu = 0.15, g = 0.04, rho = 0.2)
beta * X_C * (Y_C + Y_A) - (mu + g + rho) * Y_C
dX_Adt <- function(X_C, Y_C, X_A, Y_A, beta = 0.1, g = 0.04)
g * X_C - beta * X_A * (Y_C + Y_A)
dY_Adt <- function(X_A, Y_C, Y_A, beta = 0.1, g = 0.04, rho = 0.2)
beta * X_A * (Y_C + Y_A) + g * Y_C - rho * Y_A
x <- eode(dX_Cdt = dX_Cdt, dY_Cdt = dY_Cdt, dX_Adt = dX_Adt, dY_Adt = dY_Adt)
x
### An ODE system: 4 equations
### equations:
### dX_Cdt = nu * (X_A + Y_A) - beta * X_C * (Y_C + Y_A) - (mu + g) * X_C
### dY_Cdt = beta * X_C * (Y_C + Y_A) - (mu + g + rho) * Y_C
### dX_Adt = g * X_C - beta * X_A * (Y_C + Y_A)
### dY_Adt = beta * X_A * (Y_C + Y_A) + g * Y_C - rho * Y_A
### variables: X_C Y_C X_A Y_A
### parameters: nu = 0.15, beta = 0.1, mu = 0.15, g = 0.04, rho = 0.2
### constraints: X_A<1000 X_A>0 X_C<1000 X_C>0 Y_A<1000 Y_A>0 Y_C<1000 Y_C>0
由于我们没有在模型中指定模型变量的范围，ecode包自动将变量范围设置在
(
0
,
1000
)
(0,1000)
(0,1000)内。此时调用plot函数，
plot(x)
### Set X_A = 500, Y_A = 500 for mapping in two axis
输出结果如图所示。plot函数自动限制了
X
A
=
500
,
Y
A
=
500
X_A = 500, Y_A = 500
XA=500,YA=500，并以
X
C
,
Y
C
X_C, Y_C
XC,YC为横、纵坐标系绘制相速矢量场。需要注意的是，该相速矢量图仅仅表示
X
A
,
Y
A
X_A,Y_A
XA,YA为固定值时的相速矢量场，该矢量场位于相空间内部的一个平面。
当常微分方程组含有多个模型变量时，如果不给
plot任何参数，则plot函数默认会以常微分方程中前两个变量为坐标轴绘制相速矢量场，而其余变量将会被赋上一个固定值，其值等于该模型变量范围的中值。
二、自定义模型变量的值
如果不希望以
X
A
=
500
,
Y
A
=
500
X_A = 500, Y_A = 500
XA=500,YA=500为限制条件，则可以在plot函数中添加set_covar参数，
plot(x, set_covar = list(X_A = 10, Y_A = 20))
此时，plot函数给出的是
X
A
=
10
,
Y
A
=
20
X_A=10, Y_A=20
XA=10,YA=20株时，以
X
C
,
Y
C
X_C,Y_C
XC,YC为坐标轴作出的相速矢量场。
如果想要固定 X C , Y C X_C, Y_C XC,YC，以 X A , Y A X_A, Y_A XA,YA为坐标轴作相速矢量场，则
plot(x, set_covar = list(X_C = 10, Y_C = 20))
此为
X
C
=
10
,
Y
C
=
20
X_C=10, Y_C=20
XC=10,YC=20时，以
X
A
,
Y
A
X_A, Y_A
XA,YA为坐标轴的相速矢量场。
三、自定义模型变量的范围
上一节中已经提到，可以采用eode_set_constraint重新设置模型变量的范围。例如，以下代码将
X
C
,
Y
C
,
X
A
,
Y
A
X_C, Y_C, X_A, Y_A
XC,YC,XA,YA的范围位置为
(
0
,
5
)
(0,5)
(0,5)，
x <- eode_set_constraint(x, new_constraint = c("X_C<5","Y_C<5","X_A<5","Y_A<5"))
x
### An ODE system: 4 equations
### equations:
### dX_Cdt = nu * (X_A + Y_A) - beta * X_C * (Y_C + Y_A) - (mu + g) * X_C
### dY_Cdt = beta * X_C * (Y_C + Y_A) - (mu + g + rho) * Y_C
### dX_Adt = g * X_C - beta * X_A * (Y_C + Y_A)
### dY_Adt = beta * X_A * (Y_C + Y_A) + g * Y_C - rho * Y_A
### variables: X_C Y_C X_A Y_A
### parameters: nu = 0.15, beta = 0.1, mu = 0.15, g = 0.04, rho = 0.2
### constraints: X_A<5 X_A>0 X_C<5 X_C>0 Y_A<5 Y_A>0 Y_C<5 Y_C>0
接下来，我们尝试固定 X A = 2 , Y A = 2 X_A = 2, Y_A = 2 XA=2,YA=2，以 X C , Y C X_C, Y_C XC,YC为坐标轴，绘制相速矢量场，
plot(x, set_covar = list(X_A = 2, Y_A = 2))
可以看到，该常微分方程组似乎在
X
C
,
Y
C
X_C, Y_C
XC,YC的值很小时存在使
d
X
C
/
d
t
,
d
Y
C
/
d
t
=
0
dX_C/dt, dY_C/dt=0
dXC/dt,dYC/dt=0的点。
四、一维相速矢量场
如果一个常微分方程只有一个模型变量，或者在含有
n
n
n个模型变量的常微分方程组中，有
(
n
−
1
)
(n-1)
(n−1)个变量的值都被固定了，那么plot函数就会绘制一维的相速矢量场。
现在我们尝试固定
X
A
=
2
,
Y
A
=
2
,
X
C
=
2
X_A = 2, Y_A = 2, X_C = 2
XA=2,YA=2,XC=2。这样，只剩下模型变量
X
A
X_A
XA的值没有被固定。plot函数将以
X
A
X_A
XA为唯一的坐标轴，绘制一维的相速矢量场，
plot(x, set_covar = list(X_A = 2, Y_A = 2, X_C = 2))
其中，每个相速矢量的值代表的是
d
Y
C
/
d
t
dY_C/dt
dYC/dt在某一相点的对应值。
五、相速矢量
当所有模型变量都被赋值时，plot函数将会作出某一相点所对应的相速矢量。在相空间中，相速矢量的起点在其对应的相点，长度代表相点在相点在该处运动的速率。例如，在本节案例中，位于相点
(
X
C
0
,
Y
C
0
,
X
A
0
,
Y
A
0
)
(X_{C0},Y_{C0},X_{A0},Y_{A0})
(XC0,YC0,XA0,YA0)的相速矢量的值为，
(
d
X
C
/
d
t
,
d
Y
A
/
d
t
,
d
X
A
/
d
t
,
d
Y
A
/
d
t
)
∣
X
C
=
X
C
0
,
Y
C
=
Y
C
0
,
X
A
=
X
A
0
,
Y
A
=
Y
A
0
(dX_C/dt,dY_A/dt,dX_A/dt,dY_A/dt)|_{X_{C}=X_{C0}, Y_{C}=Y_{C0}, X_{A}=X_{A0}, Y_{A}=Y_{A0}}
(dXC/dt,dYA/dt,dXA/dt,dYA/dt)∣XC=XC0,YC=YC0,XA=XA0,YA=YA0
调用plot时，plot函数会画出相速矢量在各个维度上的值，
plot(x, set_covar = list(X_A = 2, Y_A = 2, X_C = 2, Y_C = 2))
该图给出了相点
(
2
,
2
,
2
,
2
)
(2,2,2,2)
(2,2,2,2)所对应的相速矢量，其中横坐标的每一个标签代表相速矢量在某一维度上的分解值，即
d
X
C
/
d
t
∣
X
C
=
2
,
d
Y
C
/
d
t
∣
Y
C
=
2
,
d
X
A
/
d
t
∣
X
A
=
2
,
d
Y
A
/
d
t
∣
Y
A
=
2
dX_C/dt|_{X_C=2}, dY_C/dt|_{Y_C=2}, dX_A/dt|_{X_A=2}, dY_A/dt|_{Y_A=2}
dXC/dt∣XC=2,dYC/dt∣YC=2,dXA/dt∣XA=2,dYA/dt∣YA=2。
如何引用
Wu, H. (2023). ecode: An R package to investigate community dynamics in ordinary differential equation systems. bioRxiv, 2023-06.文章来源：https://www.toymoban.com/news/detail-621543.html
原文见bioRxiv。文章来源地址https://www.toymoban.com/news/detail-621543.html
到了这里，关于常微分方程建模R包ecode（二）——绘制相速矢量场的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
