Numerical integration of population models satisfying conservation laws: NSFD methods

Abstract

Population models arising in ecology, epidemiology and mathematical biology may involve a conservation law, i.e. the total population is constant. In addition to these cases, other situations may occur for which the total population, asymptotically in time, approach a constant value. Since it is rarely the situation that the equations of motion can be analytically solved to obtain exact solutions, it follows that numerical techniques are needed to provide solutions. However, numerical procedures are only valid if they can reproduce fundamental properties of the differential equations modeling the phenomena of interest. We show that for population models, involving a dynamical conservation law the use of nonstandard finite difference (NSFD) methods allows the construction of discretization schemes such that they are dynamically consistent (DC) with the original differential equations. The paper will briefly discuss the NSFD methodology, the concept of DC, and illustrate their application to specific problems for population models.

Keywords:

• Population models
• Conservation laws
• Numerical integration
• NSFD methods

1. Introduction

Interacting population models are used to analyze, understand, and predict the dynamics of phenomena in the natural and engineering sciences. In particular, they provide useful descriptions of systems in ecology, epidemiology and the biosciences. For the case where space effects are not taken into consideration, interacting population models are generally represented as systems of coupled nonlinear ordinary differential equations (ODE) 1–4, i.e.

where for a given system the x _i (t) are the individual population numbers or densities and f _i (X) are given functions, . An explicit example is the following SIR model which includes vital dynamics 1 4

where (S, I, R) stand, respectively, for the susceptible, infective, and removed populations; k is the birth rate; β is related to the rate that a susceptible individual becomes infective; and γ is the rate an average individual leaves from being infective. The μ parameter gives the death rate, assumed here to be the same for all three populations. Adding the three components of equation (3) gives

For the usual situation where S(0) > 0, I(0) > 0, and R(0) = 0, we have and the solution to equation (4) is

From this we conclude that the total population goes asymptotically in time to the value

Another way of expressing this situation is to say that for any finite t, the total population is changing, but eventually it goes to the constant value given in equation (6).

In the work to follow, we present interacting population systems for which the total population satisfies a scalar, first-order ODE which takes the form

where, see equation (2),

Under this condition, we define equation (7) as the conservation law associated with the system of coupled ODEs given by equations (1) and (2). The remarkable feature of the systems that we have studied to date is that the conservation law for all of them take the form

A valid mathematical model for interacting populations must satisfy a condition of positivity, i.e. the various subpopulations must never take on negative values 4,

Several theorems exist on the conditions that the functions f _i (X) must satisfy to have boundedness and positivity. A very good discussion is given by Thieme 4 in section 1 of Appendix A.

Our interest in examining interacting population models, as given by equations (1) and (2), is to formulate discretized versions that can be used for their numerical integration. Except for a small class of models for which exact analytical solutions can be calculated 5, essentially all of the existing models must be numerically integrated if we are to understand the full dynamics of these models for particular ranges of the parameters and/or initial values 6 7. A measure of the validity of a numerical integration scheme, for the cases examined in this work, is whether important constraint conditions and/or equations are satisfied by the numerical solutions. We show that our particular nonstandard finite difference (NSFD) discretizations 6 always exactly satisfy the associated conservation law. This is in contrast to the use of most standard methods for which the conservation law only holds approximately 7. Thus, the main purpose of this paper is to show that dynamical consistent finite difference schemes can be constructed for interacting population models having an associated conservation law.

In section 2, we give a brief general summary of the NSFD methodology 6, 8. A number of interacting population models are considered in section 3 where we demonstrate the power of NSFD schemes to give dynamically consistent discrete schemes. Finally, in the last section we discuss our results and mention some related issues.

2. NSFD methodology

We begin with the definition of dynamic consistency 9 for a single (scalar) ODE.

Definition

Let

be a finite difference model for the equation

where λ represents m-parameters characterizing the system, i.e. ; h = Δ t and t _k = hk; and u _k is the approximation to u(t _k ). Let the ODE and/or its solutions have some property P. Then the finite difference scheme, equation (11), is dynamic consistent (DC) with the ODE, with respect to property P, if it and/or its solutions also have property P.

