A Nonstandard Dynamically Consistent Numerical Scheme Applied to Obesity Dynamics

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.


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.population into six subpopulations.In this model, N t denotes the proportion of normalweight children, L t denotes the proportion of latent children, that is, children with the habit of high consumption of unhealthy food but who are still normal weight, S t denotes the proportion of children with overweight, O t denotes the proportion of obese children, D S t denotes the proportion of overweight children on diet, and D O t denotes the proportion of obese children on diet.The model is represented by the following nonlinear system of ordinary differential equations:  In this model, since the constant population is normalized to unity one gets for all time t, 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 t.

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 The discretization is based on the approximations of the temporal derivatives by a first-order generalized forward scheme.If f t ∈ C 1 R , we define an equivalent derivative as where ψ h 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

3.3
An example of a function ψ h that satisfies the above conditions is ψ h 1−e −h see 4, 19 .

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 N n , L n , S n , O n , D n S , and D n O the approximations of N nh , L nh , S nh , O nh , D S nh , and D O nh , respectively, for n 0, 1, 2, . . ., and h the time step size of the scheme.Using 3.2 , the first-order numerical scheme to obtain solutions N, L, S, O, D S , and D O 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 as with M ≥ max{α, σ} in order to guarantee the positivity condition.However, positivity can also be guaranteed taking the traditional ψ h h if hM ≤ 1. Therefore throughout this paper when hM ≤ 1, we use ψ h h.
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 .

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: 2 satisfy conservation law, that is, constant population.From system 3.4 one gets that It follows by induction that if

3.13
Now, observe that the subpopulations must be calculated from scheme 3.5 -3.10 in the following order: O , D n S , and O n ; 4 obtain N n 1 from 3.7 using D n 1 S , N n , S n , O n , and L n ; 5 find L n 1 from 3.8 with N n 1 , S n , O n , and L n ; 6 get S n 1 from 3.9 using S n , L n 1 , and D n 1 S ; 7 compute O n 1 from 3.10 with the values of O n , S n 1 and D n 1 O .
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 O h 2 , confirming that the constructed NSFD scheme is first-order accurate.

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.

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 λ 4 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.

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

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

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

2
μ: inflow rate population in the system or inversely proportional to average time spent by an individual in the system, 3 γ L : rate at which a latent individual moves to the overweight subpopulation, 4 γ S : 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 γ D : rate at which an obese individual on diet becomes an overweight individual on diet.

Figure 1 :
Figure 1: Profiles of subpopulations N t , L t , S t , O t , D S t , and D O t .Solutions are obtained by the proposed NSFD scheme with ψ h h, where h 0.2 and initial conditions are N 0 0.462, L 0 0.194, S 0 0.2176, O 0 0.09, D S 0 0.0249, and D O 0 0.0115.

Figure 2 :
Figure 2: Profiles of subpopulations N t , L t , S t , O t , D S t , and D O t obtained with the NSFD scheme with ψ h h, h 0.2 and initial conditions are N 0 0.999, L 0 0.001, S 0 0, O 0 0, D S 0 0, and D O 0 0.

Figure 3 :
Figure 3: Profiles of subpopulations N t , L t , S t , O t , D S t , and D O t .Solutions are obtained by Euler dashed and proposed NSFD scheme line with ψ h h and h 0.2 for system 2.1 .

Figure 4 :
Figure 4: Profiles of subpopulations N t , L t , S t , O t , D S t , and D O t .Solutions are obtained by Euler scheme dashed and developed NSFD scheme line with ψ h h and h 0.25 for system 2.1 .

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

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

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