Research Article  Open Access
Rafael J. Villanueva, Abraham J. Arenas, Gilberto GonzálezParra, "A Nonstandard Dynamically Consistent Numerical Scheme Applied to Obesity Dynamics", Journal of Applied Mathematics, vol. 2008, Article ID 640154, 14 pages, 2008. https://doi.org/10.1155/2008/640154
A Nonstandard Dynamically Consistent Numerical Scheme Applied to Obesity Dynamics
Abstract
The obesity epidemic is considered a health concern of paramount importance in modern society. In this work, a nonstandard finite difference scheme has been developed with the aim to solve numerically a mathematical model for obesity population dynamics. This interacting population model represented as a system of coupled nonlinear ordinary differential equations is used to analyze, understand, and predict the dynamics of obesity populations. The construction of the proposed discrete scheme is developed such that it is dynamically consistent with the original differential equations model. Since the total population in this mathematical model is assumed constant, the proposed scheme has been constructed to satisfy the associated conservation law and positivity condition. Numerical comparisons between the competitive nonstandard scheme developed here and Euler's method show the effectiveness of the proposed nonstandard numerical scheme. Numerical examples show that the nonstandard difference scheme methodology is a good option to solve numerically different mathematical models where essential properties of the populations need to be satisfied in order to simulate the real world.
1. Introduction
Systems of nonlinear ordinary differential equations are used to model different kind of diseases in mathematical epidemiology (see [1]). In these models, the variables commonly represent populations of susceptible, infected, recovered, transmitted disease vectors, and so forth. Since the system describes the evolution of different classes of populations, a solution over the time is required. The ideal scenario is to obtain analytical solutions, but in most of the cases, this is not possible to achieve. Therefore, it is necessary to turn to numerical methods to obtain numerical simulations.
Traditional schemes, like forward Euler, RungeKutta, and other numerical methods used to solve nonlinear initial value problems, sometimes fail generating oscillations, bifurcations, chaos, and false steady states [2]. Methods that use adaptive step size, also may fail [3]. One approach to avoid this class of numerical instabilities is the construction of nonstandard finite difference schemes. This technique, developed by Mickens [4, 5], has brought a lot of applications where the nonstandard methods have been applied to various problems in science and engineering [6–14] in which the numerical solutions preserve properties of the solution of the continuous model, and additionally, in some cases, it is possible to use large time step sizes. Furthermore, it has been shown that the nonstandard finite difference schemes, generated using the rules created by Mickens [4], generally provide accurate numerical solutions to differential equations.
The aim of this paper is to develop a nonstandard finite difference (NSFD) scheme free of numerical instabilities in order to obtain the numerical solution of a mathematical model of infant obesity with constant population size presented in [15]. This interacting population model represented as a system of coupled nonlinear ordinary differential equations exhibits two steady states: a trivial steady state called obesity free equilibrium (OFE) and a nontrivial steady state called obesity endemic equilibrium (OEE). The general philosophy for constructing NSFD schemes is to obtain numerical solutions dynamically consistent with the underlying continuous model. This means that all of the critical, qualitative properties of the solutions to the system of differential equations should also be satisfied by the solutions of the discrete scheme [16].
The design of a nonstandard scheme is not a straightforward task, in fact, many schemes may be developed for one model and several can fail to converge. Therefore the innovative part in this work is the construction of a NSFD scheme, such that is consistent with the properties of continuous dynamic model, as positivity and conservation of the population. All the numerical simulations with different parameter values suggest that the NSFD scheme developed here preserves numerical stability in larger regions in comparison to the smaller regions of other standard numerical methods. However, theoretical justification of this fact is not possible due to unmanageable analytic expression of the eigenvalues of the Jacobian matrix corresponding to the linearization at the equilibrium points. It is important to remark that most of the previous applications of nonstandard methods to ODE's have been done to one, two, and three equations, where theoretical analysis is easier.
Since subpopulations must never take on negative values, the proposed NSFD scheme is designed to satisfy positivity condition. Furthermore, due to the fact that in the mathematical model the population is assumed constant, the NSFD scheme has been developed in order to always exactly satisfy the associated conservation law [17].
The organization of this paper is as follows. In Section 2, the mathematical model for the evolution of overweight and obesity population is presented. The construction of the proposed NSFD numerical scheme is carried out in Section 3. In Section 4, a basic analysis of the stability of the mathematical model is developed using data of the region of Valencia, Spain. Additionally, numerical analysis of the NSFD numerical scheme is done. Numerical simulations for the NSFD scheme are reported also in this section in order to show its computational advantages. Discussion and conclusions are presented in Section 5.
2. Mathematical Model
In [15], a dynamic obesity model for the evolution of infant overweight and obesity population was introduced. This continuous model is based on the partition of the infant population into six subpopulations. In this model, denotes the proportion of normalweight children, denotes the proportion of latent children, that is, children with the habit of high consumption of unhealthy food but who are still normal weight, denotes the proportion of children with overweight, denotes the proportion of obese children, denotes the proportion of overweight children on diet, and denotes the proportion of obese children on diet. The model is represented by the following nonlinear system of ordinary differential equations: where the timeinvariant parameters of the model are (1): social transmission rate to high consumption of unhealthy food due to social environment pressure,(2): inflow rate population in the system or inversely proportional to average time spent by an individual in the system,(3): rate at which a latent individual moves to the overweight subpopulation,(4): rate at which an overweight individual becomes an obese individual by continuous consumption of unhealthy food,(5): rate at which an overweight individual on diet becomes a normalweight individual,(6): rate at which an overweight individual stops or reduces unhealthy food consumption, that is, the individual is on diet,(7): rate at which an overweight individual on diet fails, that is, the individual resumes a high unhealthy food consumption,(8): rate at which an obese individual stops or reduces unhealthy food consumption,(9): rate at which an obese individual on diet fails,(10): rate at which an obese individual on diet becomes an overweight individual on diet. In this model, since the constant population is normalized to unity one gets for all time ,We define (2.2) as the conservation law associated with the system (2.1) and this equation must be satisfied by the NSFD numerical scheme for any time .
3. Nonstandard Finite Difference Discretization
In order to avoid dynamic inconsistencies such as oscillations and numerical instabilities in this section, an NSFD scheme for solving numerically system (2.1) is constructed. A numerical scheme is called NSFD if at least one of the following conditions is satisfied [8, 18]:(1) nonlocal approximation is used [4, 5, 7];(2) discretization of derivative is not traditional and a nonnegative functioncalled a denominator function is used [19]. The discretization is based on the approximations of the temporal derivatives by a firstorder generalized forward scheme. If , we define an equivalent derivative aswhere is a nonnegative realvalued function defined in (3.1). Note that the above definition is consistent with the traditional definition of the derivative. IndeedAn example of a function that satisfies the above conditions is (see [4, 19]).
3.1. Construction of the NSFD Scheme
In this subsection, the main goal is to construct a dynamically consistent numerical discrete scheme for solving system (2.1). Notice that all the model variables (subpopulations) and parameters are positive. Therefore, in order to obtain a dynamically consistent discrete scheme, we must ensure that the resulting discrete solutions are all positive which is necessary to avoid schemedependent instabilities. Furthermore, since the population is assumed constant, the NSFD scheme should satisfy the associated conservation law [17]. These aspects will be taken into account in the construction of the numerical scheme below.
Let us denote by , , , , , and the approximations of , , , , , and , respectively, for and the time step size of the scheme. Using (3.2), the firstorder numerical scheme to obtain solutions , , , , , and of the model (2.1) is given by the following implicit scheme: and after rearranging, we obtain the following explicit form: where the denominator function is chosen aswith in order to guarantee the positivity condition. However, positivity can also be guaranteed taking the traditional if . Therefore throughout this paper when , we use .
Remark 3.1. Notice that the populations must be calculated in a particular order. This sequential form of calculation is a generic feature of NSFD methods [4, 16, 19].
Remark 3.2. The incorporation of (2.2) to system (2.1) introduces more complexity in order to satisfy that, if the same term occurs in more than one differential equation, then it must be modeled discretely the same way in all equations [4, 17].
3.2. Properties and Computation in the NSFD Scheme
It can be seen from system (3.5)–(3.10) that we have constructed an NSFD scheme for the model (2.1) having the following properties:(1) positivity. If , , , , , , then , , , , , , for all since that ;(2) satisfy conservation law, that is, constant population. From system (3.4) one gets that It follows by induction that ifNow, observe that the subpopulations must be calculated from scheme (3.5)–(3.10) in the following order:(1) select values , , , , , , such that ;(2) determine from (3.5)) using and ;(3) calculate from (3.6) with , , and ;(4) obtain from (3.7) using , , , , and ;(5) find from (3.8) with , , , and ;(6) get from (3.9) using , , and ;(7) compute from (3.10) with the values of , and . Note that the sequential form of calculation is a generic feature of NSFD schemes [17]. In addition, it can be seen that the main part of the local truncation error associated with each equation of system (3.2) is of order , confirming that the constructed NSFD scheme is firstorder accurate.
4. Numerical Results and Dynamic Consistency
In this section, some theoretical properties of the mathematical model of overweight and obesity population dynamics proposed in [15] for a particular set of values of the parameters are studied, in order to know a priori the correct behavior that needs to come out from the numerical solution corresponding to the NSFD numerical scheme (3.4). Furthermore, this section is devoted to show numerically the dynamic consistency and numerical advantages of the developed NSFD scheme and some comparisons with other wellknown numerical methods. The chosen values for the set of parameters of model (2.1) are the ones used in [15], and are shown in Table 1 in a week scale.

