Dynamically Consistent NSFD Discretization of Some Productive-Destructive Population Models Satisfying Conservations Laws

Show more

1. Introduction

The dynamics of a population consisting of N interacting sub-populations can be described by the initial value problem

$\frac{\text{d}x}{\text{d}t}=f\left(x\right)$, $x\left({t}_{0}\right)={x}_{0}$, (1.1)

where $x={\left({x}^{1},{x}^{2},\cdots ,{x}^{n}\right)}^{\text{T}}$, $t\ge 0$, $f\left(x\right)={\left[{f}^{1},{f}^{2},\cdots ,{f}^{n}\right]}^{\text{T}}:{\mathbb{R}}^{n}\to {\mathbb{R}}^{n}$ is differentiable, and ${x}^{1}\left(0\right)\ge 0,{x}^{2}\left(0\right)\ge 0,\cdots ,{x}^{n}\left(0\right)\ge 0$.

The function $f\left(x\right)$ for such models is typically non-linear, so that the system (1.1) does not have analytic solutions, thus leaving numerical methods the only way that these models can be simulated; these can generally be expressed in the form

${D}_{h}\left({x}_{k}\right)={F}_{h}\left(f;{x}_{k}\right)$ (1.2)

where

${D}_{h}\left({x}_{k}\right)\approx {\frac{\text{d}x}{\text{d}t}|}_{t={t}_{k}},{x}_{k}\approx x\left({t}_{k}\right),h>0,{t}_{k}={t}_{o}+kh$, and ${F}_{h}\left(f;{x}_{k}\right)$ is a discrete approximation of $f\left(x\right)$.

It is well known that standard discretizations (1.2) of systems such as (1.1) result in un-reliable simulations [1] [2] [3]; in general, standard methods leads to discrete models with spurious solutions and instabilities that are not present in the original system. It has been observed ( [1] [4] ) that it is the lack of dynamic consistency with respect to some essential property of the continuous system, in the sense of Definition 1.1 below ( [4] ), that leads to numerical instabilities in the discrete approximations.

Definition 1.1 ( [4] ) Let the system (1.1) and its solution have a property P. The discrete system (1.2) as a model of (1.1) is said to be dynamically consistent with (1.1) with respect to the property P, if (1.2) and its solution also possess the property P.

To address the problem of spurious solutions and instabilities, the nonstandard finite difference (NSFD) method for discretizing continuous ODEs has been developed [1], whose hallmark is the production of numerical schemes that are dynamically consistent with the original ODE system. Therefore, salient features of the model, such as positivity of solutions, obeying conservation or sub-conservation laws, and other known structures of the original system play a prominent role in the NSFD methodology in the construction of dynamically consistent numerical models for such systems. However, a method that guarantees dynamic consistency with respect to one property may not guarantee dynamic consistency with respect to another. The question addressed in the present communication is therefore this: Are discrete models resulting from general productive-destructive (PD) NSFD construction methods applied to conservative PD systems dynamically consistent with respect to conservation laws?

A general NSFD construction method was advanced in [5] - [11] by Dimitrov, Koujaharov and Wood for PD systems with certain features in the interaction terms. Their method, which exploits the structure of the interaction terms to guarantee positivity, has been used as a basis for constructions of NSFD schemes for population models in [12] [13] [14] [15]. It is designed to produce numerical models that are dynamically consistent with the differential system with respect to positivity of solutions, equilibria and their linear stability, but not with respect to conservation laws. Another general method was proposed by Mickens in [16] and illustrated with applications to several population models by Mickens and Washington in [17] for the systematic construction of dynamically consistent models for systems obeying conservation or sub-conservation laws. The method uses a discretization of the conservation law to construct the denominator function and then exploits a positivity trick similar to the Koujaharov group method to complete the construction. The dynamic consistency with respect to conservation laws of discrete models resulting from the PD general method is investigated in the present article for conservative PD systems, and examples are presented to illustrate the theoretical considerations.

