# A Nonstandard Dynamically Consistent Numerical Scheme Applied to Obesity Dynamics

Authors: Rafael J. Villanueva; Abraham J. Arenas; Gilberto González-Parra
Journal: Journal of Applied Mathematics
 (2008)

---

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


---

## Body

## 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 Section2, 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, N(t) denotes the proportion of normal-weight
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, DS(t) denotes the proportion of overweight children
on diet, and DO(t) denotes the proportion of obese children on
diet. The model is represented by the following nonlinear system of ordinary
differential equations:(2.1)N'(t)=μ+εDS(t)−μN(t)−βN(t)[L(t)+S(t)+O(t)],L'(t)=βN(t)[L(t)+S(t)+O(t)]−[μ+γL]L(t),S'(t)=γLL(t)+φDS(t)−[μ+γS+α]S(t),O'(t)=γSS(t)+δDO(t)−[μ+σ]O(t),DS'(t)=γDDO(t)+αS(t)−[μ+ε+φ]DS(t),DO'(t)=σO(t)−[μ+γD+δ]DO(t), 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)
γ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. In this model, since the constant population is
normalized to unity one gets for all time t,(2.2)N(t)+L(t)+S(t)+O(t)+DS(t)+DO(t)=1.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.


## 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)ψ(h)=h+O(h2)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 f(t)∈C1(ℝ),
we define an equivalent derivative as(3.2)df(t)dt=f(t+h)−f(t)ψ(h)+O(ψ(h))ash⟶0,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)df(t)dt=limh→0{f(t+h)−f(t)ψ(h)+O(φ(h))}=limh→0f(t+h)−f(t)hlimh→0hψ(h)+limh→0O(ψ(h))=f˙(t).An example of a function ψ(h) that satisfies the above conditions is ψ(h)=1−e−h (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 byNn, Ln, Sn, On, DSn,
and DOn the approximations of N(nh), L(nh), S(nh), O(nh), DS(nh),
and DO(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, DS,
and DO of the model (2.1) is given by the following
implicit scheme:(3.4)DOn+1−DOnψ(h)=σOn−(δ+μ+γD)DOn+1,DSn+1−DSnψ(h)=−(μ+ϵ+φ)DSn+1+αSn+γDDOn+1,Nn+1−Nnψ(h)=μ−βNn+1(Sn+On+Ln)+ϵDSn+1−μNn+1,Ln+1−Lnψ(h)=βNn+1(Sn+On+Ln)−(γL+μ)Ln+1,Sn+1−Snψ(h)=γLLn+1+φDSn+1−αSn−(γS+μ)Sn+1,On+1−Onψ(h)=γSSn+1−σOn−μOn+1+δDOn+1, and after rearranging, we obtain
the following explicit form:(3.5)DOn+1=DOn+ψ(h)σOn1+ψ(h)(δ+μ+γD),(3.6)DSn+1=DSn+ψ(h)(αSn+γDDOn+1)1+ψ(h)(ϵ+μ+φ),(3.7)Nn+1=Nn+ψ(h)(μ+ϵDSn+1)1+ψ(h)(μ+β(Sn+On+Ln)),(3.8)Ln+1=Ln+ψ(h)βNn+1(Sn+On+Ln)1+ψ(h)(γL+μ),(3.9)Sn+1=Sn(1−ψ(h)α)+ψ(h)(γLLn+1+φDSn+1)1+ψ(h)(γS+μ),(3.10)On+1=On(1−ψ(h)σ)+ψ(h)(γSSn+1+δDOn+1)1+ψ(h)μ,  where the denominator function
is chosen as(3.11)ψ(h)=1−e−hMM,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].


## 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.
IfDOn>0, DSn>0, Nn>0, Ln>0, Sn>0, On>0, then DOn+1>0, DSn+1>0, Nn+1>0, Ln+1>0, Sn+1>0, On+1>0, for all n=0,1,2,…, since that 1−ψ(h)M>0;(2)
satisfy conservation law, that is, constant
population. From system (3.4) one gets that(3.12)DOn+1+DSn+1+Nn+1+Ln+1+Sn+1+On+1=DOn+DSn+Nn+Ln+Sn+On+μψ(h)1+μψ(h). It follows  by induction that if(3.13)DOn+DSn+Nn+Ln+Sn+On=1,thenDOn+1+DSn+1+Nn+1+Ln+1+Sn+1+On+1=1∀n=0,1,….Now, observe that the subpopulations must be calculated from scheme
(3.5)–(3.10) in the following order:(1)
select valuesDOn>0, DSn>0, Nn>0, Ln>0, Sn>0, On>0, such that DOn+DSn+Nn+Ln+Sn+On=1;(2)
determineDOn+1 from (3.5)) using DOn and On;(3)
calculateDSn+1 from (3.6) with DOn+1, DSn, and On;(4)
obtainNn+1 from (3.7) using DSn+1, Nn, Sn, On,
and Ln;(5)
findLn+1 from (3.8) with Nn+1, Sn, On,
and Ln;(6)
getSn+1 from (3.9) using Sn, Ln+1,
and DSn+1;(7)
computeOn+1 from (3.10) with the values of On, Sn+1 and DOn+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 𝒪(h2),
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.
ParameterValueμ0.00641025γS0.0030853α0.004068σ0.0044379γD0.0004189γL0.0089ε0.0027768φ0.12735δ0.15974β0.0222949