4.1. Numerical Stability of the Mathematical Model
A linear stability analysis of system (2.1) for a particular set of values of the parameters is developed here in order to check numerically dynamic consistency between the continuous model and the NSFD discrete scheme. A more general stability analysis has not been performed since general parameters gives unmanageable analytic expressions. Nevertheless, several sets of values of the parameters are used to suggest numerically the aforementioned dynamic consistency.
It is well known that an equilibrium point is an asymptotically stable node if and only if the real part of the eigenvalues of the Jacobian associated system evaluated at the equilibrium points are negative.
System (2.1) has two steady states; OFE and OEE. For the particular set of values of the parameters presented in Table 1 for system (2.1), the following equilibrium points are obtained:
The eigenvalues of the Jacobian evaluated at the OFE point and at the positive OEE point are shown in Table 2. The OFE point is unstable, since the eigenvalue is positive. On the other hand, the OEE point is locally asymptotically stable, due to the fact that all real parts of the eigenvalues of the Jacobian of system (2.1) evaluated at the OEE point are negative. Therefore, a dynamic consistent numerical scheme should reproduce this numerical behavior of the continuous model. In Section 4.2, this fact will be checked.

4.2. Numerical Simulations
This subsection is devoted to show numerically the dynamic consistency and numerical advantages of the developed NSFD scheme. One basic property that should satisfy the NSFD numerical scheme (3.5)–(3.10) in order to be dynamically consistent is that their fixed points should be the same equilibrium points as the continuous model (2.1). The equilibrium points of the NSFD numerical scheme (3.4) are obtained by setting to zero the lefthand sides of system (3.4). It is easy to prove that the NSFD numerical scheme (3.4) has the same two steady states; that the continuous model for any set of parameter values.
Numerical simulation for different sets of values of the parameters was performed to investigate the dynamic consistency of the developed NSFD numerical scheme (3.4). In Figure 1, it can be seen that the NSFD scheme (3.4) converges to the same OEE point of the continuous model using the particular set of values of the parameters shown in Table 1.
(a)
(b)
(c)
(d)
(e)
(f)
Additionally in Figure 2, it is observed that despite the initial condition is close to the OFE point, the numerical solution moves further away from this equilibrium point, suggesting numerically the instability of the OFE point and the consistency of the NSFD numerical scheme with the continuous model (2.1).
(a)
(b)
(c)
(d)
(e)
(f)
4.3. Numerical Comparisons of the NSFD Scheme
In order to study the effect of the time step size in the NSFD numerical scheme and dynamic consistency related to local stability, the spectral radius of the Jacobian corresponding to the linearization of NSFD scheme and Euler schemes are computed for different time steps for system (2.1). Using general parameters values, one gets a general Jacobian. However, the functional relationship between the spectral radius and the time step size gives an unmanageable complex analytic expression. Therefore, in order to investigate numerically the dynamic consistency of the NSFD scheme, the particular set of values of the parameters presented in Table 1 (but in year scale) is used for numerical results. However, other sets of values of the parameters have been used to confirm the previous numerical results.
It is well known that a necessary and sufficient condition for the convergence of a fixed point scheme is that the spectral radius of the Jacobian satisfies . In Table 3, the spectral radius of the Jacobian at the OEE point of the NSFD and Euler schemes for different time steps are presented. Table 3 compares the convergence properties of the Euler scheme with that of NSFD scheme when used to integrate system (2.1) subject to the same initial conditions and parameter values. It is clear from Table 3 that the NSFD scheme is more competitive in terms of numerical stability. In all of these simulations, the NSFD scheme is seen to be monotonically convergent to the correct endemic equilibrium point. However, Euler fails to converge for a time step size of . The RungeKutta fourthorder scheme improves the result obtained by Euler scheme, but fails to converge for time step sizes when the particular set of values of the parameters presented in Table 1 is used. As expected in other numerical simulations with several different values of the parameters, the RungeKutta fourthorder method outperform Euler scheme, since convergence is achieved for larger time step sizes.

