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, Runge-Kutta, 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 normal-weight 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 time-invariant 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 normal-weight 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 first-order generalized forward scheme. If , we define an equivalent derivative aswhere is a nonnegative real-valued 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 scheme-dependent 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 first-order 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 first-order 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 well-known 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.

Table 1: Parameters (week scale) of the model given for the nonlinear ordinary differential equation system (3.4) for the region of Valencia, Spain.

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.

Table 2: The eigenvalues of the Jacobian of system (2.1) evaluated at the OFE point and at the OEE point.

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 left-hand 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.

Figure 1: Profiles of subpopulations and . Solutions are obtained by the proposed NSFD scheme with , where and initial conditions are , , , , , and .

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).

Figure 2: Profiles of subpopulations , and obtained with the NSFD scheme with , and initial conditions are and .

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 Runge-Kutta fourth-order 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 Runge-Kutta fourth-order method outperform Euler scheme, since convergence is achieved for larger time step sizes.

Table 3: Spectral radius for different time step size of the Euler and NSFD numerical schemes.

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 real-world 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.

Figure 3: Profiles of subpopulations , and . Solutions are obtained by Euler (dashed) and proposed NSFD scheme (line) with and for system (2.1).

Figure 4: Profiles of subpopulations , and . Solutions are obtained by Euler scheme (dashed) and developed NSFD scheme (line) with and for system (2.1).

Figure 5: Solutions for subpopulation obtained by NSFD scheme and routines ODE45, ODE23s, ODE113, and ODE23 from MATLAB package (dashed), using .

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 Runge-Kutta fourth-order schemes. The NSFD scheme also outperforms MATLAB package routines since negative values for populations are obtained with these routines. Therefore, for real-world 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.