Lyapunov direct method for investigating stability of nonstandard finite difference schemes for metapopulation models

TL;DR

This paper develops Lyapunov-based methods to analyze the stability of nonstandard finite difference schemes for metapopulation models, ensuring the discrete models retain key properties of the continuous systems.

Contribution

It introduces a Lyapunov function approach to verify stability of NSFD schemes, simplifying the process compared to previous methods.

Findings

NSFD schemes preserve properties of metapopulation models

Lyapunov method simplifies stability analysis

Numerical examples confirm theoretical results

Abstract

In this paper nonstandard finite difference (NSFD) schemes of two metapopulation models are constructed. The stability properties of the discrete models are investigated by the use of a generalization of Lyapunov stability theorem. Due to this result we have proved that the NSFD schemes preserve all properties of the metapopulation models. Numerical examples confirm the obtained theoretical results of the properties of the constructed difference schemes. The method of Lyapunov functions proves to be much simpler than the standard method for studying stability of the discrete metapopulation model in our very recent paper.

Equations

$$\frac{d{p}_0}{dt}=e(p_1+p_2)-\lambda p_0,\frac{d{p}_1}{dt}=\lambda p_0-\beta p_1{p}_2+\delta p_2-e{p}_1,\frac{d{p}_2}{dt}=\beta p_1{p}_2-(\delta +e)p_2,$$

$$\frac{d{p}_1}{dt}=\lambda (1-p_1-p_2)-\beta p_1{p}_2+\delta p_2-e{p}_1,\frac{d{p}_2}{dt}=p_2(\beta p_1-\delta -e).$$

$$\big{(}p_{1}(0),p_{2}(0)\big{)}\in D_{2}:=\Big{\{}(x,y)\in\mathbb{R}^{2}:0\leq x,y;x+y\leq 1\Big{\}}.$$

$$P_{1}^{*}=\big{(}\dfrac{\lambda}{\lambda+e},0\big{)},\qquad P_{2}^{*}=(\dfrac{\delta+e}{\beta},1-\dfrac{e}{\lambda+e}-\dfrac{\delta+e}{\beta})=(\dfrac{\delta+e}{\beta},\dfrac{\lambda}{\lambda+e}-\dfrac{\delta+e}{\beta}).$$

$$\begin{array}{rl}& \frac{dI}{dt}={\beta }_{I}SI-e_{I}I+fL-gI,\frac{dS}{dt}=e_{I}I-{\beta }_{I}SI+fR-gS,\\ & \frac{dL}{dt}=gI-fL-e_{L}L+{\beta }_{L}RI,\frac{dR}{dt}=gS-fR+e_{L}L-{\beta }_{L}RI.\end{array}$$

$$\big{(}I(0),S(0),L(0),R(0)\big{)}\in D_{4}:=\big{\{}(I,S,L,R)\in\mathbb{R}^{4}:0\leq I,S,L,R;I+S+L+R=1\Big{\}}.$$

$$I_1^{\ast }=0,S_1^{\ast }=\frac{f}{f+g},L_1^{\ast }=0,R_1^{\ast }=\frac{g}{g+f}.$$

$$\begin{split}&I_{2}^{*}:=I^{*}=\dfrac{-b+\sqrt{b^{2}-4ac}}{2a},\quad S_{2}^{*}:=S^{*}=\dfrac{b+2\beta_{I}\beta_{L}\dfrac{f}{f+g}-\sqrt{b^{2}-4ac}}{2a},\\ &R_{2}^{*}:=R^{*}=\dfrac{g}{f+g}-\dfrac{\beta_{I}}{f}{I^{*}}^{2}+\big{(}\dfrac{\beta_{I}}{f+g}-\dfrac{g+e_{I}}{f}\big{)}I^{*},\\ &L_{2}^{*}:=L^{*}=1-I^{*}-S^{*}-R^{*}=\dfrac{g}{f+g}-R^{*}=\dfrac{\beta_{I}}{f}{I^{*}}^{2}-\big{(}\dfrac{\beta_{I}}{f+g}-\dfrac{g+e_{I}}{f}\big{)}I^{*},\end{split}$$

$$\begin{array}{rl}a=& {\beta }_I{\beta }_L,b={\beta }_I(f+e_L)+{\beta }_L(e_I+g)-{\beta }_I{\beta }_L\frac{f}{f+g},\\ c=& (f+e_L)(e_I-{\beta }_I\frac{f}{f+g})+g(e_L-{\beta }_L\frac{f}{f+g}).\end{array}$$

$$\begin{array}{rl}\frac{x_{k+1}-x_k}{\phi (h)}& =-c_1(\lambda +e)x_k-c_2(\lambda +e)x_{k+1}+c_3(\delta -\lambda )y_k+c_4(\delta -\lambda )y_{k+1}\\ & -c_5\beta x_k{y}_k-c_6\beta x_{k+1}{y}_k-c_7\beta x_k{y}_{k+1}-c_8\beta x_{k+1}{y}_{k+1}+\lambda ,\\ \frac{y_{k+1}-y_k}{\phi (h)}& =-c_1(\lambda +e)y_k-c_2(\lambda +e)y_{k+1}+c_3(\lambda -\delta )y_k+c_4(\lambda -\delta )y_{k+1}\\ & +c_5\beta x_k{y}_k+c_6\beta x_{k+1}{y}_k+c_7\beta x_k{y}_{k+1}+c_8\beta x_{k+1}{y}_{k+1},\end{array}$$

$$c_1+c_2=1,c_3=1,c_4=0,c_5+c_6=1,c_7=c_8=0,\phi (h)=h+\mathcal{O}(h^2).$$

