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, 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:𝑁^′(𝑡)=𝜇+𝜀𝐷_𝑆,𝐿(𝑡)−𝜇𝑁(𝑡)−𝛽𝑁(𝑡)𝐿(𝑡)+𝑆(𝑡)+𝑂(𝑡)_′(𝑡)=𝛽𝑁(𝑡)𝐿(𝑡)+𝑆(𝑡)+𝑂(𝑡)−𝜇+𝛾_𝐿𝑆𝐿(𝑡),_′(𝑡)=𝛾_𝐿𝐿(𝑡)+𝜑𝐷_𝑆(𝑡)−𝜇+𝛾_𝑆𝑂+𝛼𝑆(𝑡),_′(𝑡)=𝛾_𝑆𝑆(𝑡)+𝛿𝐷_𝑂𝐷(𝑡)−[𝜇+𝜎]𝑂(𝑡),_′𝑆(𝑡)=𝛾_𝐷𝐷_𝑂(𝑡)+𝛼𝑆(𝑡)−[𝜇+𝜀+𝜑]𝐷_𝑆𝐷(𝑡),_′𝑂(𝑡)=𝜎𝑂(𝑡)−𝜇+𝛾_𝐷+𝛿𝐷_𝑂(𝑡), (2.1) 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 𝑡,𝑁(𝑡)+𝐿(𝑡)+𝑆(𝑡)+𝑂(𝑡)+𝐷_𝑆(𝑡)+𝐷_𝑂(𝑡)=1. (2.2)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 function𝜓(ℎ)=ℎ+𝒪ℎ² (3.1)called 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 as𝑑𝑓(𝑡)=𝑑𝑡𝑓(𝑡+ℎ)−𝑓(𝑡)𝜓(ℎ)+𝒪𝜓(ℎ)asℎ⟶0, (3.2)where 𝜓(ℎ) is a nonnegative real-valued function defined in (3.1). Note that the above definition is consistent with the traditional definition of the derivative. Indeed𝑑𝑓(𝑡)𝑑𝑡=lim_ℎ→0𝑓(𝑡+ℎ)−𝑓(𝑡)𝜓(ℎ)+𝒪𝜑(ℎ)=lim_ℎ→0𝑓(𝑡+ℎ)−𝑓(𝑡)ℎlim_ℎ→0ℎ𝜓(ℎ)+lim_ℎ→0𝒪=̇𝜓(ℎ)𝑓(𝑡). (3.3)An example of a function 𝜓(ℎ) that satisfies the above conditions is 𝜓(ℎ)=1−𝑒^−ℎ (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 𝑛=0,1,2,…, 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:𝐷_𝑂^𝑛+1−𝐷^𝑛_𝑂𝜓(ℎ)=𝜎𝑂^𝑛−𝛿+𝜇+𝛾_𝐷𝐷_𝑂^𝑛+1,𝐷_𝑆^𝑛+1−𝐷^𝑛_𝑆𝜓(ℎ)=−(𝜇+𝜖+𝜑)𝐷_𝑆^𝑛+1+𝛼𝑆^𝑛+𝛾_𝐷𝐷_𝑂^𝑛+1,𝑁^𝑛+1−𝑁^𝑛𝜓(ℎ)=𝜇−𝛽𝑁^𝑛+1𝑆^𝑛+𝑂^𝑛+𝐿^𝑛+𝜖𝐷_𝑆^𝑛+1−𝜇𝑁^𝑛+1,𝐿^𝑛+1−𝐿^𝑛𝜓(ℎ)=𝛽𝑁^𝑛+1𝑆^𝑛+𝑂^𝑛+𝐿^𝑛−𝛾_𝐿+𝜇𝐿^𝑛+1,𝑆^𝑛+1−𝑆^𝑛𝜓(ℎ)=𝛾_𝐿𝐿^𝑛+1+𝜑𝐷_𝑆^𝑛+1−𝛼𝑆^𝑛−𝛾_𝑆+𝜇𝑆^𝑛+1,𝑂^𝑛+1−𝑂^𝑛𝜓(ℎ)=𝛾_𝑆𝑆^𝑛+1−𝜎𝑂^𝑛−𝜇𝑂^𝑛+1+𝛿𝐷_𝑂^𝑛+1, (3.4) and after rearranging, we obtain the following explicit form:𝐷_𝑂^𝑛+1=𝐷^𝑛_𝑂+𝜓(ℎ)𝜎𝑂^𝑛1+𝜓(ℎ)𝛿+𝜇+𝛾_𝐷𝐷_𝑆𝑛+1,=𝐷^𝑛_𝑆+𝜓(ℎ)𝛼𝑆^𝑛+𝛾_𝐷𝐷_𝑂^𝑛+1𝑁1+𝜓(ℎ)(𝜖+𝜇+𝜑)_𝑛+1,=𝑁^𝑛+𝜓(ℎ)𝜇+𝜖𝐷_𝑆^𝑛+11+𝜓(ℎ)𝜇+𝛽𝑆^𝑛+𝑂^𝑛+𝐿^𝑛𝐿_𝑛+1,=𝐿^𝑛+𝜓(ℎ)𝛽𝑁^𝑛+1𝑆^𝑛+𝑂^𝑛+𝐿^𝑛1+𝜓(ℎ)𝛾_𝐿𝑆+𝜇_𝑛+1,=𝑆^𝑛1−𝜓(ℎ)𝛼+𝜓(ℎ)𝛾_𝐿𝐿^𝑛+1+𝜑𝐷_𝑆^𝑛+11+𝜓(ℎ)𝛾_𝑆𝑂+𝜇_𝑛+1,=𝑂^𝑛1−𝜓(ℎ)𝜎+𝜓(ℎ)𝛾_𝑆𝑆^𝑛+1+𝛿𝐷_𝑂^𝑛+11+𝜓(ℎ)𝜇, (3.5)(3.6)(3.7)(3.8)(3.9)(3.10) where the denominator function is chosen as𝜓(ℎ)=1−𝑒^−ℎ𝑀𝑀, (3.11)with 𝑀≥max{𝛼,𝜎} in order to guarantee the positivity condition. However, positivity can also be guaranteed taking the traditional 𝜓(ℎ)=ℎ if ℎ𝑀≤1. Therefore throughout this paper when ℎ𝑀≤1, 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 𝐷^𝑛_𝑂>0, 𝐷^𝑛_𝑆>0, 𝑁^𝑛>0, 𝐿^𝑛>0, 𝑆^𝑛>0, 𝑂^𝑛>0, then 𝐷_𝑂^𝑛+1>0, 𝐷_𝑆^𝑛+1>0, 𝑁^𝑛+1>0, 𝐿^𝑛+1>0, 𝑆^𝑛+1>0, 𝑂^𝑛+1>0, for all 𝑛=0,1,2,…, since that 1−𝜓(ℎ)𝑀>0;(2) satisfy conservation law, that is, constant population. From system (3.4) one gets that𝐷_𝑂^𝑛+1+𝐷_𝑆^𝑛+1+𝑁^𝑛+1+𝐿^𝑛+1+𝑆^𝑛+1+𝑂^𝑛+1=𝐷^𝑛_𝑂+𝐷^𝑛_𝑆+𝑁^𝑛+𝐿^𝑛+𝑆^𝑛+𝑂^𝑛+𝜇𝜓(ℎ).1+𝜇𝜓(ℎ) (3.12) It follows by induction that if𝐷^𝑛_𝑂+𝐷^𝑛_𝑆+𝑁^𝑛+𝐿^𝑛+𝑆^𝑛+𝑂^𝑛𝐷=1,then_𝑂𝑛+1+𝐷_𝑆^𝑛+1+𝑁^𝑛+1+𝐿^𝑛+1+𝑆^𝑛+1+𝑂^𝑛+1=1∀𝑛=0,1,…. (3.13)Now, observe that the subpopulations must be calculated from scheme (3.5)–(3.10) in the following order:(1) select values 𝐷^𝑛_𝑂>0, 𝐷^𝑛_𝑆>0, 𝑁^𝑛>0, 𝐿^𝑛>0, 𝑆^𝑛>0, 𝑂^𝑛>0, such that 𝐷^𝑛_𝑂+𝐷^𝑛_𝑆+𝑁^𝑛+𝐿^𝑛+𝑆^𝑛+𝑂^𝑛=1;(2) determine 𝐷_𝑂^𝑛+1 from (3.5)) using 𝐷^𝑛_𝑂 and 𝑂^𝑛;(3) calculate 𝐷_𝑆^𝑛+1 from (3.6) with 𝐷_𝑂^𝑛+1, 𝐷^𝑛_𝑆, and 𝑂^𝑛;(4) obtain 𝑁^𝑛+1 from (3.7) using 𝐷_𝑆^𝑛+1, 𝑁^𝑛, 𝑆^𝑛, 𝑂^𝑛, and 𝐿^𝑛;(5) find 𝐿^𝑛+1 from (3.8) with 𝑁^𝑛+1, 𝑆^𝑛, 𝑂^𝑛, and 𝐿^𝑛;(6) get 𝑆^𝑛+1 from (3.9) using 𝑆^𝑛, 𝐿^𝑛+1, and 𝐷_𝑆^𝑛+1;(7) compute 𝑂^𝑛+1 from (3.10) with the values of 𝑂^𝑛, 𝑆^𝑛+1 and 𝐷_𝑂^𝑛+1. 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:OFE=(1,0,0,0,0,0),OEE=(0.2937,0.2971,0.2708,0.1267,0.0080,0.0033). (4.1)

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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Profiles of subpopulations 𝑁(𝑡),𝐿(𝑡),𝑆(𝑡),𝑂(𝑡),𝐷𝑆(𝑡), and 𝐷𝑂(𝑡). Solutions are obtained by the proposed NSFD scheme with 𝜓(ℎ)=ℎ, where ℎ=0.2 and initial conditions are 𝑁(0)=0.462, 𝐿(0)=0.194, 𝑆(0)=0.2176, 𝑂(0)=0.09, 𝐷𝑆(0)=0.0249, and 𝐷𝑂(0)=0.0115.

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)