The remainder of this article is organized as follows. In Section 2, a synopsis of the NSFD methodology is given and relevant definitions are recalled concerning conservation laws and PD systems. Section 3 presents a summary of the respective general NSFD construction methods for PD and conservative systems, and the dynamic inconsistencies from the PD system method under consideration are discussed. Two examples, one satisfying a direct conservation law and the other a generalized conservation law, are presented in Section 4 to illustrate the inconsistencies and present dynamically consistent schemes using the conservation laws.

2. NSFD Methodology and Some Salient Features Used in NSFD Construction

In this section, a synopsis of the NSFD methodology is given and relevant definitions are recalled concerning some salient feature of systems used in population modeling. An important feature of the NSFD methodology is that the essential features of the continuous model must be incorporated into the discrete model since it is the failure of the discrete model to capture essential underlying features that lead to numerical instabilities in the resulting discrete solutions.

2.1. NSFD Methodology

There are two key requirements that form the foundation for the NSFD discretization of ordinary differential equations, which are well documented elsewhere [1] [3], and we only briefly outline them here.

The first NSFD requirement is that non-linear terms of the dependent functions must be modelled by non-local discrete representations on the computational grid. The following examples are typical in practice:

${u}^{2}\to {x}_{k}{x}_{k+1}$ ;

${u}^{2}\to \left(\frac{{x}_{k}+{x}_{k+1}}{2}\right){x}_{k}$ ;

${u}^{3}\to \frac{2{x}_{k}^{2}{x}_{k+1}^{2}}{{x}_{k}+{x}_{k+1}}$.

The second NSFD requirement is that the standard discrete model of the first derivative

$\frac{\text{d}u}{\text{d}t}\to \frac{{u}_{k+1}-{u}_{k}}{h}$,

must be replaced with a more general representation,

$\frac{\text{d}u}{\text{d}t}\to \frac{{x}_{k+1}-{x}_{k}}{\varphi \left(h\right)}$,

where the denominator function $\varphi \left(h\right)$ has the following properties:

1) $\varphi \left(h\right)=h+O\left({h}^{2}\right)$ as $h\to 0$ ;

2) $\varphi \left(h\right)$ is an increasing function of h;

3) $\varphi \left(h\right)$ may be a function of one or more parameters in the differential equations.

2.2. Some Salient Features Used in the NSFD Methodology

This section briefly outlines some salient features of ordinary differential equations systems that may be used in the development of general NSFD construction methods for certain classes of the system (1.1). A fundamental requirement for system (1.1) as a model for the study of populations is that all solutions are positive [11] [17], and hence the next definition.

Definition 2.1. ( [11] [17] ) The solution $x\left(t\right)$ of (1.1) is said to satisfy a positivity condition if (and only if)

${x}^{i}\left(0\right)\ge 0\Rightarrow {x}^{i}\left(t\right)\ge 0$ for $t>0,1\le i\le n$. (2.1)

This condition guarantees that populations start off non-negative and remain non-negative, since negative populations are generally meaningless, and has been used in all general construction methods. Systems with features contained in the following two definitions arise often enough in population and physical modeling that general discretization methods have been developed for them utilizing these features.

Definition 2.2 ( [11] ) The system (1.1) is said to be productive-destructive (PD) if it can be expressed in the form

$\frac{\text{d}x}{\text{d}t}=P\left(x\right)-D\left(x\right),{x}_{0}\in {\mathbb{R}}_{+}^{n}$. (2.2)

where $P\left(x\right)\ge 0$, and $D\left(x\right)\ge 0$, are positive input and output functions.

Definition 2.3 ( [16] [17] ) The system (1.1) is said to obey a direct or generalized conservation law, respectively, if it satisfies the condition

$\frac{\text{d}N}{\text{d}t}=0,$ (2.3a)

or

$\frac{\text{d}N}{\text{d}t}\le a-bN$ (2.3b)

where $a\ge 0,b\ge 0$ are constants and $N={x}^{1}+{x}^{2}+\cdots +{x}^{n}$.