Note that the lack of DC with respect to some essential property generally gives rise to numerical instabilities in the computed solutions of the finite difference scheme 6 9. The importance of DC is to use it to aid in the construction of NSFD schemes 8 9. In this regard, the NSFD methodology is based on the two fundamental principles:

• (i) The replacement of the discrete first-order derivative by a structure having the form 6

  where depends on both the parameters, λ, and the time step-size and φ has the property

• (ii) The general replacement of linear and nonlinear functions of u(t) by ‘nonlocal’ discrete representations. Particular examples, for first-order ODEs, are 6 9

The function is called the denominator function 6. Recently, we have created a procedure for calculating denominator functions for positivity preserving systems 10. The basic idea is to consider equations (1) and (2). If in the ith ODE, we have a linear term involving the ith x-variable, i.e.

then in the NSFD scheme for the differential equation, the denominator function is given by the expression

The full details and related explanations on the NSFD methodology are given in Mickens 6 9 10. Applications of this numerical integration technique to a broad range of problems in the natural, engineering, and bio-medical sciences is illustrated in Mickens 11.

3. Applications

The following comments are of general applicability for the construction of NSFD schemes for population models satisfying a conservation law.

• (a) The requirement that the positivity condition holds determines the nonlocal structure of the discrete representation of the first ODE considered in any system of ODEs 6 8 11.

• (b) If the same term occurs in more than one differential equation, then it must be modeled discretely the same way in all of the equations. Thus, terms appearing in earlier ODEs dictate the discrete form of these same terms occurring in later equations 6 10.

• (c) Since a conservation law exists, thus connecting all of the dependent variables together as given by equations (7) and (8), the denominator functions for all of the finite difference schemes must be the same function 10.

• (d) While the discrete schemes for the system of ODEs are NSFD representations and therefore are approximations to the theoretically existing exact schemes, the conservation law discretization must be an exact finite difference scheme 6.

We now present and illustrate the construction of NSFD schemes for three systems of equations relevant to the biosciences.

3.1 SIR model

Consider an SIR model with vital dynamics, i.e. containing births and deaths 12,

Let and add these equations to obtain

First, observe that β IS appears in the first and second ODEs. The negative sign of this term in equation (18a) implies that we replace it by a NSFD discretization such as

so as to have positivity for S _k+1. Likewise, the μ S term in this equation is replaced by

Also, from the fact that equation (18a) has, on the right-hand side, a term, −μ S, it follows from the discussion leading to equation (17) that the denominator function must be

Putting all of this together gives for equation (18a), the following NSFD scheme

Solving for S _k+1, we obtain

Both equations (18b) and (18c) contain the term γ I with opposite signs. For equation (18b) the positivity condition requires that (−γ I) be replaced by 9. Thus, the NSFD scheme for equation (18b) is

where again, because of positivity, the (−μ I) term is discretized by the expression (−μ I _k+1). Note that the β IS term discrete representation is determined by what was used in equation (23). Solving for I _k+1 gives

Without presenting the details, it is clear that equation (18c) must have the discrete representation

or

Adding equations (23), (25) and (27), and using

gives

which is the exact finite difference scheme 6 9 for equation (19).

In summary, inspection of equations (24), (26), (28) and (30) shows that we have constructed a NSFD scheme for the SIR model having the following properties:

• (i) positivity; S _k ≥0, I _k ≥0, , I _k+1≥0, R _k+1≥0;

• (ii) bounded solutions for since is bounded;

• (iii) an explicit representation for the denominator function φ(μ, h).

Observe that the populations must be calculated in a particular order. First, S _k+1 is determined from S _k and I _k , see equation (24); second, I _k+1 is obtained using knowledge of S _k+1 and I _k , see equation (26); finally, R _k+1 is found from I _k+1 and R _k . This sequential form of calculation is a generic feature of NSFD methods 6 9 10.

3.2 Brauer-van den Driessche SIR model

This model takes the form (To Reluga, personal communication, November 2005):