Table 3 compares the convergence properties of the Euler scheme with that of the NSFD scheme for a specific set of values of the parameters obtained from a realworld application. However, several different values of the parameters have been used and numerical results suggest that the NSFD scheme developed here preserves numerical stability in larger regions in comparison to the smaller regions of the Euler numerical scheme. However, theoretical justification of this fact is not possible due to unmanageable analytic expression of the eigenvalues of the Jacobian matrix corresponding to the linearization at the equilibrium points. It is important to remark that most of the previous applications of nonstandard methods to ODE's have been done to one, two, and three equations, where theoretical analysis is easier. In Figure 3, it can be seen that the NSFD scheme converges and the Euler scheme oscillates incorrectly but converges to the OEE point for a time step size . In Figure 4, it is depicted that the NSFD scheme converges to the correct OEE point for a time step size , but in this case, Euler fails to converge. Finally, Figure 5 shows that the NSFD scheme converges to the OEE point taking only positive values. Nevertheless, it can be observed that the several routines from MATLAB package produce incorrect negative values.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
5. Discussion and Conclusions
In this paper, we applied the NSFD methodology to develop a scheme to solve numerically a mathematical model for obesity population dynamics. This interacting population model represented as a system of coupled nonlinear ordinary differential equations can be used to analyze, understand, and predict the dynamics of overweight and obese populations.
The advantages of the NSFD scheme developed here are that they preserve numerical stability in larger regions for the time step size in comparison to the smaller regions of numerical stability of the approximations obtained by the other standard numerical methods. Furthermore, the proposed scheme satisfies positivity condition and the conservation law that any good scheme for mathematical models of population dynamics must fulfill. The construction of these nonstandard schemes is not a straightforward task, in fact, many schemes may be developed for one model and several can fail to converge. Nevertheless, the principle of dynamic consistency can be used with great efficiency to place restrictions on the design of NSFD schemes in order to obtain efficient schemes, as has been shown in this work.
The NSFD scheme proposed here was analyzed and tested in several numerical simulations. All the numerical simulations performed in this work show that the NSFD scheme converges to the OEE point for any time step size and satisfies a condition of positivity required for populations. Finally, it is important to remark that the developed NSFD scheme is easy to use and convergence is achieved for larger time step sizes than the Euler or the RungeKutta fourthorder schemes. The NSFD scheme also outperforms MATLAB package routines since negative values for populations are obtained with these routines. Therefore, for realworld applications in nature and society, where numerical values are required it is important to be cautious about which numerical scheme is more suitable to solve the mathematical model since unreal results can be obtained.
References
 F. Brauer and C. CastilloChávez, Mathematical Models in Population Biology and Epidemiology, vol. 40 of Texts in Applied Mathematics, Springer, New York, NY, USA, 2001. View at: Zentralblatt MATH  MathSciNet
 J. D. Lambert, Numerical Methods for Ordinary Differential Systems: The Initial Value Problem, John Wiley & Sons, Chichester, UK, 1991. View at: Zentralblatt MATH  MathSciNet
 S. M. Moghadas, M. E. Alexander, B. D. Corbett, and A. B. Gumel, “A positivitypreserving Mickenstype discretization of an epidemic model,” Journal of Difference Equations and Applications, vol. 9, no. 11, pp. 1037–1051, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, River Edge, NJ, USA, 1994. View at: Zentralblatt MATH  MathSciNet
 R. E. Mickens, Application of Nonstandard Finite Difference Schemes, World Scientific, River Edge, NJ, USA, 2000.
 S. M. Moghadas, M. E. Alexander, and B. D. Corbett, “A nonstandard numerical scheme for a generalized Gausetype predatorprey model,” Physica D, vol. 188, no. 12, pp. 134–151, 2004. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Anguelov and J. M.S. Lubuma, “Nonstandard finite difference method by nonlocal approximation,” Mathematics and Computers in Simulation, vol. 61, no. 3–6, pp. 465–475, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Anguelov and J. M.S. Lubuma, “Contributions to the mathematics of the nonstandard finite difference method and applications,” Numerical Methods for Partial Differential Equations, vol. 17, no. 5, pp. 518–543, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 W. Piyawong, E. H. Twizell, and A. B. Gumel, “An unconditionally convergent finitedifference scheme for the SIR model,” Applied Mathematics and Computation, vol. 146, no. 23, pp. 611–625, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. B. Gumel, “A competitive numerical method for a chemotherapy model of two HIV subtypes,” Applied Mathematics and Computation, vol. 131, no. 23, pp. 329–337, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 D. T. Dimitrov and H. V. Kojouharov, “Nonstandard finitedifference schemes for general twodimensional autonomous dynamical systems,” Applied Mathematics Letters, vol. 18, no. 7, pp. 769–774, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. Jansen and E. H. Twizell, “An unconditionally convergent discretization of the SEIR model,” Mathematics and Computers in Simulation, vol. 58, no. 2, pp. 147–158, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 K. C. Patidar, “On the use of nonstandard finite difference methods,” Journal of Difference Equations and Applications, vol. 11, no. 8, pp. 735–758, 2005. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 B. M. ChenCharpentier, D. T. Dimitrov, and H. V. Kojouharov, “Combined nonstandard numerical methods for ODEs with polynomial righthand sides,” Mathematics and Computers in Simulation, vol. 73, no. 1–4, pp. 105–113, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 L. Jódar, F. J. Santonja, and G. GonzálezParra, “Modeling dynamics of infant obesity in the region of Valencia, Spain,” Computers & Mathematics with Applications, vol. 56, no. 3, pp. 679–689, 2008. View at: Publisher Site  Google Scholar  MathSciNet
 R. E. Mickens, “Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations,” Journal of Difference Equations and Applications, vol. 11, no. 7, pp. 645–653, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. E. Mickens, “Numerical integration of population models satisfying conservation laws: NSFD methods,” Journal of Biological Dynamics, vol. 1, no. 4, pp. 427–436, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 J. M.S. Lubuma and K. C. Patidar, “Contributions to the theory of nonstandard finite difference methods and applications to singular perturbation problems,” in Advances in the Applications of Nonstandard Finite Difference Schemes, R. E. Mickens, Ed., pp. 513–560, World Scientific, Hackensack, NJ, USA, 2005. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 R. E. Mickens, “Calculation of denominator functions for nonstandard finite difference schemes for differential equations satisfying a positivity condition,” Numerical Methods for Partial Differential Equations, vol. 23, no. 3, pp. 672–691, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
Copyright
Copyright © 2008 Rafael J. Villanueva et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.