To correctly represent the same dynamics as the original system (1.1), a dynamically consistent numerical method (1.2) of (1.1) must, according to Definition 1.1, therefore be required to satisfy discrete conditions analogous to those satisfied by (1.1). Our interest in the present article is the construction of NSFD models of PD systems that are dynamically consistent with respect to positivity and a conservation law as per Definition 2.3. As has been summarized in [16] [17], the resulting numerical schemes must satisfy the following conditions:

1) The positivity condition must be satisfied for each dependent variable in the discretization;

2) If a term appears in more than one equation of the system, then its discrete representation must be the same in all appearances;

3) The discrete model of the system must exactly satisfy the same conservation law, in discrete form, as the continuous system;

4) The denominator function for the discrete first derivatives must have explicit functional form and must be the same for all equations in the system.

3. General NSFD Construction Methods

In this section we present a summary of the general NSFD construction method for PD developed in [5] - [11] and for conservative (or semi-conservative) systems developed in [16] [17], and discuss the inconsistency with the conservation laws of the NSDF schemes from the PD systems method.

3.1. PD Systems NSFD Construction

In [9] [10] [11] a general NSFD method was developed for a class of PD systems that satisfy the positivity as well as the following condition:

$\frac{\text{d}{x}^{i}}{\text{d}t}={P}^{i}\left(x\right)-{F}^{i}\left(x\right){x}^{i},1\le i\le n$, (3.1)

where the functions ${P}^{i}\left(x\right),{F}^{i}\left(x\right)$ are non-negative.

That general construction method for (1.1) yields the numerical model

$\frac{{x}_{k+1}^{i}-{x}_{k}^{i}}{\varphi \left(h\right)}={f}^{i}\left({x}_{k}\right){\omega}_{k}^{i},1\le i\le n$, (3.2a)

where

${\omega}_{k}^{i}=\{\begin{array}{l}1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}{f}^{i}\left({x}_{k}\right)\ge 0\\ \frac{{x}_{k+1}^{i}}{{x}_{k}^{i}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}{f}^{i}\left({x}_{k}\right)<0\end{array}$ (3.2b)

For a PD system of the form (3.1), this yields the numerical model

$\frac{{x}_{k+1}^{i}-{x}_{k}^{i}}{\varphi \left(h\right)}={P}^{i}\left({x}_{k}\right)-{F}^{i}\left({x}_{k}\right){x}_{k+1}^{i},1\le i\le n$ (3.3)

The model (3.3), while it satisfies the positivity condition, may fail to satisfy the underlying conservation law if some of the terms in ${F}^{i}\left(x\right){x}^{i}$ are also present in ${P}^{j}\left(x\right),i\ne j$ since such a term will have non-identical discretizations.

If a system satisfies the generalized conservation law (2.3b), that is, $\frac{\text{d}N}{\text{d}t}\le a-bN$, such as in [10] (Example 1 in Section 4), then the discretization of said system must satisfy the discrete conservation analog of (2.3b), which is

$\frac{{N}_{k+1}-{N}_{k}}{{\varphi}_{1}\left(h,b\right)}\le a-b{N}_{k}$, ${\varphi}_{1}\left(h,b\right)=\frac{1}{b}\left(1-{\text{e}}^{-bh}\right)$ (3.4a)

or

$\frac{{N}_{k+1}-{N}_{k}}{{\varphi}_{2}\left(h,b\right)}\le a-b{N}_{k+1}$, ${\varphi}_{2}\left(h,b\right)=\frac{1}{b}\left({\text{e}}^{bh}-1\right)$ (3.4b)

However, because of the non-identical discretization of some terms of the form ${c}_{i}{x}^{i}{x}^{j}$ which appear in both the positive and negative components of (1.1), instead of resulting in (3.4), addition of the components of (1.1) result in an expression of the form

$\frac{{N}_{k+1}-{N}_{k}}{\varphi \left(h,b\right)}\le a+{\displaystyle {\sum}_{i=1}^{m}{c}_{i}\left({x}_{k}^{i}-{x}_{k+1}^{i}\right)}-b{N}_{k+1},m<n$. (3.5)