## 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:(4.1)OFE=(1,0,0,0,0,0),OEE=(0.2937,0.2971,0.2708,0.1267,0.0080,0.0033).The eigenvalues of the Jacobian evaluated at the OFE
point and at the positive OEE point are shown in Table2. 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.Table 2
The eigenvalues of the Jacobian of
system (2.1) evaluated at the OFE point and at the OEE point.
EigenvalueAt OFE pointAt the OEE pointλ1−0.170976−0.170977λ2−0.140644−0.140622λ3−0.0151898−0.0143034+0.00134272iλ40.0134083−0.0143034−0.00134272iλ5−0.00921949−0.00909934λ6−0.00641026−0.00641026


## 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.Profiles of subpopulationsN(t),L(t),S(t),O(t),DS(t), and DO(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, DS(0)=0.0249,
and DO(0)=0.0115.
(a)(b)(c)(d)(e)(f)Additionally in Figure2, 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).Profiles of subpopulationsN(t),L(t),S(t),O(t),DS(t),
and DO(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,DS(0)=0, and DO(0)=0.
(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 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
Spectral radius for different time
step sizeh of the Euler and NSFD numerical schemes.
Time stephρ (Euler)Eulerρ (NSFD)NSFD (ψ(h)=h)0.010.996667Convergence0.996678Convergence0.10.9935Convergence0.967742Convergence0.20.933333Convergence0.9375Convergence0.251.22302Divergence0.923077Convergence0.53.44604Divergence0.857143Convergence17.89207Divergence0.75Convergence1087.920Divergence0.230769ConvergenceTable3 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.Profiles of subpopulationsN(t),L(t),S(t),O(t),DS(t),
and DO(t).
Solutions are obtained by Euler (dashed) and proposed NSFD scheme (line) with ψ(h)=h and h=0.2 for system (2.1).
(a)(b)(c)(d)(e)(f)Profiles of subpopulationsN(t),L(t),S(t),O(t),DS(t),
and DO(t).
Solutions are obtained by Euler scheme (dashed) and developed NSFD scheme (line)
with ψ(h)=h and h=0.25 for system (2.1).
(a)(b)(c)(d)(e)(f)Solutions for subpopulationDO obtained by NSFD scheme and routines ODE45,
ODE23s, ODE113, and ODE23 from MATLAB package (dashed), using σ=0.0044379×10−5.
(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 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.



---

*Source: 640154-2008-12-30.xml*