A class of nonstandard numerical methods for autonomous dynamical systems

Accepted Manuscript A class of nonstandard numerical methods for autonomous dynamical systems Daniel T. Wood, Hristo V. Kojouharov PII: DOI: Reference:

S0893-9659(15)00191-3 http://dx.doi.org/10.1016/j.aml.2015.06.008 AML 4810

To appear in:

Applied Mathematics Letters

Received date: 15 April 2015 Revised date: 2 June 2015 Accepted date: 2 June 2015 Please cite this article as: D.T. Wood, H.V. Kojouharov, A class of nonstandard numerical methods for autonomous dynamical systems, Appl. Math. Lett. (2015), http://dx.doi.org/10.1016/j.aml.2015.06.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to view linked References

A class of nonstandard numerical methods for autonomous dynamical systems Daniel T. Wooda,∗, Hristo V. Kojouharova a Department

of Mathematics, The University of Texas at Arlington, P.O. Box 19408, Arlington, TX, 76019-0408, USA

Abstract A nonstandard finite difference method for solving autonomous dynamical systems is constructed. The proposed numerical method is computationally efficient and easy to implement. It is designed so that it preserves positivity of solutions and the local behavior of the dynamical system near equilibria. Keywords: finite difference, nonstandard, elementary stable, positivity 2010 MSC: 37N30, 65L05, 65L12, 65P40 1. Introduction Consider an n-dimensional dynamical system of the form dx = f (x); dt

(1)

x(0) ∈ Rn+ ,

whose solutions are of the form x = [x1 , x2 , . . . , xn ]T : R+ → Rn+ and f is a C 2 function that satisfies conditions which guarantee that Rn+ is positively invariant (e.g. see Appendix B of [1]). The set Rn+ is defined as Rn+ := {(x1 , x2 , . . . , xn ) ∈ Rn : x1 ≥ 0, x2 ≥ 0, . . . , xn ≥ 0},

and it is assumed that the system (1) has a finite number of equilibria, each of which is hyperbolic. There have been several successful attempts at designing numerical methods which preserve positivity of solutions and the local behavior, but only for special classes of systems of ordinary differential equations which exhibit certain structure (e.g. [2, 3, 4, 5, 6, 7, 8, 9]). Here, a finite difference method is constructed in Section 2 which preserves the aforementioned qualitative properties for the system (1). The method is analyzed in Section 3 and a numerical example is presented in Section 4 to show the utility of the method. In Section 5, some concluding remarks are given. ∗ Corresponding

author Email addresses: [email protected] (Daniel T. Wood), [email protected] (Hristo V. Kojouharov)

Preprint submitted to Elsevier

June 2, 2015

2. Construction of the Numerical Method A finite difference method to approximate the system (1) can be written as Dh (xk ) = Fh (f ; xk ), (2) dx , xk ≈ x(tk ), Fh (f ; xk ) approximates f (x) in the syswhere Dh (xk ) ≈ dt t=tk tem (1) and tk = t0 + kh, where h > 0. The positive and elementary stable properties [10] of a numerical method are defined below: Definition 1. The finite difference method (2) is called positive, if, for any value of the step size h, and x0 ∈ Rn+ its solution remains positive, i.e., xk ∈ Rn+ for all k ∈ N. Definition 2. The finite difference method (2) is called elementary stable if, for any value of the step size h, its fixed points x ¯ are the same as the equilibria of the differential system (1) and the local stability properties of each x ¯ are the same for both the differential system and the difference method. For a numerical method to accurately approximate the system (1), it must always exhibit a nonnegative solution when the initial value is nonnegative. The need for elementary stability is derived from the fact that any approximation to the system (1) should exhibit the same local stability around equilibria, and should not introduce any spurious fixed points. In order to design a numerical method which is both positive and elementary stable, the nonstandard finite difference (NSFD) method framework, first introduced by Mickens [11], is adopted. The following concise definition by Anguelov and Lubuma [10] is used: Definition 3. The one-step finite-difference scheme (2) for solving System (1) is a NSFD method if at least one of the following conditions is satisfied: • Dh (xk ) =

xk+1 − xk , where ϕ(h) = h + O(h2 ) is a non-negative function; ϕ(h)