Since the middle term in (3.5) can be positive, that is,

$a+{\displaystyle {\sum}_{i=1}^{m}{c}_{i}}\left({x}_{k}^{i}-{x}_{k+1}^{i}\right)>0$ (3.6)

is possible, then in this case Equations (3.4) do not necessarily hold. The discretization (3.3) therefore does not necessarily result in (3.4), and thus is not dynamically consistent with respect to the conservation law (2.3b).

If a PD system satisfies the exact conservation law (2.3a), this may be because the equation takes the form

$\frac{\text{d}N}{\text{d}t}=aN-bN=\left(a-b\right)=0$

that is, $a=b$, which is a non-conditional identity, such as in [14] (Example 2 in Section 4). In this case, the discrete conservation analog of (2.3a) is

$\frac{{N}_{k+1}-{N}_{k}}{{\varphi}_{1}\left(h\right)}=a{N}_{k}-b{N}_{k}=\left(a-b\right){N}_{k}=0$ (3.7a)

or

$\frac{{N}_{k+1}-{N}_{k}}{{\varphi}_{2}\left(h\right)}=a{N}_{k+1}-b{N}_{k+1}=\left(a-b\right){N}_{k+1}=0$ (3.7b)

for suitable denominator functions ${\varphi}_{1}\left(h\right),{\varphi}_{2}\left(h\right)$. However, in this case, the non-identical discretization of some terms of the form ${c}_{i}{x}^{i}{x}^{j}$ which appear in both the positive and negative components of (1.1), instead of resulting in (3.4), addition of the components of the system result in an expression of the form

$\begin{array}{c}\frac{{N}_{k+1}-{N}_{k}}{\varphi \left(h,b\right)}=a{N}_{k+1}+{\displaystyle {\sum}_{i=1}^{m}{c}_{i}}\left({x}_{k}^{i}-{x}_{k+1}^{i}\right)-b{N}_{k+1}\\ ={\displaystyle {\sum}_{i=1}^{m}{c}_{i}}\left({x}_{k}^{i}-{x}_{k+1}^{i}\right)+\left(a-b\right){N}_{k+1},\end{array}$ (3.8)

It can be seen that the righthand side of (3.8) is only conditionally equal to zero, and hence does not preserve (2.3a), which is unconditional. Since the middle term in (3.8) can be positive, that is, it is possible that

$a+{\displaystyle {\sum}_{i=1}^{m}{c}_{i}}\left({x}_{k}^{i}-{x}_{k+1}^{i}\right)>0$, (3.9)

the Equations (3.7) do not necessarily hold. The discretization (3.3) can therefore fail to be dynamically consistent with respect to the conservation law (2.3b) since it does not necessarily result in (3.7); that is, it does not preserve the conservation law.

3.2. Conservative Systems NSFD Construction

In [16] [17] a general NSFD method was developed for a class of systems satisfying the conservation laws (2.3a) and (2.3b). For systems satisfying (2.3b), the Mickens-Washington begins with the exact solution (3.4a) or (3.4b) to obtain the denominator function and then and uses positivity to complete the construction. First, each equation is expressed in the following form, which is a special case of (3.1):

$\frac{\text{d}{x}^{i}}{\text{d}t}={P}^{i}\left(x\right)-{f}_{1}^{i}\left(x\right)-{f}_{2}^{i}\left(x\right){x}^{i}-b{x}^{i},1\le i\le n$, (3.10)

where $b>0$ and each of ${P}^{i}\left(x\right),{f}_{1}^{i}\left(x\right),{f}_{2}^{i}\left(x\right)$ are non-negative.

The resulting general discretization is

$\begin{array}{l}\frac{{x}_{k+1}^{i}-{x}_{k}^{i}}{\varphi}={P}^{i}\left({x}_{k},{x}_{k+1}\right)-{f}_{1}^{i}\left({x}_{k},{x}_{k+1}\right)\frac{{x}_{k+1}^{i}}{{x}_{k}^{i}}-{f}_{2}^{i}\left({x}_{k}\right){x}_{k+1}^{i}-b{x}^{i},\\ \varphi =\frac{1}{b}\left(1-{\text{e}}^{-bh}\right),\end{array}$ (3.11)

