宁波建设银行网站首页,微商商城系统,wordpress 引用网页,十大不收费看盘软件排名在建立常微分方程后#xff0c;我们常常会问的一个问题是#xff1a; 这个常微分方程系统中是否存在平衡点#xff1f;如果存在#xff0c;这个平衡点是否稳定#xff1f; 本节将展示如何通过ecode包来回答这个问题。
我们首先利用ecode包建立一个Lotka–Volterra竞争模型…在建立常微分方程后我们常常会问的一个问题是 这个常微分方程系统中是否存在平衡点如果存在这个平衡点是否稳定 本节将展示如何通过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). dtdx1f1(x1,x2,...,xn),dtdx2f2(x1,x2,...,xn),...dtdxmfn(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 xx_0 xx0随后不断进行如下运算 x n 1 x n − f ( x n ) f ′ ( x n ) , n ∈ N x_{n1}x_n-\frac{f(x_n)}{f(x_n)},\quad n∈\textbf N xn1xn−f′(xn)f(xn),n∈N 此时 x n 1 x_{n1} xn1相较于 x n x_n xn而言更接近于方程 f ( x ) 0 f(x)0 f(x)0的真实解因为 x n 1 x_{n1} xn1在 x n x_n xn的基础上沿着函数 f ( x ) f(x) f(x)的梯度下降。当 ∣ x n 1 − x n ∣ |x_{n1}-x_n| ∣xn1−xn∣小于一个精度值 ε ε ε时Newton迭代法停止 x n 1 x_{n1} xn1为方程 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_{n1} xn1可能陷入函数 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_{n1}\textbf x_n - \textbf J^{-1}(\textbf x_n)\textbf f(\textbf x_n) xn1xn−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 0x,y1000 0x,y1000因而Newton迭代法失败。
由此可见利用eode_get_cripoi函数探寻一个常微分方程的平衡点时需要仔细设置初始解并尝试多个不同的初始解从而得到比较满意的结果。