• Fh (f ; xk ) = g(xk , xk+1 , h), where g(xk , xk+1 , h) is a non-local approximation of the right-hand side of System (1). It should be noted that a more robust and technical definition is given by Lubuma and Patidar in [12]. Using the above NSFD framework, the new numerical method which approximates the system (1) is constructed as follows: xik+1 − xik = f i (xk )ωki , φ(h)

2

(3)

where ωki

=

  

1, if f i (xk ) ≥ 0

xik+1 , xik

if f i (xk ) < 0

(4)

for i = 1, 2, . . . , n and φ(h) = h + O(h2 ) is a denominator function which is to be determined in the next section. The numerical treatment of the right-hand side function in the Method (3) is motivated by some of the discretization ideas of J. Benz et al. [13]. In [13], the authors propose a numerical method that works for both conservative and non-conservative systems by using a term similar to equation (4) which switches based on the conservativity-property of the system. Remark 1. There are a variety of functions that satisfy the condition φ(h) = h + O(h2 ). One of the most widely used denominator functions in the NSFD framework is φ(h) = 1 − e−h . However, stronger requirements are required here to ensure that the method is elementary stable. 3. Analysis of the Numerical Method It can be easily shown that the numerical method (3) is a first-order method [10]. The zero-stability of the method follows from the fact that ωk → 1 as h → 0, which implies that Method (3) is convergent. 3.1. Positivity The numerical method (3) can be written in an explicit form as follows:  i  xk + φ(h)f i (xk ), if f i (xk ) ≥ 0 i (xik )2 (5) xk+1 = , if f i (xk ) < 0  i xk − φ(h)f i (xk )

for i = 1, 2, . . . , n. To show that the method is positive, observe that when xik ≥ 0 and f i (xk ) ≥ 0, then xik+1 = xik + φ(h)f i (xk ) ≥ 0. Also, when xik ≥ 0 and f i (xk ) < 0, then (xik )2 (xik )2 ≥ 0 and −φ(h)f i (xk ) > 0, thus xik+1 = i ≥ 0. xk − φ(h)f i (xk ) 3.2. Elementary Stability In order to ensure that the numerical method (3) is elementary stable, a special denominator function, φ(h), must be selected. Let ϕ(h) be a real-valued function which satisfies 0 < ϕ(h) < 1 for all h > 0 and ϕ(h) = h + O(h2 ). A relationship between the eigenvalues of the Jacobians of the differential system (1) and the numerical method (3) at the equilibrium points can be used to guarantee that both (1) and (3) share the same local stability.

3

Denote by A the Jacobian of the differential system (1) evaluated at an equilibrium, x ¯. Then the Jacobian, B, of the numerical method (3) evaluated at x ¯ can be written as B = I + ϕ(h)A. The following relationship holds:    n B−I 1 det (A − λI) = det − λI = det (B − [1 + ϕ(h)λ]I) . ϕ(h) ϕ(h)

(6)

Therefore λ is an eigenvalue of A if and only if 1 + ϕ(h)λ is an eigenvalue of B. Now, take   |λ|2 Q ≥ max (7) 2|Reλ|

across all equilibria, x ¯. Then, by choosing the denominator function φ(h) =

ϕ(qh) , q

(8)

where q > Q, the numerical method (3) can be shown to be elementary stable as follows: Equations (7) and (8), along with the fact that 0 < ϕ(Qh) < 1, imply that for any λ, 2|Reλ| φ(h) < |λ|2

which yields

−2φ(h)|Reλ| + φ(h)2 |λ|2 < 0.

If the system (1) is locally stable at x ¯, then Re(λ) < 0. Using (6), which gives a relation between the eigenvalues of differential system (1) and the numerical method (3), the eigenvalues, 1 + φ(h)λ, of B satisfy p |1 + φ(h)λ| = 1 − 2φ(h)|Reλ| + φ(h)2 |λ|2 < 1. Next, if the system (1) is unstable at x ¯, then there exists an eigenvalue λ0 such that Re(λ0 ) > 0. This implies that p |1 + φ(h)λ0 | = 1 + 2φ(h)Reλ0 + φ(h)2 |λ0 |2 > 1.

This result demonstrates that the numerical method (3) shares the same local stability at each equilibrium as the original system (1) and is, thus, elementary stable. Remark 2. The inequality (7) used to find Q has been previously employed in other articles (e.g. most recently in [2]). Also, instead of finding the eigenvalues directly, one can make use of the Routh-Hurwitz conditions to more easily find a suitable constant Q (see [9]). Remark 3. Method (3) can easily be modified to work for dynamical systems with negative solutions, by switching the inequalities in equation (4). 4