which is consistent with (3.4a); the ${x}_{k+1}$ in ${P}^{i}\left({x}_{k},{x}_{k+1}\right)$ and ${f}_{1}^{i}\left({x}_{k},{x}_{k+1}\right)$ accounts for the terms which also appear in some ${f}_{2}^{j}\left({x}_{k}\right){x}_{k+1}^{j},j\ne i$, and thus require the same discretization.

For systems satisfying (2.3a), the Mickens-Washington begins with the exact solution of (3.7a), with a denominator function chosen to ensure that the conservation law is obeyed, yielding a construction similar to (3.11).

4. Examples

In this section, two examples are presented to illustrate the pitfalls which may lead to dynamic inconsistencies for conservative systems in the PD system NSFD general construction method of Dimitrov and co-workers. Dynamically consistent NSFD schemes are presented using the Mickens-Washington general construction method for conservative systems, which begins with the exact discretization of the conservation or a sub-conservation law to obtain the discrete denominator for the system and is completed by applying positivity.

*Example *1

The first example is a model for the study the role of incubation in AIDS which was used in [11] to illustrate the PD system NSFD construction method.

$\frac{\text{d}{S}^{p}}{\text{d}t}=c\text{\Lambda}-\left(1-{\alpha}_{s}\right)\beta \frac{I{S}^{p}}{N}-\mu {S}^{p}$ (4.1a)

$\frac{\text{d}S}{\text{d}t}=\left(1-c\right)\text{\Lambda}-\beta \frac{IS}{N}-\mu S$ (4.1b)

$\frac{\text{d}I}{\text{d}t}=\beta \frac{IS}{N}+\left(1-{\alpha}_{s}\right)\beta \frac{I{S}^{p}}{N}-\left(\mu +d\right)I$ (4.1c)

where $N={S}^{p}+S+I$.

The system (4.1) satisfies the generalized conservation law (2.3b) since total population $N={S}^{p}+S+I$ satisfies

$\frac{\text{d}N}{\text{d}t}=c\text{\Lambda}-dI-\mu N\le c\text{\Lambda}-\mu N$. (4.2)

The numerical model obtained in [17] using (3.2) is

$\frac{{S}_{k+1}^{p}-{S}_{k}^{p}}{\varphi \left(h\right)}=c\text{\Lambda}-\left(1-{\alpha}_{s}\right)\beta \frac{{I}_{k}{S}_{k+1}^{p}}{{N}_{k}}-\mu {S}_{k+1}^{p}$ (4.3a)

$\frac{{S}_{k+1}-{S}_{k}}{\varphi \left(h\right)}=\left(1-c\right)\text{\Lambda}-\beta \frac{{I}_{k}{S}_{k+1}}{{N}_{k}}-\mu {S}_{k+1}$ (4.3b)

$\frac{{I}_{k+1}-{I}_{k}}{\varphi \left(h\right)}=\beta \frac{{I}_{k}{S}_{k}}{{N}_{k}}+\left(1-{\alpha}_{s}\right)\beta \frac{{I}_{k}{S}_{k}^{p}}{{N}_{k}}-\left(\mu +d\right){I}_{k+1}$ (4.3c)

Because of the non-identical discretization of the terms $\frac{I{S}^{p}}{N}$ in (4.3a), (4.3c) and of $\frac{IS}{N}$ in (4.3b), (4.3c), adding (4.3a) through (4.3c) yields

$\frac{{N}_{k+1}-{N}_{k}}{\varphi \left(h\right)}=\Lambda -\beta \frac{{I}_{k}}{{N}_{k}}\left(\left({S}_{k+1}^{p}-{S}_{k}^{p}\right)+\left({S}_{k+1}-{S}_{k}\right)\right)-d{I}_{k+1}-\mu {N}_{k+1}$ (4.4a)

so that condition (3.4b), the discrete equivalent of (4.2), is true if and only if