Figure 2 

Profiles of subpopulations 𝑁(𝑡),𝐿(𝑡),𝑆(𝑡),𝑂(𝑡),𝐷𝑆(𝑡), and 𝐷𝑂(𝑡) obtained with the NSFD scheme with 𝜓(ℎ)=ℎ, ℎ=0.2 and initial conditions are 𝑁(0)=0.999,𝐿(0)=0.001,𝑆(0)=0,𝑂(0)=0,𝐷𝑆(0)=0, and 𝐷𝑂(0)=0.

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 𝜌<1. 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 ℎ=0.25. The Runge-Kutta fourth-order scheme improves the result obtained by Euler scheme, but fails to converge for time step sizes ℎ>0.33 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 ℎ=0.2. In Figure 4, it is depicted that the NSFD scheme converges to the correct OEE point for a time step size ℎ=0.25, 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)

Figure 3 

Profiles of subpopulations 𝑁(𝑡),𝐿(𝑡),𝑆(𝑡),𝑂(𝑡),𝐷𝑆(𝑡), and 𝐷𝑂(𝑡). Solutions are obtained by Euler (dashed) and proposed NSFD scheme (line) with 𝜓(ℎ)=ℎ and ℎ=0.2 for system (2.1).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Profiles of subpopulations 𝑁(𝑡),𝐿(𝑡),𝑆(𝑡),𝑂(𝑡),𝐷𝑆(𝑡), and 𝐷𝑂(𝑡). Solutions are obtained by Euler scheme (dashed) and developed NSFD scheme (line) with 𝜓(ℎ)=ℎ and ℎ=0.25 for system (2.1).

(a)

(b)

(c)

(d)

Figure 5 

Solutions for subpopulation 𝐷𝑂 obtained by NSFD scheme and routines ODE45, ODE23s, ODE113, and ODE23 from MATLAB package (dashed), using 𝜎=0.0044379×10−5.

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.