4. Numerical Simulations For illustrative purposes, the numerical method (3) is applied to the following two dimensional system: dx = y − x; dt (9) dy = sin(xy) − y; dt with initial conditions x(0) = 2, y(0) = 0.5. The system (9) has one trivial equilibrium at x = 0, y = 0, which is locally asymptotically stable. Here, we have selected the denominator function for the NSFD method (3) as φ(h) = (1 − exp(−qh))/q, where q = 0.6 satisfies q > Q = 0.5. The value Q is calculated following the procedure described in Subsection 3.2. In Figure 1, the first order NSFD method (3) is compared with the Heun’s method. One can see that they both converge for the small time step h = 0.01. However, the Heun’s method does not preserve positivity and does not approach the equilibrium for the larger time step h = 2, while the NSFD method (3) still preserves all the desired properties. Remark 4. Since the system (9) cannot be written as a productive-destructive system, the other NSFD methods available in the literature (e.g. [2, 9]) cannot be applied to preserve both the positivity of solutions and the local behavior of the dynamical system near equilibria.

0.6

NSFD Heun’s Equilibrium

NSFD Heun’s Equilibrium

2

1

y

y

0.4 0 0.2 −1 0 0

0.5

1

x

1.5

2

2.5

0

2

4

x

6

8

10

12

Figure 1: Numerical approximations of the system (9) using the NSFD method and the Heun’s method with time-step h = 0.01 (left) and h = 2 (right).

5. Conclusions A first-order nonstandard finite difference method is constructed for solving autonomous dynamical systems with positive solutions. It preserves two important qualitative dynamical properties, namely, positivity of solutions and elementary stability. The method is both computationally efficient and easy to implement, and it can be used to solve a broad range of problems in science and engineering. 5

References [1] H. Smith, P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, Cambridge, 1995. [2] R. Anguelov, Y. Dumont, J. Lubuma, M. Shillor, Dynamically consistent nonstandard finite difference schemes for epidemiological models, Journal of Computational and Applied Mathematics 255 (2014) 161 – 182. [3] D. Dimitrov, H. Kojouharov, Positive and elementary stable nonstandard numerical methods with applications to predator-prey models, Journal of Computational and Applied Mathematics 189 (1–2) (2006) 98–108. [4] D. Dimitrov, H. Kojouharov, Dynamically consistent numerical methods for general productive-destructive systems, Journal of Difference Equations and Applications 17 (12) (2011) 1721–1736. [5] H. Obaid, R. Ouifki, K. Patidar, An unconditionally stable nonstandard finite difference method applied to a mathematical model of HIV infection, International Journal of Applied Mathematics and Computer Science 23 (2) (2013) 357–372. [6] L. Roeger, Dynamically consistent discrete-time Lotka-Volterra competition models, Discrete and Continuous Dynamical Systems (2009) 650–658. [7] L. Roeger, General nonstandard finite-difference schemes for differential equations with three fixed-points, Computers & Mathematics with Applications 57 (3) (2009) 379 – 383. [8] A. Suryanto, W. Kusumawinahyu, I. Darti, I. Yanti, Dynamically consistent discrete epidemic model with modified saturated incidence rate, Computational and Applied Mathematics 32 (2) (2013) 373–383. [9] D. Wood, D. Dimitrov, H. Kojouharov, A nonstandard finite difference method for n-dimensional productive-destructive systems, Journal of Difference Equations and Applications 21 (3) (2015) 240–254. [10] R. Anguelov, J. Lubuma, Contributions to the mathematics of the nonstandard finite difference method and applications, Numerical Methods for Partial Differential Equations 17 (5) (2001) 518–543. [11] R. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, Singapore, 1994. [12] J.-S. Lubuma, K. Patidar, Contributions to the theory of non-standard finite-difference methods and applications to singular perturbation problems, World Scientific, Sinagpore, 2005. [13] J. Benz, A. Meister, P. Zardo, A conservative, positivity preserving scheme for advection-diffusion-reaction equations in biochemical applications, in: Proceedings of Symposia in Applied Mathematics, American Mathematical Society, 2009, p. 399. 6