$\left({S}_{k+1}^{p}-{S}_{k}^{p}\right)+\left({S}_{k+1}-{S}_{k}\right)+d{I}_{k+1}\ge 0$. (4.4b)

Therefore, the conservation law (4.2), an unconditional identity, is preserved only conditionally by the scheme (4.3), which is not consistent with the original system.

Applied to the system (4.1), the Mickens-Washington method (3.11) results in

$\frac{{S}_{k+1}^{p}-{S}_{k}^{p}}{\varphi \left(h\right)}=c\Lambda -\left(1-{\alpha}_{s}\right)\beta \frac{{I}_{k}{S}_{k+1}^{p}}{{N}_{k}}-\mu {S}_{k}^{p}$ (4.5a)

$\frac{{S}_{k+1}-{S}_{k}}{\varphi \left(h\right)}=\left(1-c\right)\Lambda -\beta \frac{{I}_{k}{S}_{k+1}}{{N}_{k}}-\mu {S}_{k}$ (4.5b)

$\frac{{I}_{k+1}-{I}_{k}}{\varphi \left(h\right)}=\beta \frac{{I}_{k}{S}_{k+1}}{{N}_{k}}+\left(1-{\alpha}_{s}\right)\beta \frac{{I}_{k}{S}_{k+1}^{p}}{{N}_{k}}-d{I}_{k+1}-\mu {I}_{k}$ (4.5c)

where $\varphi =\varphi \left(h,\tau \right)=\frac{1}{\mu}\left(1-{\text{e}}^{-\mu h}\right)$.

The denominator function $\varphi $ in (4.5) is determined from the exact discretization of the conservation law $\frac{\text{d}N}{\text{d}t}=c\text{\Lambda}-\mu N$, which is $\frac{{N}_{k+1}-{N}_{k}}{\varphi}=c\text{\Lambda}-\mu {N}_{k}$.

The system (4.5) satisfies the positivity requirement and unconditionally satisfies

$\frac{{N}_{k+1}-{N}_{k}}{\varphi \left(h\right)}=\text{\Lambda}-d{I}_{k+1}-\mu {N}_{k}\le \text{\Lambda}-\mu {N}_{k}$ (4.6)

which preserves the conservation law (4.2) and is therefore dynamically consistent with respect to the conservation law of the original system. Moreover, while implicit in construction, the model (4.5) can be solved for explicit implementation to obtain in the order ${S}_{k+1}^{p},{S}_{k+1},{I}_{k+1}$.

*Example *2

The next example is a model studied in [6] to illustrate the PD system NSFD construction method.

$\frac{\text{d}S}{\text{d}t}=\mu N-\beta \frac{SI}{N}-\left(\mu +\tau \right)S+cI+\delta V$ (4.7a)

$\frac{\text{d}I}{\text{d}t}=\beta \frac{IS}{N}-\left(\mu +c\right)I$ (4.7b)

$\frac{\text{d}V}{\text{d}t}=\tau S-\left(\mu +\delta \right)V$ (4.7c)

The total population $N=S+I+V$ in model (4.7) satisfies the direct conservation law

$\frac{\text{d}N}{\text{d}t}=\mu N-\mu N=0$, (4.8)

The PD method of the Dimitrov group results in

$\frac{{S}_{k+1}-{S}_{k}}{\varphi}=\mu {N}_{k}-\beta \frac{{I}_{k}{S}_{k+1}}{{N}_{k}}-\mu {S}_{k+1}-\tau {S}_{k+1}+c{I}_{k}+\delta {V}_{k}$ (4.9a)

$\frac{{I}_{k+1}-{I}_{k}}{\varphi}=\beta \frac{{I}_{k}{S}_{k}}{{N}_{k}}-\left(c+\mu \right){I}_{k+1}$ (4.9b)

$\frac{{V}_{k+1}-{V}_{k}}{\varphi}=\tau {S}_{k}-\left(\delta +\mu \right){V}_{k+1}$. (4.9c)

The discrete representations $\mu N\to \mu {N}_{k},cI\to c{I}_{k}$ in Equation (4.9a) and