$$\begin{array}{rl}& x_{k+1}=\frac{x_k-\phi c_1(\lambda +e)x_k+\phi (\delta -\lambda )y_k-\phi c_5\beta x_k{y}_k+\lambda \phi }{1+\phi c_2(\lambda +e)+\phi c_6\beta y_k},\\ & y_{k+1}=\frac{[1+\phi c_2(\lambda +e)+\phi c_6\beta y_k][y_k-\phi c_1(\lambda +e)y_k+\phi (\lambda -\delta )y_k+\phi c_5\beta x_k{y}_k]}{[1+\phi c_2(\lambda +e)+\phi c_6\beta y_k][1+c_2(\lambda +e)]}\\ & +\frac{\phi c_6\beta y_k[x_k-\phi c_1(\lambda +e)x_k+\phi (\delta -\lambda )y_k-\phi c_5\beta x_k{y}_k+\lambda \phi ]}{[1+\phi c_2(\lambda +e)+\phi c_6\beta y_k][1+c_2(\lambda +e)]}.\end{array}$$

$$c_1\le 0,c_2\ge 0.$$

$$c_5\le 0,c_6\ge 0,c_2\ge c_6,c_1\le -\frac{\delta }{\lambda +e},$$

$$V(x,y)=\big{(}x+y-\dfrac{\lambda}{\lambda+e}\big{)}^{2},\qquad(x,y)\in D_{2},$$

$$\begin{split}\Delta V(x_{k},y_{k})&=V(x_{k+1},y_{k+1})-V(x_{k},y_{k})\\ &=\big{(}x_{k+1}+y_{k+1}-x_{k}-y_{k}\big{)}\big{(}x_{k+1}+y_{k+1}+x_{k}+y_{k}-2\dfrac{\lambda}{\lambda+e}).\end{split}$$

$$G^{*}=\big{\{}(x,y)\in D_{2}:x+y=\dfrac{\lambda}{\lambda+e}\big{\}},$$

$$x_{k+1}=x_{k}+\dfrac{\varphi\beta x_{k}^{2}-\varphi\Big{[}\dfrac{\beta\lambda}{\lambda+e}+\big{(}\delta+e\big{)}\Big{]}x_{k}+\varphi\lambda\dfrac{\delta+e}{\lambda+e}}{1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\big{(}\dfrac{\lambda}{\lambda+e}-x_{k}\big{)}}.$$

$$P_{1,x}^{\ast }=\frac{\lambda }{\lambda +e},P_{2,x}^{\ast }=\frac{\delta +e}{\beta },$$

$$c_2>c^{\ast }:=\frac{\frac{\beta \lambda }{\lambda +e}+\delta +e+c_6\beta (2-\frac{\lambda }{\lambda +e})}{\lambda +e},$$

$$f(x,h):=x+\dfrac{\varphi\beta x^{2}-\varphi\Big{[}\dfrac{\beta\lambda}{\lambda+e}+\big{(}\delta+e\big{)}\Big{]}x+\varphi\lambda\dfrac{\delta+e}{\lambda+e}}{1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\big{(}\dfrac{\lambda}{\lambda+e}-x\big{)}}.$$

$$\begin{split}&\dfrac{\partial f(x,h)}{\partial x}=\dfrac{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}\big{]}\varphi\big{[}c_{2}(\lambda+e)+\beta c_{6}\dfrac{\lambda}{\lambda+e}-2\beta c_{6}x-(\dfrac{\beta\lambda}{\lambda+e}+\delta+e)\big{]}}{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}-\varphi c_{6}\beta x\big{]}^{2}}\\ &+\dfrac{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}\big{]}(1+2\varphi\beta x)}{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}-\varphi c_{6}\beta x\big{]}^{2}}+\dfrac{\varphi^{2}\beta^{2}(c_{6}^{2}-c_{6})x^{2}+\varphi^{2}c_{6}\beta\lambda\dfrac{\delta+e}{\lambda+e}}{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}-\varphi c_{6}\beta x\big{]}^{2}}.\end{split}$$

$$\dfrac{\partial f(x,h)}{\partial x}>\dfrac{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}\big{]}\varphi(\lambda+e)(c_{2}-c^{*})}{\big{[}1+\varphi c_{2}(\lambda+e)+\varphi\beta c_{6}\dfrac{\lambda}{\lambda+e}-\varphi c_{6}\beta x\big{]}^{2}}>0.$$

$$u_{k+1}=u_{k}-\dfrac{u_{k}^{2}-\varphi\big{[}\dfrac{\beta\lambda}{\lambda+e}-(\delta+e)\big{]}u_{k}}{1+\varphi c_{2}(\lambda+e)+\varphi c_{6}\beta u_{k}}.$$

$$f(x)>x\Longleftrightarrow g(x):=\varphi\beta x^{2}-\varphi\Big{[}\dfrac{\beta\lambda}{\lambda+e}+\big{(}\delta+e\big{)}\Big{]}x+\varphi\lambda\dfrac{\delta+e}{\lambda+e}>0.$$