where are positive parameters. Defining M to be

then M(t) satisfies the following first-order differential equation

which is the conservation law for this SIR model.

A trivial calculation shows that the exact finite difference scheme for equation (33) is

where the denominator function is

Following the rules, indicated at the beginning of this section, the following NSFD scheme is obtained for equations (31a), (31b), and (31c):

Observe that equations (36a), (36b) and (36c) are three linear equations in terms of (S _k+1, I _k+1, R _k+1). Solving for these quantities gives

and from equation (34)

The expressions for are

with

Putting all this together, it can be concluded, just as for the previous SIR case in section 3.1, that all solutions for are nonnegative and bounded, and that the discrete version of the conservation law is exactly satisfied.

3.3 Spatial spread of rabies

In his book, Mathematical Biology, Murray introduces a SIR model having diffusion taking place among the infective population; see section 20.3 of Murray 13. The three equations take the form

where (r, a, D) are positive parameters. Note that partial time derivatives appear in the three equations since all dependent variables now depend on x and t, i.e. S = S(x, t), etc. Adding these equations gives

Note that in the absence of diffusion the total population is constant, i.e. for this case, a conservation law exists

The NSFD scheme for equations (41a), (41b) and (41c) is given by the following expressions:

where the following notation is used

and and have the properties 6 8

For the time being, we do not specify the denominator functions φ(z) and ψ (z).

If equations (44a), (44b) and (44c) are solved, respectively, for , , and , the following expressions are obtained:

where is defined to be

Inspection of equations (47a), (47b), (47c) shows that can be determined from and ; can be calculated from knowledge of , , , and ; and, finally, follows from knowing and . Thus, the calculation of the solutions is done sequentially. Note that for any values

an immediate consequence is and . However, to ensure that , we require 9

A choice often selected is to use

where γ is a nonnegative number. Making this substitution in equation (47b) gives a positivity-preserving scheme for the calculation of , i.e.

How are and to be determined? First, it should be noted that in the absence of diffusion, D = 0, then is

(See for example section 3.1, equation (22).) The simplest possibility for is

Using the definition of R, equation (54), and the results from equations (53) and (54), an elementary calculation gives

This formula provides a functional relationship between the space and time step-sizes. Such relations are generally expected for partial differential equations 6,9.

In summary, to calculate a numerical solution to equations (41a), (41b), and (41c), the following procedure should be carried out:

• (a) First, select the parameters (r, a, D).

• (b) Second, choose a value for Δ x and the parameter γ≥0, and determine Δ t from equation (55a).

• (c) Third, calculate, sequentially from, respectively, equations (47a), (52), and (47c).

4. Discussion

The results obtained in section 3 demonstrate the power of NSFD schemes to provide dynamical consistent (DC) finite difference discretizations of systems of nonlinear, coupled ODEs modeling interacting populations having a conservation law. In contrast to the physical and engineering sciences where a conservation law implies that some quantity is constant, in this work we modify the usual definition and extend it to the case where the total population is time-dependent, but goes to a limiting finite value as t→∞. Such systems and the corresponding NSFD discretizations resolve two issues. First, the discrete schemes obey the positivity requirement for the various subpopulations. Second, these schemes only produce bounded numerically computed solutions. In addition, no spurious fixed-points are introduced. Also, the example in section 3.3 shows that these techniques can easily be extended to the situation where diffusion also occurs in the system.

The critical aspect in the construction of NSFD schemes is determining which features of the original differential equations should be incorporated into the discrete model. Clearly, if some essential aspect of the equations is not also a feature of the discrete finite difference models, then numerical instabilities will appear.

Finally, with regard to systems in the biosciences, many satisfy a conservation law. The books by Thieme 4, Murray 13, and Diekmann and Heesterbeek 13 provide important examples. Two other relevant resources are the papers of Liu et al. 15, 16.

Acknowledgements

This research was supported in part by grants from DOE and MBRS-SCORE Program at Clark Atlanta University.

Additional information

Notes on contributors

Ronald E. Mickens

Email: [email protected]