$\mu S\to \mu {S}_{k+1}$, $\left(c+\mu \right)I\to \left(c+\mu \right){I}_{k+1}$, $\mu V\to \mu {V}_{k+1}$ in Equations (4.9a)-(4.9c) result in

$\frac{{N}_{k+1}-{N}_{k}}{\varphi}=f\left({S}_{k},{V}_{k},{I}_{k},{S}_{k},{V}_{k+1},{I}_{k+1}\right)+\mu \left({N}_{k}-{N}_{k+1}\right)$ (4.10)

which is only conditionally equal to zero; this is not consistent with the conservation law (4.8), which is an unconditional identity.

The following numerical model of (4.7) is the result of the Mickens-Washington procedure for such systems:

$\frac{{S}_{k+1}-{S}_{k}}{\varphi}=\mu \left({I}_{k+1}+{V}_{k+1}\right)-\beta \frac{{I}_{k}{S}_{k+1}}{\stackrel{\u02dc}{{N}_{k}}}-\tau {S}_{k}+c{I}_{k+1}+\delta {V}_{k+1}$ (4.11a)

$\frac{{I}_{k+1}-{I}_{k}}{\varphi}=\beta \frac{{I}_{k}{S}_{k+1}}{\stackrel{\u02dc}{{N}_{k}}}-\left(c+\mu \right){I}_{k+1}$ (4.11b)

$\frac{{V}_{k+1}-{V}_{k}}{\varphi}=\tau {S}_{k}-\left(\delta +\mu \right){V}_{k+1}$, (4.11c)

where $\stackrel{\u02dc}{{N}_{k}}={S}_{k}+{I}_{k}+{V}_{k+1}$

and $\varphi =\varphi \left(h,\tau \right)=\frac{1}{\tau}\left(1-{\text{e}}^{-\tau h}\right)$ (4.12)

The denominator function in (4.11) is determined from noting that when $I=V=0$, then

$\frac{\text{d}S}{\text{d}t}=\mu N-\left(\mu +\tau \right)S=-\tau S$, whose exact solution $S\left(t\right)=S\left(0\right)\mathrm{exp}\left(-\tau t\right)$ has the exact finite difference scheme $\frac{{S}_{k+1}-{S}_{k}}{\varphi}=-\tau {S}_{k}$, with $\varphi $ given by (4.12).

The model (4.11) unconditionally satisfies the discrete equivalent of conservation law (4.8) as well as the positivity requirement; while implicit in construction, it can be solved for explicit implementation to obtain the variables in the order ${V}_{k+1},{I}_{k+1},{S}_{k+1}$.

5. Conclusions

The dynamic consistency with respect to conservation laws of discrete models of conservative PD systems resulting from an NSFD general construction method for PD systems has been investigated. Two general methods of constructing NSFD numerical models are recalled, where one method is designed for PD systems and the other for systems satisfying conservation laws. It is shown that the NSFD method for general PD systems does not produce numerical models that are dynamically consistent with respect to conservation laws. This is significant since the failure of the discrete model to exactly satisfy the conservation laws of the original continuous system leads to discrete solutions that do not possess the qualitative features of the continuous model, which in turn may lead to numerical instabilities. As it fails to preserve conservation laws of the systems it models, the PD method considered here may therefore be unreliable as a basis for the general construction of numerical models of systems of differential equations satisfying such conservation laws, since such failure may lead to numerical instabilities. In conclusion therefore, if additional structure is imposed on PD dynamical systems in the form of conservation laws, then the NSFD methods designed for PD systems might not preserve the conservation laws, and alternative methods should be used in such cases.

Two examples of important population models are presented to illustrate the theory discussed, where one example satisfies the direct conservation law and the other a generalized conservation law. For these examples, the numerical model based on the PD system NSFD general construction method is presented and the dynamic inconsistency with the underlying conservation law is made explicit. In each case an NSFD model is also presented which is dynamically consistent with respect to positivity of solutions and the applicable conservation law.

References