$$\{\begin{array}{l}f(x)>x⟺0<x<P_{2,x}^{\ast },\\ \\ f(x)<x⟺x>P_{2,x}^{\ast }.\end{array}$$

$$c_{5}\leq 0,\quad c_{6}\geq 0,\quad c_{2}\geq\max\big{\{}c_{6},c^{*}\big{\}},\quad c_{1}\leq-\dfrac{\delta}{\lambda+e}.$$

$$\begin{array}{rl}\frac{S_{k+1}-S_k}{\phi }& =e_I{I}_k-{\beta }_I{S}_{k+1}{I}_k+f{R}_k-g{S}_k,\\ \frac{I_{k+1}-I_k}{\phi }& ={\beta }_I{S}_{k+1}{I}_k-e_I{I}_k+f{L}_K-g{I}_k,\\ \frac{R_{k+1}-R_k}{\phi }& =g{S}_k-f{R}_k+e_L{L}_k-{\beta }_L{R}_{k+1}{I}_k,\\ \frac{L_{k+1}-L_k}{\phi }& =g{I}_k-f{L}_k-e_L{L}_k+{\beta }_L{R}_{k+1}{I}_k.\end{array}$$

$$\phi (h)<\frac{1}{f+g},\forall h>0.$$

$$a_{k+1}=(1-\phi f-\phi g)a_k+\phi f,b_{k+1}=(1-\phi f-\phi g)b_k+\phi g.$$

Full text

{classcode}

37M05, 39A10, 65L12, 65L20

1 Introduction

Many phenomena and processes in physics, chemistry, biology, ecology, finance, environment etc. are modelled by ordinary or partial differential equations [1, 5, 6, 7, 17, 21, 25, 26, 35, 43, 44]. The solutions of these equations often possess specific properties, such as positivity, monotonicity, periodicity, stability and some invariant properties. In general, these differential equations are very complicated and it is difficult, even impossible, to find their exact solutions. Therefore, the study of numerical methods and simulation for the solution of these differential equations is very important. Among the numerical methods for differential equations the finite difference method is most popular. The theory of this method for ordinary and partial differential equations is developed rather completed in [24, 34, 42, 45]. We call the difference schemes in these books and all related sources standard finite difference schemes. In many nonlinear problems the standard difference schemes exhibit a serious drawback which is called ”numerical instabilities” [28, 29, 31]. Under this concept we have in mind the phenomena when the discrete models, for example, the difference schemes, do not preserve properties of the corresponding differential equations. In [28, 29, 30, 31] Mickens showed many examples and analysed the numerical instabilities when using standard difference schemes. In general, standard difference schemes preserve the properties of the differential equations only in the case if the discretization parameter $h$ is sufficiently small. Therefore, when studying dynamical models in large time intervals the selection of small time steps will requires very large computational effort, so these discrete models are inefficient. Besides, for some special dynamical problems standard difference schemes cannot preserve the properties of the problems for any step sizes.

In order to overcome the numerical instabilities phenomena in 1989 Mickens [27] introduced the concept Nonstandard Finite Difference (NSDF) schemes and after that has developed NSDF methods in many works, such as [28, 29, 30, 31]. According to Mickens, NSDF schemes are those constructed following a set of five basic rules. The NSDF schemes preserve main properties of the differential counterparts, such as positivity, monotonicity, periodicity, stability and some other invariants including energy and geometrical shapes. It should be emphasized that NSFD schemes can preserve all properties of the continuous models for any discretization parameters. The discrete models with these properties are called dynamically consistent [8, 4, 12, 16, 32, 40, 41].

For the last two decades NSDF methods have attracted attention from many researchers and achieved significant results [3, 11, 13, 14, 22, 23, 33, 36, 37, 39, 40, 46]. The property of stability of the set of equilibria of differential equations is one of these results because it plays the essential role in the study of asymptotical behaviour of the solutions of differential equations. The construction of difference schemes, which preserve the stability of the equilibrium points, is important in numerical simulation of differential equations. The difference schemes with this stability property is called elementary stable schemes [3, 9, 10, 46]. There are many works concerning the elementary stable schemes. The typical results are for general dynamical systems [9, 10] and for other specific systems [13, 38, 46] etc. One popular approach to the elementary stability is the investigation of Jacobian matrice of the discrete models at the equilibria, namely, determination of conditions ensuring that all eigenvalues of Jacobian matrice have moduli less or equal to 1. This is the necessary and sufficient condition for hyperbolic equilibrium points to be locally stable [1, 21]. The mentioned above approach has the following weaknesses and limitations:

It is applicable when all the equilibrium points are hyperbolic. To our best knowledge, at present no results on NSFD schemes preserving the stability of non-hyperbolic equilibrium points are available. 2. 2.

Even when all the equilibrium points are hyperbolic, the determination of the conditions ensuring all the eigenvalues of Jacobian to be in the unit ball, is very difficult. Theoretically, it is possible to use the Jury’s conditions [1, 15, 21] for finding these conditions, but in practice, this is extremely complex in many cases, for example, when the system of equations has large dimension or contains several parameters. 3. 3.

The consideration of Jacobian only guarantees the local stability meanwhile many models have the global stability.

In order to overcome the above shortcomings of the approach with the use of Jacobian it is useful to use the Lyapunov stability theory for proving the stability of both hyperbolic and non-hyperbolic equilibrium points [1, 15]. The drawback of this approach is that not always it is possible to find appropriate Lyapunov’s function. Nevertheless, in many specific problems it is easy to find associated Lyapunov’s functions. Then the stability of the equilibrium points may be established without the study of Jacobian matrices of discrete models. This is a perspective approach for many discrete models.

For illustrating this approach in this paper we consider two metapopulation models, one was proposed by Keymer [20] and another by Amarasekare [2]. This is the models with complex properties. By using a generalization of Lyapunov’s stability theorem [18] we construct NSDF schemes preserving stability properties of the models. This way is much simpler than the using Jacobian because it doesn’t require complicated calculations and difficult techniques.

The paper is organized as follows. In Section 2 we recall the two models of metapopulations. In Sections 3 and 4 we construct NSDF schemes dynamically consistent with the continuous models. Next, in Section 5 we report some numerical experiments for validating the obtained theoretical results. Finally, some concluding remarks are given in Section 6.

2 Mathematical models of metapopulations

2.1 Keymer’s metapopulation model

Consider the metapopulation model proposed by Keymer in 2000 [20]. It is described by the system of three nonlinear differential equations

[TABLE]

where $p_{0}$, $p_{1}$ and $p_{2}$ denote the proportion of uninhabitable patches, the proportion of the habitable patches that are not occupied and the proportion of habitable patches that are occupied, respectively, $\lambda$ is the rate of patch creation, $e$ is the rate of patch destruction, $\delta$ is the rate of population extinction and $\beta$ is the rate of propagule reproduction. Because of $p_{0}+p_{1}+p_{2}=1$, the system (1) can be reduced to the two equations

[TABLE]

From the biological meaning of the model we shall consider the initial conditions $p_{1}(0),p_{2}(0)$ satisfying

[TABLE]

The mathematical analysis in [1, 20] shows that the model (2) possesses the following properties:

$(P_{1})$

The monotone convergence of the sum $s(t):=p_{1}(t)+p_{2}(t)$:

With the initial conditions satisfying (3) the sum of the solutions $s(t):=p_{1}(t)+p_{2}(t)$ monotonically converges to $s^{*}:=\lambda/(\lambda+e)$.

$(P_{2})$

Boundedness

All the solution $p_{1}(t),p_{2}(t)$ with the initial conditions satisfying (3) also satisfy (3). In other words, the set $D_{2}$ is positive invariant.

$(P_{3})$

Local asymptotic stability

The model (2) has two equilibria

[TABLE]

Set $\mathcal{R}_{0}=\dfrac{\beta\lambda}{(\lambda+e)(\delta+e)}$. This is a threshold for persistence indicating the number of propagules needed during the species and the patch lifetime for the species to persist for the model.

According to [1, Chapter 6] the first equilibrium $P_{1}^{*}$ is locally asymptotically stable if $\mathcal{R}_{0}<1$ and the second equilibrium $P_{2}^{*}$ is locally asymptotically stable if $\mathcal{R}_{0}>1$.

$(P_{4})$

Global stability

Due to the Poincare - Bendixson theory [1, Chapter 5-6] it is possible to show that if $\mathcal{R}_{0}>1$ then the second equilibrium is globally asymptotically stable, and if $\mathcal{R}_{0}<1$ then the first equilibrium is globally asymptotically stable.

$(P_{5})$

Non-periodic solution

Applying Dulax’s criterion it is possible to show that the system (2) does not have periodic solutions in the domain $D=\{0<p_{1}+p_{2}<1\}$.

Clearly, this is a model with complex properties. Very recently, in [8] by using NSFD methods we have successfully constructed a discrete metapopulation model dynamically consistent with the continuous counterpart. It means that the discrete model preserves all the properties $(P_{1})-(P_{5})$ of the model (1). The construction of the NSFD scheme is very complicated. The main difficulty is in the stability properties $(P_{3})$ and $(P_{4})$. For constructing NSFD scheme preserving these properties it is needed to use difficult techniques and cumbersome calculations.

In Section 3 by using the Lyapunov’s stability theory we construct NSFD scheme preserving the global stability property of the model (1), so due to it the difference scheme preserves all the properties of the model. As will see later, the method of Lyapunov’s function is much simpler and shorter than the method for proving stabilities in (1).

2.2 Amarasekare-Possingham’s metapopulation model

Consider the metapopulation model proposed by Amarasekare and Possingham in 2001 [2]. It is described by the system of four nonlinear differential equations

[TABLE]

Here $f$ is the disturbance frequency and $g$, the rate of habitat succession. Quantities $e_{I}$ and $e_{L}$ represent local extinction rates, and $\beta_{I}$ and $\beta_{L}$ the per patch colonization rates of infected and latent patches, respectively. The total number of patches in the system is assumed to be constant such that $I+S+L+R=P$. Alternatively, $I,S,L$ and $R$ can be thought of as the frequency of each patch type in the landscape in which case $I+S+L+R=1$.

From the biological meaning of the model we shall consider the initial conditions $I(0),S(0)$, $L(0),R(0)$ satisfying

[TABLE]

The mathematical analysis shows that the model (4) possesses the following properties:

$(P_{1})$

The monotone convergence of the sum $a(t):=I(t)+S(t)$ and $b(t):=L(t)+R(t)$:

For any initial conditions satisfying (5) the sum $a(t):=I(t)+S(t)$ monotonically converges to $a^{*}:=f/(f+g)$, and $b(t):=L(t)+R(t)$ monotonically converges to $b^{*}:=g/(f+g)$.

$(P_{2})$

Boundedness:

All the solutions $I(t),S(t),L(t),R(t)$ with the initial conditions satisfying (5) also satisfy (5). In other words, the set $D_{4}$ is positive invariant.

$(P_{3})$

Local asymptotic stability (see [2]):

The model (5) has two equilibria (lying on the boundary or inside of $D_{4}$) $E_{i}^{*}=(I_{i}^{*},S_{i}^{*},L_{i}^{*},R_{i}^{*}),i=1,2$, where

[TABLE]

[TABLE]

[TABLE]

Set $\mathcal{R}_{0}=1-c$. Then $\mathcal{R}_{0}$ is threshold parameter, i.e., if $\mathcal{R}_{0}>1$ then the model has positive equilibrium point, otherwise positive equilibrium point does not exist. The first equilibrium point is locally asymptotically if and only if $c>0\,\,(\mathcal{R}_{0}<1)$ and the second equilibrium point is locally asymptotically if and only if $c<0\,\,(\mathcal{R}_{0}>1)$ .

3 NSFD scheme for the model (2)

In this section we construct NSFD scheme for the model (2) so that the obtained difference scheme preserves all dynamic properties of the original continuous model for any discretization parameter or step size $h>0$. According to Mickens, a finite difference scheme is called nonstandard if at least one of the following conditions is satisfied [28, 29, 30, 31]:

• •

A nonlocal approximation is used.

• •

The discretization of the derivative is not traditional and uses a function $0<\varphi(h)=h+\mathcal{O}(h^{2})$.

For simplicity of presentation, in the model (2) we use the notations $x(t)$ and $y(t)$ instead of $p_{1}(t)$ and $p_{2}(t)$, respectively. We reconsider the family of difference schemes of the form [8]

[TABLE]

where the parameters $c_{i}$ and the function $\varphi(h)$ satisfy

[TABLE]

The explicit form of the difference scheme is defined by [8, Theorem 2.4]

[TABLE]

3.1 Properties $(P_{1})$ and $(P_{2})$

The results of the difference schemes preserving the properties $(P_{1})$ and $(P_{2})$ of the model (2) are stated in the following theorems [8]:

Theorem 3.1**.**

The difference scheme (9)-(10) preserves Property $(P_{1})$ of the model (2) if

[TABLE]

Theorem 3.2**.**

Consider the scheme (9)-(10). Under the assumptions

[TABLE]

the scheme (9)-(10) preserves Property $(P_{2})$ of the model (2).

3.2 Stability properties $(P_{3})$ and $(P_{4})$

The main difficulty in the construction of difference schemes preserving the properties of the model (2) is in Properties $(P_{3})$ and $(P_{4})$, especially in the property of global stability $(P_{4})$. In this section by using a generalization of Lyapunov’s stability theorem [18, Theorem 3.3] we show the stability properties of the equilibrium points in a simple and easy way.

For this purpose we consider the function

[TABLE]

where $D_{2}$ is defined by (3). We shall show that the function $V(x,y)$ satisfies all the conditions of [18, Theorem 3.3] on $D_{2}$. Obviously, $V(x,y)$ is continuous on $D_{2}$. Moreover,

(i) $V(x,y)\geq 0$ for any $(x,y)\in D_{2}$ and $V(P_{i}^{*})=0$ $(i=1,2)$.

(ii)

[TABLE]

Since Property $(P_{1})$ of the model (2) is preserved it follows that $\Delta V(x_{k},y_{k})\leq 0$ for any $(x_{k},y_{k})\in D_{2}$. Besides $\Delta V(x_{k},y_{k})=0$ if and only if $x_{k}+y_{k}=\dfrac{\lambda}{\lambda+e}$. Therefore,

[TABLE]

where $G^{*}$ is the largest positively invariant set containing in $G=\{(x,y):\Delta V(x,y)=0\}$.

(iii) On the other hand, due to the preservation of Property $(P_{2})$ of the model all the solutions of (9) are bounded. Hence, all the conditions (1), (2) and (4) of [18, Theorem 3.3] are satisfied. In order to prove the global stability of (2) it remains to show that

If $\mathcal{R}_{0}<1$ then $P_{1}^{*}$ is $G^{*}-$ globally asymptotically stable, i.e., $P_{1}^{*}$ is $G^{*}$- asymptotically stable and $G^{*}$- globally attractive. 2. 2.

If $\mathcal{R}_{0}>1$ then $P_{2}^{*}$ is $G^{*}-$ globally asymptotically stable, i.e., $P_{2}^{*}$ is $G^{*}$- asymptotically stable and $G^{*}$- globally attractive.

Since Property $(P_{1})$ of the model (2) is preserved for any initial conditions belonging to $G^{*}$, i.e., $x_{0}+y_{0}=\dfrac{\lambda}{\lambda+e}$, the solutions of (9)-(10) satisfy $x_{k}+y_{k}=\lambda/(\lambda+e)$. Substituting $y_{k}=\lambda/(\lambda+e)-x_{k}$ into (11) we obtain the scheme depending only on $x_{k}$ of the form

[TABLE]

This is the discretization of the model (2) on the set \Big{\{}\big{(}x(t),y(t)\big{)}\in D_{2}:x(t)+y(t)=\lambda/(\lambda+e)\Big{\}}. Denote by

[TABLE]

$x-$components of the equilibrium points $P_{i}^{*}\ (i=1,2)$. This is also two equilibrium points of (15). Therefore, the proof of $G*$- global stability of the equilibrium points $P_{i}^{*}\ (i=1,2)$ is equivalent to showing that

If $\mathcal{R}_{0}<1$ then $P_{1,x}^{*}$ is the globally asymptotically stable equilibrium point of (15) in D_{2}^{*}:=\Big{\{}x\in\mathbb{R}:0\leq x\leq\dfrac{\lambda}{\lambda+e}\Big{\}}. 2. 2.

If $\mathcal{R}_{0}>1$ then $P_{2,x}^{*}$ is the globally asymptotically stable equilibrium point of (15) in $D_{2}^{*}$.

Of course, it is possible to show the local stability of $P_{1,x}^{*}$ and $P_{2,x}^{*}$ via the Jacobian matrice of the scheme (15), and doing this is not difficult. However, the scheme (15) depends only on $x_{k}$. Therefore, it is wise to use directly the results of Lubuma Anguelov [3, Theorem 3] to show the elementary stability of (15). Since the set of equilibrium points is preserved it suffices to show the monotone dependence on initial value of the scheme (15).

Theorem 3.3**.**

Consider the scheme (9)-(10) under the conditions of Theorem 3.2. If additionally assume that

[TABLE]

then the scheme (15) is elementary stable, i.e., $P_{1,x}^{*}$ and $P_{2,x}^{*}$ are locally asymptotically stable. Hence, the local stability Property $(P_{3})$ of the model (2) is preserved.

Proof.

Consider the function

[TABLE]

It is easy to obtain

[TABLE]

Since $c_{6}>1$ there holds $c_{6}^{2}>c_{6}$. On the other hand, in view of $0\leq x\leq 1$ from (19) and (17) it follows

[TABLE]

Therefore, (15) is monotonically depends on the initial value. The proof of the theorem is complete by using Theorem [3, Theorem 3] . ∎

The global attractiveness is obtained with the use of following theorem.

Theorem 3.4**.**

Consider the scheme (9)-(10). Under the assumptions of Theorems 3.2 and 3.3 then

If $\mathcal{R}_{0}<1$ then $P_{1,x}^{*}$ is globally stable equilibrium point of (15) in $D_{2}^{*}$. 2. 2.

If $\mathcal{R}_{0}>1$ then $P_{2,x}^{*}$ is globally stable equilibrium point of (15) in $D_{2}^{*}$.

Proof.

(i) Case $\mathcal{R}_{0}<1$.

If $\mathcal{R}_{0}<1$ then $P_{2,x}^{*}=\dfrac{\delta+e}{\beta}>\dfrac{\lambda}{\lambda+e}$. Therefore $P_{2,x}^{*}=\dfrac{\delta+e}{\beta}$ cannot be a globally stable point of (15) in $D_{2}^{*}$. We shall show that $P_{1,x}^{*}$ is a globally stable point of (15).

Set $u_{k}={\lambda}/{(\lambda+e)}-x_{k}$. Noticing that on $D_{2}^{*}$ there holds $u_{k}\geq 0$ for any $k$, from (15) we obtain

[TABLE]

Since $\mathcal{R}_{0}<1$ there holds $\dfrac{\beta\lambda}{\lambda+e}-(\delta+e)<0$. It follows $u_{k+1}<u_{k}$ for any $k$. The sequence \big{\{}u_{k}\big{\}} is decreasing and bounded from below, consequently, it is convergent. Therefore, the sequence \big{\{}x_{k}\big{\}} also is convergent. Since \big{\{}x_{k}\big{\}} cannot converge to $P_{2,x}^{*}=\dfrac{\delta+e}{\beta}$ it implies that $P_{1,x}^{*}=\dfrac{\lambda}{\lambda+e}$ is globally attractive point of (15) in $D_{2}^{*}$.

(ii) Case $\mathcal{R}_{0}>1$.

For simplicity, here instead of $f(x,h)$ we write $f(x)$. First, from (18) we have

[TABLE]

The function $g(x)$ has two positive roots $P_{2,x}^{*}<P_{1,x}^{*}$ (because $\mathcal{R}_{0}>1$). The equation $g^{\prime}(x)=0$ has a unique positive root $x^{*}=\dfrac{P_{1,x}^{*}+P_{2,x}^{*}}{2}$. From the investigating the behaviour of $g(x)$ we obtain

[TABLE]

Besides, by Theorem 3.3 the function $f(x)$ is increasing. Therefore, from the result [19, Problem 2.5.38] it follows that $P_{2,x}$ is a globally attractive point of (15). The proof is complete. ∎

Remark 1*.*

The function $V(x,y)$ constructed in the proof of the stability properties of the set of equilibrium points does not satisfy the conditions of the classical Lyapunov stability theorem. It satisfies a generalization of this theorem, namely [18, Theorem 3.3]. The properties $(P_{1}),(P_{2})$ play an important role in the construction of the function $V(x,y)$ for ensuring the stability of the NSFD scheme.

3.3 Non-periodicity of solution

Non-periodicity of solution of $(P_{5})$ follows from the proved fact that the equilibrium points are globally asymptotically stable.

Summarizing the above results we obtain the difference schemes preserving the properties of the model (2).

Theorem 3.5**.**

The difference schemes (9)-(10) preserve Properties $(P_{1})-(P_{5})$ of the model (2) if

[TABLE]

Remark 2*.*

Theorems 2.22, 2.23, 2.24 in [8] and Theorem 5 above give sufficient conditions for the NSFD schemes to preserve the properties of the metapopulation model (2). Similar to the scheme in Theorem 2.22 [8], the scheme in Theorem 5 contains 4 parameters $c_{1},c_{2},c_{5},c_{6}$. The conditions of the both theorems are satisfied for sufficiently small negative numbers $c_{1},c_{5}$ and sufficiently large positive numbers $c_{2},c_{6}$, and in this case the constructed difference schemes coincide. Nevertheless, Theorem 2.22 [8] only guarantees that the scheme preserves Properties $(P_{1})-(P_{3})$ of the model although numerical experiments show that Properties $(P_{4})$ and $(P_{5})$ also are preserved. Moreover, the derivation of the conditions in Theorem 2.22 is very hard because it is done via the analysis of the Jacobian for proving the local stability property $(P_{3})$. Using the approach in [8] it is difficult to find the conditions for the scheme to preserve Property $(P_{4})$ of the model.

The NSFD schemes satisfying Theorems 2.23, 2.24 in [8] and Theorem 5 preserve Properties $(P_{1})-(P_{5})$ of the model. Due to the complexity of the model the scheme in Theorem 2.23 contains only two parameters $c_{1}$ and $c_{2}$, and the scheme in Theorem 2.24 contains only one parameter being the denominator function $\varphi$. The scheme satisfying Theorem 5 does not satisfy Theorem 2.23 and Theorem 2.24, and conversely. However, the scheme in Theorem 5 still contains 4 parameters, and this allows to construct schemes with other properties, for example, the second order accuracy property. Besides, it should be emphasized that the derivation of conditions in Theorem 2.23 and Theorem 2.24 is very complicated, while the conditions of Theorem 5 are easily derived by using a generalization of Lyapunov stability theorem.

4 NSFD schemes for the model (4)

We propose NSFD schemes for the model (4) in the form

[TABLE]

Our task now is to determine the conditions for the function $\varphi(h)$ so that the scheme (21) preserves Properties $(P_{1})-(P_{3})$ of the model (4).

4.1 Monotone convergence

Theorem 4.1**.**

The scheme (21) preserves Property $(P_{1})$ of the model (4) if the function $\varphi(h)$ satisfies

[TABLE]

Proof.

Set $a_{k}:=S_{k}+I_{k}$, $b_{k}:=R_{k}+L_{k}$. In (21) adding consecutively the first equation with the second one, the third equation with the fourth one we obtain

[TABLE]

From here we have

[TABLE]

Since $\varphi(h)$ satisfies (22) then $1-\varphi f-\varphi g\in(0,1)$. From here it follows the proof of the theorem. ∎

4.2 Boundedness

Theorem 4.2**.**

The scheme (21) preserves Property $(P_{2})$ of the model (4) if the function $\varphi(h)$ satisfies

[TABLE]

Proof.

We prove the theorem by induction. First, adding the equations of (21) side-by-side we obtain

[TABLE]

Therefore, if $S_{k}+I_{k}+R_{k}+L_{k}=1$ then

[TABLE]

On the other hand, it is easy to transform the scheme (21) to the explicit form

[TABLE]

Hence, if $S_{k},I_{k},R_{k},L_{k}\geq 0$ and $\varphi(h)$ satisfy (24) then $S_{k+1},I_{k+1},R_{k+1},L_{k+1}\geq 0$. From this fact and (25) it follows the proof of the theorem. ∎

4.3 Stability properties

Since the discrete model (21) consists of four equations, the study of its stability via the set of eigenvalues of the Jacobian is very hard, even when we use the condition $S_{k}+I_{k}+R_{k}+L_{k}=1$ to reduce the system to a system of three equations. In the latter case the estimate of the Jacobian with the help of the Jury criterion remains complicated. As in Section 3 we shall use a generalization of Lyapunov stability theorem for getting NSFD schemes, which preserve Property $(P_{3})$ of the continuous model without complex computations. For this reason consider the function

[TABLE]

We shall show that the function $V$ satisfies all the conditions of [18, Theorem 3.2] on $D_{4}$. Clearly, $V$ is continuous on $D_{4}$, where $D_{4}$ is defined by (5). Moreover $V(E_{i}^{*})=0\;(i=1,2)$ and

[TABLE]

On the other hand we have

[TABLE]

Since Property $(P_{2})$ of the model is preserved, obviously $\Delta V(I_{k},S_{k},L_{k},R_{k})\leq 0$ for any $(I_{k},S_{k},L_{k},R_{k})\in D_{4}$. In this case we have

[TABLE]

In view of the fact that Property $(P_{1})$ of the model is preserved all the solutions of (21) are bounded.

Thus, the conditions (1), (2) and (4) of [18, Theorem 3.2] are satisfied. It remains only to show the $G^{*}$- local stability of the equilibrium points $E_{i}^{*}\,(i=1,2)$.

Notice that, for any initial conditions \big{(}I_{0},S_{0},L_{0},R_{0}\big{)} belonging to $G^{*}$ the solution \big{(}I_{k},S_{k},L_{k},R_{k}\big{)} also belongs $G^{*}$, that is, $I_{k}+S_{k}=f/(f+g)$ and $L_{k}+R_{k}=g/(f+g)$. Using this relation we reduce (21) to two equations depending on $I_{k},L_{k}$

[TABLE]

Put

[TABLE]

It is easy to see that the equilibrium points $E_{i}^{*}\,(i=1,2)$ are $G^{*}$-locally asymptotically stable if and only if $e_{i}^{*}(i=1,2)$ are locally asymptotically stable equilibrium points of the system (27). The Jacobian of (27) at $E^{*}=(I^{*},L^{*})$ are defined by

[TABLE]

For convenience in the future use we restate the Jury criterion for system of two difference equations.

Lemma 4.3**.**

The equilibrium point $E^{*}$ is locally asymptotically stable if and only if the Jacobian matrix $J(E^{*})$ satisfies [1, p.64, Theorem 2.10]

$\det(J(E^{*}))<1$** 2. 2.

$1-trace(J(E^{*}))+\det(J(E^{*}))>0$** 3. 3.

$1+trace(J(E^{*}))+\det(J(E^{*}))>0$**

**(i). The first equilibrium point **

Lemma 4.4**.**

If $c>0$ then

[TABLE]

where $c$ is defined by (8).

Proof.

Suppose the opposite, i.e., $\gamma\leq 0$. Then

[TABLE]

From (8) we obtain

[TABLE]

This is contrary to the assumption $c>0$. Thus, the lemma is proved. ∎

At the equilibrium point $e_{1}^{*}$, from (29) we have

[TABLE]

Theorem 4.5**.**

In the case $c>0$, if $\varphi(h)$ is a function satisfying

[TABLE]

then $e_{1}^{*}$ is a locally asymptotically stable equilibrium point of (27), i.e., $G^{*}$- locally asymptotically stable.

Proof.

From (31) we obtain

[TABLE]

Therefore, $\det(J(e_{1}^{*}))<1$ if and only if $\varphi<\gamma/c$. On the other hand

[TABLE]

Hence, if $\varphi<2/\gamma$ then $1+trace(J(e_{1}^{*}))+\det(J(e_{1}^{*}))>0$. Finally, we have

[TABLE]

Thus, all three conditions of Lemma 4.3 are satisfied. It implies that required to prove. ∎

**(ii). The second equilibrium point

**

Lemma 4.6**.**

Consider the polynomial of degree $n$ in the variable $\varphi$ with the coefficients $a_{i}\in\mathbb{R}$:

[TABLE]

Then, if $a_{0}>0$ then there exists a number $\varphi_{0}>0$ such that

[TABLE]

Proof.

The lemma is straightforward deduced from the definition of the limit. ∎

Lemma 4.7**.**

If $c<0$ then

[TABLE]

Proof.

From (8) we obtain

[TABLE]

So, if $c<0$ we have $b^{2}-4ac>\beta_{I}^{2}(f+e_{L})^{2}$. Having in mind this from (7) we obtain (33). The lemma is proved. ∎

The Jacobian at the second equilibrium defined from (29) has the form

[TABLE]

where

[TABLE]

Its trace has the form

[TABLE]

where

[TABLE]

Theorem 4.8**.**

If $c<0$ there exists a number $\varphi_{0}=\varphi_{0}(e_{I},e_{L},\beta_{I},\beta_{L},f,g)>0$ such that for any function $\varphi(h)$ satisfying

[TABLE]

the equilibrium point $e_{2}^{*}$ is $G^{*}$- asymptotically stable.

Proof.

We shall use Lemma 4.3 to prove the theorem. For this reason we show that all the conditions of the lemma are satisfied.

First from (34) we see that $\det(J(e_{2}^{*}))<1$ is equivalent to

[TABLE]

Now we show that

[TABLE]

Using the estimate (33) we have

[TABLE]

therefore from (40) we obtain $A_{1}>\beta_{L}I^{*}+f+e_{L}>0$. Hence, according to Lemma 4.6 there exists a number $\varphi_{1}>0$ such that for any $\varphi\in(0,\varphi_{1})$ there holds $\lambda_{1}(\varphi)>0$, i.e., $\det(J(E_{2}^{*}))<1$ for any $\varphi\in(0,\varphi_{1})$.

On the other hand, combining (34) and (36) we obtain $1-trace(J(e_{2}^{*}))+\det(J(e_{2}^{*}))>0$ if and only if

[TABLE]

We shall show that

[TABLE]

Making elementary calculations in combination with (7) and (8) we obtain

[TABLE]

Hence, according to Lemma 4.6 there exists a number $\varphi_{2}>0$ such that for any $\varphi\in(0,\varphi_{2})$ there holds $\lambda_{2}(\varphi)>0$, i.e., $1-trace(J(e_{2}^{*}))+\det(J(e_{2}^{*}))>0$ .

Finally, combining (34) and (36) we obtain $1+trace(J(e_{2}^{*}))+\det(J(e_{2}^{*}))>0$ if and only if

[TABLE]

Hence, according to Lemma 4.6 there exists a number $\varphi_{3}>0$ such that for any $\varphi\in(0,\varphi_{3})$ there holds $\lambda_{3}(\varphi)>0$, i.e., $1+trace(J(e_{2}^{*}))+\det(J(e_{2}^{*}))>0$ .

Set

[TABLE]

Then for any $\varphi<\varphi_{0}$ all the three conditions of Lemma 4.3 are satisfied. Thus, the theorem is proved. ∎

Summarizing the above results we obtain the following theorem of the nonstandard difference schemes preserving Properties $(P_{1})-(P_{3})$ of the model (4).

Theorem 4.9**.**

Consider the model (4). Then

In the case $c>0$, the scheme (21) preserves Properties $(P_{1})-(P_{3})$ of the model (4) if

[TABLE]

where $c$ and $\gamma$ are defined by (8) and (30), respectively. 2. 2.

In the case $c<0$, the scheme (21) preserves Properties $(P_{1})-(P_{3})$ of the model (4) if

[TABLE]

where $\varphi_{0}$ is defined as in Theorem 4.8.

Remark 3*.*

There are many ways for selecting the function $\varphi$ satisfying the conditions of Theorem 4.9, for example,

[TABLE]

5 Numerical simulations

Numerical examples presented in this section show that the obtained theoretical results of the NSFD preserving the properties of the metapopulation model are valid.

5.1 Numerical simulations for the model (2)

Very recently, in [8] we performed numerical simulations for the model (2) in particular cases. The result is that the standard difference schemes, for example, the typical four stages Runge-Kutta do not preserve the positivity of the solution. The two stage Runge-Kutta method and the explicit Euler method have solutions oscillating near equilibrium points with increasing amplitude. Therefore, if considering the models on large time intervals, the amplitude of oscillation will be large. The numerical solutions in this case cannot preserve properties of the model (2). In general, all standard difference schemes such as Runge-Kutta and Taylor methods preserve the properties of the continuous model for sufficiently small grid sizes.

Example 5.1**.**

Case $\mathcal{R}_{0}<1$.

Consider the model (2) with the parameters

[TABLE]

In this case $\mathcal{R}_{0}=0.4<1$. The model has two equilibrium points $P_{1}^{*}=(0.25,\,0)$ and $P_{2}^{*}$ = (0.625,   -0.375), among them $P_{1}^{*}$ is a globally asymptotically stable equilibrium point while $P_{2}^{*}$ is unstable equilibrium point.

For the scheme (9) - (10) we choose parameters $c_{i}$ $(i=1,2,5,6)$ satisfying Theorem 3.5, where

[TABLE]

The numerical solutions obtained by the schemes in this case are depicted in Figure 1, where each pair of blue and red curves corresponds to a solution $(x_{k},y_{k})$. In this case the properties of the model (2) are preserved.

Example 5.2**.**

Case $\mathcal{R}_{0}>1$.

Consider the model (2) with the parameters

[TABLE]

In this case $\mathcal{R}_{0}=3.75>1$. The model has two equilibrium points $P_{1}^{*}=(0.75,\,0)$ and $P_{2}^{*}$ = (0.2,   0.55), among them $P_{2}^{*}$ is a globally asymptotically stable equilibrium point and $P_{1}^{*}$ is unstable equilibrium point.

For the scheme (9)-(10) we choose the parameters $c_{i}$ $(i=1,2,5,6)$ satisfying Theorem 3.5, where

[TABLE]

The numerical solutions obtained by the schemes in this case are depicted in Figure 2. Clearly, the properties of the model (2) are preserved.

5.2 Numerical simulation for the model (4)

Example 5.3**.**

Case $c>0$.

Consider the model (4) with the parameters

[TABLE]

In this case $c=0.7031>0$. The model has two equilibrium points

[TABLE]

where $P_{1}^{*}$ is a locally asymptotically stable equilibrium point, and $P_{2}^{*}$ is unstable equilibrium point.

The numerical solutions obtained by the four stages Runge-Kutta method and the explicit Euler method are depicted in Figures 3 and 4, respectively. Clearly, the properties of the model are not preserved.

The numerical solutions obtained by these differences schemes are not positive, their boundedness is destroyed. The four stages Runge-Kutta generate spurious fixed points depending on the grid sizes. The explicit Euler method gives the numerical solutions oscillating near the equilibrium points. In general, the standard difference schemes preserve the properties of the continuous model only when the grid size is sufficiently small.

We shall use NSFD schemes (21) with the denominator function $\varphi(h)$ defined (43) in Theorem 4.9. In this case, since $c=0.7031$ and $\gamma=2$ we have $\varphi^{*}=0.8$. Therefore, we choose

[TABLE]

The numerical solutions of the scheme (21) in this case are depicted in Figures 5 and 6. All the properties of the model are preserved for any grid size $h$.

Example 5.4**.**

Case $c<0$.

Consider the model (4) with the parameters

[TABLE]

In this case we have $c=-0.2163<0$. The model has two equilibrium points

[TABLE]

where $P_{1}^{*}$ is an unstable equilibrium point and $P_{2}^{*}$ is a locally asymptotically stable equilibrium point.

As in the previous examples, the standard difference schemes such as Runge-Kutta or Taylor methods cannot preserve the properties of the continuous model. We shall use NSFD schemes (21) with the denominator function $\varphi(h)$ defined by (44) in Theorem 4.9. In this case we have

[TABLE]

The graphs of the functions $\lambda_{i}\;(i=1,2,3)$ are given in Figure 7. We can choose the number $\varphi_{0}$ in Theorem 4.9 equal to $1.5$. Using Theorem 4.9 we choose

[TABLE]

The numerical solutions obtained in this case are depicted in Figures 8 and 9. From the figures we see that all the properties of the model are preserved.

6 Conclusion

In this paper we have constructed the NSFD schemes preserving the properties of two metapopulation models. An essential result in the selection of the parameters of the schemes is that we used a generalization of Lyapunov stability theorem to ensure the stability properties of the the schemes. This approach is much simpler and more effective than the approach of studying Jacobians of the discrete systems which was used by ourselves in a previous work for the first metapopulation model [8]. In the future we shall develop this approach for constructing discrete models preserving the properties of continuous models for other problems.

Acknowledgments

This work is supported by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under the grant number 102.01-2014.20.