[1] Mickens, R.E. (1994) Nonstandard Finite Difference Models of Differential Equations. World Scientific, Singapore. https://doi.org/10.1142/2081

[2] Mickens, R.E. (2005) Advances in the Applications of Nonstandard Finite Difference Schemes. World Scientific, Singapore. https://doi.org/10.1142/5884

[3] Mickens, R.E. (2005) Nonstandard Finite Difference Schemes for Differential Equations. Journal of Difference Equations and Applications, 8, 823-847.

https://doi.org/10.1080/1023619021000000807

[4] Mickens, R.E. (2005) Dynamic Consistency: A Fundamental Principle for Constructing NSFD Schemes for Differential Equations. Journal of Difference Equations and Applications, 11, 645-653. https://doi.org/10.1080/10236190412331334527

[5] Dimitrov, D.T. and Kojouharov, H.V. (2006) Positive and Elementary Stable Nonstandard Numerical Methods with Applications to Predator-Prey Models. Journal of Computational and Applied Mathematics, 189, 98-108.
https://doi.org/10.1016/j.cam.2005.04.003

[6] Dimitrov D. and Kojouharov, H. (2007) Stability-Preserving Finite-Difference Methods for General Multi-Dimensional Autonomous Dynamical Systems. International Journal of Numerical Analysis and Modeling, 4, 280-290.

[7] Dimitrov, D.T. and Kojouharov, H.V. (2008) Nonstandard Finite-Difference Methods for Predator-Prey Models with General Functional Response. Mathematics and Computers in Simulation, 78, 1-11.
https://doi.org/10.1016/j.matcom.2007.05.001

[8] Dimitrov, D.T. and Kojouharov, H.V. (2011) Dynamically Consistent Numerical Methods for General Productive-Destructive Systems. Journal of Difference Equations and Applications, 17, 1721-1736. https://doi.org/10.1080/10236191003781947

[9] Wood, D.T., Dimitrov, D.T. and Kojouharov, H.V. (2015) A Nonstandard Finite Difference Method for n-Dimensional Productive-Destructive Systems. Journal of Difference Equations and Applications, 21, 240-254.
https://doi.org/10.1080/10236198.2014.997228

[10] Wood, D. and Kojouharov, H. (2015) A Class of Nonstandard Numerical Methods for Autonomous Dynamical Systems. Applied Mathematics Letters, 50, 78-82.

https://doi.org/10.1016/j.aml.2015.06.008

[11] Wood, D. (2015) Advancement and Applications of Nonstandard Finite Difference Methods. Ph.D. Thesis, University of Texas at Arlington.

[12] Bairagi, N. and Biswas, M. (2016) A Predator-Prey Model with Beddington-DeAngelis Functional Response: A Non-Standard Finite Difference Method. Journal of Difference Equations and Applications, 22, 581-593.

[13] Biswas, M. and Bairagi, N. (2017) Discretization of an Eco-Epidemiological Model and Its Dynamic Consistency. Journal of Difference Equations and Applications, 23, 860-877. https://doi.org/10.1080/10236198.2017.1304544

[14] Quang, D.A. and Hoang, M.T. (2017) Nonstandard Finite Difference Schemes for a General Predator-Prey System. Journal of Computational Science, 36, 101015.
https://arxiv.org/abs/1701.05663v1

https://doi.org/10.1016/j.jocs.2019.07.002

[15] Saha, P., Bairagi, N. and Biswas, M. (2019) On the Dynamic Consistency of a Discrete Predator-Prey Model. arXiv preprint arXiv:1906.02513.

[16] Mickens, R.E. (2007) Numerical Integration of Population Models Satisfying Conservation Laws: NSFD Methods. Journal of Biological Dynamics, 1, 427-436.

https://doi.org/10.1080/17513750701605598

[17] Mickens, R.E. and Washington, T.M. (2013) NSFD Discretizations of Interacting Population Models Satisfying Conservation Laws. Computers and Mathematics with Applications, 66, 2307-2316. https://doi.org/10.1016/j.camwa.2013.06.011