Analysis of a nonstandard computer method to simulate a nonlinear stochastic epidemiological model of coronavirus-like diseases.

BACKGROUND AND
OBJECTIVE: We propose a nonstandard computational model to approximate the solutions of a stochastic system describing the propagation of an infectious disease. The mathematical model considers the existence of various sub-populations, including humans who are susceptible to the disease, asymptomatic humans, infected humans and recovered or quarantined individuals. Various mechanisms of propagation are considered in order to describe the propagation phenomenon accurately.
METHODS: We propose a stochastic extension of the deterministic model, considering a random component which follows a Brownian motion. In view of the difficulties to solve the system exactly, we propose a computational model to approximate its solutions following a nonstandard approach.
RESULTS: The nonstandard discretization is fully analyzed for positivity, boundedness and stability. It is worth pointing out that these properties are realized in the discrete scenario and that they are thoroughly established herein using rigorous mathematical arguments. We provide some illustrative computational simulations to exhibit the main computational features of this approach.
CONCLUSIONS: The results show that the nonstandard technique is capable of preserving the distinctive characteristics of the epidemiologically relevant solutions of the model, while other (classical) approaches are not able to do it. For the sake of convenience, a computational code of the nonstandard discrete model may be provided to the readers at their requests.

Introduction

Mathematical epidemiology is nowadays an area of active investigation which spans from the mathematical understanding of the elementary principles describing the propagation of diseases to the computational simulation of the propagation of an epidemic at the micro and macroscopic levels [1]. To that end, various approaches may been followed. Probably, the most popular models to describe the propagation of diseases rely on systems of first-order ordinary differential equations [2]. Many such models have been applied to the investigation of infectious diseases, including susceptible-infected-susceptible nonlinear models that were employed to describe nonlinear rates of incidence along with duplicate hypotheses on the epidemic [3], susceptible-infected-recovered models that are used in considering the knowledge on the danger of a disease [4], susceptible-exposed-infected-recovered systems which forecast the spread of a disease through global laws of reaction [5], susceptible-exposed-infected-quarantined-recovered systems that take into consideration adjusted incidence and imperfect vaccinations [6], dynamical models for the epidemiology of diarrhea and their simulation considering multiple disease carriers [7], and even some fractional epidemiological models for the propagation of new computer viruses [8], just to mention some models. As we pointed out, these approaches hinge on the use of systems of ordinary differential equations. However, there are various other approaches which use stochastic modeling to propose new mathematical models for the propagation of infectious diseases among human populations [9], [10].

During the last year, the investigation in mathematical epidemiology has suffered a tremendous development in view of the outbreak of the new coronavirus COVID-19 [11]. Indeed, many effects of the propagation of this disease have surfaced on the daily development of the human activity. From the merely medical point of view, the consequences of this infection have been investigated thoroughly worldwide under a vast number of perspectives [12], [13], [14], but a complete understanding of the transmission of this disease is still lacking, and a solid mechanism to control the propagation of this infection is an open avenue of investigation. Here, it is worth pointing out that there are many other questions on the effects of this disease in human populations. As an example, researchers still discuss about the economic consequences of our current decisions to control the propagation of COVID-19 [15]. How will the micro- and macro-economy be affected by the current government decisions to control the pandemic? How will those decision impact on the migration of human population worldwide and, ultimately, how will they affect the global economy [16]? How has the social psychology changed during this crisis as a consequence of the decisions taken by the authorities [17]. Indeed, many question remain unanswered to this day, and more questions will be raised as the pandemic propagates without an affirmative solution. Still, the research in this topic will be an important avenue of research from every possible aspect in our modern human condition.

From the merely epidemiological point of view, the empirical data on the propagation of diseases have shown that a non-deterministic behavior is invariably observed [18]. In view of this fact, various stochastic models have been proposed to describe the propagation of infectious diseases in human populations. As examples, some real-time forecasting of infectious disease dynamics have been proposed for the investigation of stochastic semi-mechanistic models [19] (together with some stochastic susceptible-infectious-recovered systems with horizontal and vertical transmissions [20]), stochastic susceptible-infected-susceptible epidemic models with nonlinear incidence rate [21], stochastic susceptible-infected-recovered-susceptible epidemic models with saturated incidence rate and transfer from infectious to susceptible [22], susceptible-infected-susceptible nonlinear stochastic epidemiological models with incidence rate and double epidemic hypothesis [23], stochastic susceptible-exposed-infected-quarantined-recovered epidemic systems with quarantined-adjusted incidence and imperfect vaccination [6] and stochastic systems describing the spread of diseases across a heterogeneous social environment [24], just to mention some reports available in the specialized literature.

In the present work, we will depart from a mathematical model which may be able to describe the propagation of infectious diseases like COVID-19. The mathematical model is a deterministic system which considers the presence of various sub-populations and their nonlinear interactions. In order to incorporate a random component, we will consider a general extension of this system in the form of a stochastic differential equation. The extension will consider states and transition probabilities, and it will be provided in its most general form. Two components will be considered in the system: a deterministic and a probabilistic term. It is worth pointing out that the probabilistic component will be governed by a Brownian motion, and various deterministic parameters will be incorporated in the mathematical model. For practical purposes, a particular form of that system will be considered, and we will propose a nonstandard technique to solve the stochastic model. For comparison purposes, other standard discretizations will be studied, including a stochastic Euler and a stochastic Runge–Kutta methodologies. The deterministic form of the nonstandard approach will be fully analyzed for the existence and uniqueness of positive and bounded solutions. We will show that the scheme preserves the constant solutions, and the stability of the methodology will be thoroughly proved. Some simulations will be provided for illustration purposes in this manuscript.

This work is organized as follows. For the remainder of the present section, we present the epidemiological model under investigation. The model is a SAIR system which structures a human population into various sub-populations: individuals who are susceptible to the disease, asymptomatic individuals (that is, individuals who have the virus but do not present any symptoms and could spread the virus), infected humans (individuals who have the virus, present symptoms and may spread the virus) and recovered individuals. We will calculate the reproductive number therein, along with the equilibrium points considering disease-free and endemic populations. A general stochastic model will be derived therein, and a particular form of that system will be introduced in that stage to be investigated in the present report. In Section 2, we will present various computational methods to solve the differential problem under investigation. Three of them will be standard (and classical) discretizations, while one of them will follow a nonstandard (and non-traditional) approach. The non-stochastic form of the nonstandard method will be theoretically analyzed in Section 3. In particular, we will prove the existence and uniqueness of positive and bounded solutions of the nonstandard scheme along with its stability. Moreover, we will show that scheme is able to preserve the equilibrium solutions of the model. In that section, we will provide some simulations to show the performance of the computational model, and we will close this work with a section of concluding remarks.

As we mentioned in the previous paragraph, in this report, we will consider a population of human individuals in which an infectious disease spreads out. Throughout, we will assume that each of the individuals of the population belongs to exactly one of the following sub-populations: susceptible, asymptomatic, infected and recovered/quarantined humans. In the present work, we let represent a finite time period, and let denote the temporal variable. For each such we will let represent the total number of susceptible humans at the time . In similar fashion, will denote the number of asymptomatic humans, will be the number of infected individuals and will represent the number of recovered or quarantined humans at the time . Notice that, as functions of each of the functions and is non-negative. Moreover, the quantity provides the population size at the time for any number .

To derive a stochastic model for the spreading of the disease in our population, various assumptions will be imposed on the dynamics of each of the sub-populations. To start with, we will assume that the rate of birth of individuals is a constant and that all newborns are susceptible. For the sake of simplification, the rate of death will be also the positive number and it will be constant among all the sub-populations. In addition, let be the constant rate at which susceptible humans become asymptomatic due to contact with infected people, and let be a constant which represents the bilinear incidence rate of infected individuals. With these parameters, we will suppose that the sub-population of susceptible individuals decreases at a rate which is directly proportional to the interactions between the susceptible and infected populations, and inversely proportional to the bilinear incidence rate. More precisely, we will suppose that the following identity is satisfied:

In addition to the susceptible individuals of the population who become asymptomatic subjects, we assume that the sub-population of asymptomatic humans decreases at a rate equal to . Here, is the rate at which asymptomatic humans become infected, denotes the mortality rate of asymptomatic individuals due to the disease, and is the immunity rate of asymptomatic individuals. In light of this fact, the dynamics of the rate of change in the sub-population of asymptotic individuals is governed by the ordinary differential equation

In similar fashion, we will suppose that the sub-population of infected humans increases proportionally to the amount of asymptotic individuals, with constant of proportionality equal to . This sub-population also decreases at a rate directly proportional to the amount of infected individuals, with constant of proportionality equal to . Here, is the rate of vaccination, quarantine or treatment, and is the death rate of infected humans due to virus. Finally, the sub-population of recovered individuals increases proportionally to the number of infected individuals with a constant of proportionality and increases proportionally to the amount of asymptomatic members of the population. In summary, the following ordinary differential equations hold: For convenience, Figure 1 shows a diagram depicting the dynamics of propagation of the disease in the mathematical model (1.1)–(1.4). Moreover, Table 1 provides a summary of the physical meaning of the parameters used in this model.

Fig. 1

Diagram depicting the dynamics of propagation of the disease in the mathematical model (1.1)–(1.4).

Table 1

Table of the physical meanings of the parameters used in the stochastic model investigated in this work.

Parameter	Physical meaning
μ	Natural birth or death rate
β_0	Rate at which susceptible humans become asymptomatic
α	Bilinear incidence rate
σ	Rate at which asymptomatic humans become infected
δ	Mortality rate of asymptomatic humans due to virus
ϵ	Immunity rate of asymptotic humans
r	Rate of vaccination, quarantined or treated individuals
d	Death rate of infected humans due to virus

The reproductive number of the population system(1.1)–(1.4)is given by

The proof is straightforward. □

In this work, we will study analytically and numerically the system (1.1)–(1.4). It is obvious that this is a system of nonlinear ordinary differential equations. Some simple algebraic calculations show that the equilibrium point of these system are Here, the first equilibrium solution is the trivial solution. Meanwhile, the second equilibrium is the disease-free equilibrium solution, and the last is the disease-present equilibrium point, in which case we observe the following nomenclature:

In order to provide a stochastic generalization of this model, we will use the nomenclatureWe will consider a uniform partition of the interval consisting of subintervals, where . Define then the temporal step-size . Moreover, for the remainder, we will use the transition probabilities summarized in Table 2 . It is worth pointing out that we are employing the transition probabilities used in [25]. Under those conditions and using simple algebraic arguments, it is easy to check that the expected value and the variance transition phases are given byandrespectively. In this expression, we agree that

Table 2

Table of transitions () and their respective probabilities (), for each representing possible individual changes. These values have been taken from [25] for convenience.

i	T_i	P_i
1	[0,0,0,0]^⊤	μτ
2	[−1,1,0,0]^⊤	β_2SIτ1+αI^2
3	[−1,0,0,0]^⊤	μστ
4	[0,−1,1,0]^⊤	σAτ
5	[0,−1,0,1]^⊤	ϵAτ
6	[0,−1,0,0]^⊤	(δ+μ)Aτ
7	[0,0,−1,0]^⊤	(d+μ)Iτ
8	[0,0,−1,1]^⊤	rIτ
9	[0,0,0,−1]^⊤	μRτ

For the remainder of this work, we will define the stochastic drift and the stochastic diffusion of the process at the time respectively, by the formulas Using these conventions, the stochastic generalization of the deterministic system (1.1)–(1.4) is given byIn this context, will represent the Brownian motion, and we will consider initial conditions of the formIt is worth pointing out here that our numerical simulations will make use of the initial data .

Before closing this section, it is important to note that the present manuscript will focus on the following stochastic non-integrable extension of the deterministic system (1.1)–(1.4):In these expressions, we have obviated the dependence of and with respect to . Moreover, is a constant perturbation parameter which has been considered here to stochastically extend the initial deterministic system. If then we retrieve the original deterministic model. Otherwise, the stochastic model considers stochastic components that depend on the interaction between the susceptible and the infected members of the population.

Methods

The present section is devoted to introduce some standard and a non-standard stochastic numerical methodologies to approximate the solutions of the stochastic model (1.29). To that end, we will agree that and for each . Let and consider a uniform partition of the temporal interval with partition norm equal to . The nodes of that partition are represented byfor each . Needless to mention that for each . Moreover, we will agree that whenever and . Also, we setIt is obvious that each has a normal distribution with mean equal to zero and variance equal to 1.

The nonstandard method will be designed following the approach proposed by R. E. Mickens [26]. It is worth pointing out that such approach has the advantage of providing highly accurate approximations with a low computational cost. Also, we must point out that the preservation of various interesting features of relevant solutions may be guaranteed by using the nonstandard approach. For example, it has been used to guarantee the boundedness of the computational solutions for damped nonlinear wave equations [27], the monotonicity of numerical solutions for generalized Burgers–Huxley-type equations [28], the numerical dissipation of energy of a class of dissipative nonlinear wave equations [29], [30], the numerical preservation of positivity of a computational method to solve a susceptible-exposed-infected-recovered reaction diffusion model [31], to solve a nonlinear model for the propagation of COVID-19 using a numerically efficient and positivity-preserving computer scheme [32] and to preserve the positivity and the boundedness of some algorithms to solve fractional nonlinear epidemic models for HIV/AIDS transmission [33].

Euler method

The first numerical methodology used in this work to approximate the solutions of (1.29) is the stochastic version of the well-known Euler’s method, which is the iterative technique described by the following system of discrete difference equations, which are valid for each :

Runge–Kutta method

The second standard numerical approach used in this work is based on the fourth-order Runge–Kutta family of methods to solve systems of ordinary differential equations. In its stochastic form, this family is also a four-stage technique, and we describe each of those stages next.

Stage 1

Stage 2

Stage 3

Stage 4

Final stage

Euler–Maruyama method

The Euler–Maruyama method can be employed also to solve the stochastic model under investigation. Using the numerical conventions introduced in this work, the iterative formula of the Euler–Maruyama model is provided byfor each .

Nonstandard method

The methods introduced in the previous paragraphs are traditional discretization schemes for stochastic differential equations. The fourth method considered in this work is a nonstandard technique in the sense of R. E. Mickens [26]. In that context, we let and consider the explicit nonstandard discretization of (1.29) given by the iterative system

Results

We will focus our attention on the nonlinear computational model (2.10) and establish its main computational properties. In particular, we will prove the existence and uniqueness of non-negative deterministic solutions, the boundedness of the region of those solutions and the preservation of the equilibria of the continuous-time system. Moreover, we will prove the stability of the discrete model by examining the moduli of the eigenvalues of the discrete system.

Our first two results concern the existence and the uniqueness of positive and bounded solutions of the deterministic form of the scheme (2.10).

Positivity

If all the entries ofbelong tothen the deterministic form of the system(2.10)has a unique solutionwhere all the components ofare positive numbers, for each.

The proof uses mathematical induction. Note that the conclusion is true for in view of the hypotheses. Now, assume that the conclusion is true for some . Observe that the right-hand side of each of the equations in (2.10) are positive when and when all the parameters are likewise positive. It follows that all the entries of are positive numbers. The conclusion follows now by induction. □

Boundedness

Letbe a solution of the deterministic form of model(2.10)for non-negative initial conditions which satisfyThen, for eachthe following is satisfied:

Beforehand, notice that the solution of the discrete problem (2.10) are non-negative according to Theorem 2 Moreover, for simplification purposes, let for each . Now, the conclusion of this result is true when so let us suppose that it is valid for some . Using algebraic arguments, we can rewrite the system (2.10) as follows:Adding these four equations, simplifying algebraically and regrouping, it follows thatThe conclusion of the theorem is reached by induction. □

Equilibrium solutions

The deterministic form of the numerical model(2.10)has the same equilibrium solution of the mathematical model(1.1)–(1.4).

The proof of this result is readily reached after substituting the equilibrium points and into the deterministic form of (2.10) in the general, the disease-free and the endemic cases, respectively. □

In order to establish the stability of the nonstandard finite-difference scheme, we will require the following lemma.

Brauer et al.[34]

Letand suppose thatare the roots of. Thenforif and only if the following conditions are satisfied:

and

.

The equilibrium points of the deterministic form of the nonstandard finite-difference scheme(2.10)are stable.

For the sake of convenience, let us defined the functions by the following expressions: for each . It is clear that the finite-difference scheme (2.10) can be equivalently rewritten asThe Jacobian of this system is given byIt is easy to check thatMoreover, the following identities are satisfied:To prove the stability of the equilibrium points it suffices to show that all the eigenvalues associated to the corresponding Jacobian matrix lie within the open unit circle of with center at the origin. Consider the disease-free equilibrium point, which is given by . In that case, the expression of the Jacobian reduces significantly. Moreover, after algebraic calculations and simplifications, the characteristic polynomial associated to the Jacobian matrix takes the formwhere Obviously, has modulus less than 1. Now, algebraic arguments show that the conditions (i), (ii) and (iii) of the previous lemma are satisfied for the polynomial which implies that the eigenvalues of that polynomial have modulus less than 1 also. It follows that the eigenvalues of the Jacobian matrix associated to the equilibrium point have modulus less than 1, which implies that this equilibrium is a stable solution of (2.10). The stability of the other equilibrium points can be established in similar fashion. □

For illustration purposes, Figure 2 shows the largest eigenvalue of the Jacobian matrix (3.10) for the endemic equilibrium point as a function of . The graphs confirm that the largest eigenvalue of the Jacobian matrix has modulus less than one, for all . This is in obvious agreement with the conclusion of the previous proposition. Here, we must point out that the spectral radius is quite constant if and approximately equal to 0.983.

Fig. 2

Value of the largest eigenvalue of the Jacobian matrix (3.10) for the endemic equilibrium point as a function of . The graphs confirm that the largest eigenvalue of the Jacobian matrix has modulus less than one, for all .

Before closing this section, we provide some illustrative simulations on the performance of the computer method (2.10), especially when we compare it against the stochastic Euler and Runge–Kutta approaches. To that end, we will fix the model parameters in Table 3 (see [10]), and use the initial data described in the previous section. For the sake of concreteness, we may think of the disease being the coronavirus COVID-19. Under these circumstances, Figure 3 shows the numerical solution for the number of infected people using the stochastic nonstandard scheme, along with the stochastic Euler (top row), stochastic Runge–Kutta (middle row) and the Euler–Maruyama (bottom row) schemes [35].

Table 3

Table of the parameter values of the stochastic model investigated in this work, which have been used in the computer simulations (see [10]).

Parameter	Value
μ	0.1
β_0 (disease-free)	0.5464
β_0 (endemic)	0.9464
α	0.01
σ	0.05
δ	0.001
ϵ	0.03
d	0.02
σ_1	0.5

Fig. 3

Dynamics of the number of infected people with respect to time using the stochastic model (1.29) with the model parameters in Table 3. In all the cases, we employed the stochastic nonstandard (NSFD) scheme, along with the stochastic Euler (top row), stochastic Runge–Kutta (middle row) and the Euler–Maruyama (bottom row) schemes. Our results confirm that the nonstandard method is the only approach considered in this work which preserves the positivity and the boundedness of the numerical approximations. For convenience, the values of used to produce each figure are mentioned on each graph.

To be more precise, the three rows of Figure 3 compare the stochastic nonstandard finite-difference method proposed in this work against each of the three methods mentioned above. On the left column, we present simulations with step-sizes for which the properties of positivity and boundedness are satisfied. Meanwhile, the right columns show simulations with a different value of for which the nonstandard discretization still preserves the positivity and the boundedness, while the other numerical methods do not. These simulations are strong evidence on the capability of the nonstandard scheme to preserve the positivity and the boundedness of approximations. These features obviously represent an advantage of the nonstandard approach over the other methodologies considered in this work. Indeed, the results confirm that the only method to preserve the positivity and the boundedness of solutions is the stochastic nonstandard methodology.

For convenience, Figure 4 (a) shows the effect of the quarantine parameter on the infected individuals, and Figure 4(b) depicts the effect of the quarantine parameter on the reproduction number. Here, we have employed the same model parameters used to obtain the simulations in Figure 3. In the former graph, various values of the parameter are considered. Independently of the value of this parameter, the number of infected humans is approximately a decreasing function of time. Note that relatively small values of the quarantine parameter result in a number of infected humans which converges asymptotically to a positive value. On the other hand, the population of infected individuals tends asymptotically to zero when . Finally, notice that the bottom graph shows that the reproduction number of the population is a decreasing function of the quarantine parameter.

Fig. 4

Effect of the quarantine parameter on (a) the number of infected individuals and (b) the reproduction number. The stochastic model and parameters are the same as those used in Figure 3. For convenience, the value of used to produce figure (a) is mentioned on that graph.

Before closing this section, we must recall that most of the numerical works on mathematical epidemiology used explicit methods (like the Euler and Runge–Kutta approaches) for the discretization of stochastic epidemic models. Unfortunately, those methods fail to preserve the essential properties of the solutions (like positivity, boundedness, monotonicity and dynamical consistency), and strongly depend upon the temporal step-size. As our simulations show, those methods may produce negative and unbounded solutions which have no significance in population dynamics. Indeed, a negative value of any sub-population is meaningless. The positivity of the sub-populations is essential in all of those scenarios (the importance of these properties is discussed in detail by R. E. Mickens [26] in his works on dynamic consistency). In the present work, we have shown that our proposed nonstandard method preserves all of these properties and is unconditionally convergent. This is the main advantage of the proposed nonstandard scheme: it remains consistent with the biological nature of continuous dynamical system for all choices of the discretization parameter. Moreover, the nonstandard approach has lower computational cost when compared against the Euler and the Runge–Kutta schemes.

To support our claim that the nonstandard scheme has lower computational cost than the standard approaches, Figure 5 provides a computational efficiency analysis of those classical finite-difference schemes and the nonstandard methodology. To perform this study, we consider the same problem as that used to obtain Figure 3, for various values of the temporal step-size. For each of the methods, the exact solution was estimated calculating the numerical solution for the computational parameter value . After that, the numerical solutions were obtained for those methods and various temporal step-sizes. The error was calculated then using the standard Euclidean norm. It is obvious that the standard approaches are computationally more expensive than the nonstandard scheme investigated in this work. According to the results, the Runge–Kutta scheme is more efficient than the Euler discretization, as expected. In turn, our present nonstandard approach is actually more efficient than the Runge–Kutta technique. This is justified by the fact that the latter scheme requires many algebraic calculations at each iteration, in spite that it has convergence order equal to four. In any case, the present study revels that the nonstandard methodology is computationally more efficient.

Fig. 5

Efficiency analysis comparing the nonstandard finite-difference (NSFD) methodology (solid line) of the present work with the standard Euler (dotted) and Runge–Kutta (dashed) methods in the deterministic scenario. We used the problem investigated in Figure 3 with various values of the computational parameter to determine the absolute error (in the Euclidean norm) and the CPU times. For each of the methods, the exact solution was estimated calculating the numerical solutions for the computational parameter .

Discussion

In this manuscript, we considered a deterministic mathematical model which has the potential of being able to model the spreading of infectious diseases. In particular, the mathematical model investigated in this work may potentially describe the propagation of the new coronavirus COVID-19. The mathematical model is a deterministic system which considers the presence of various sub-populations and their nonlinear interactions. In order to incorporate a random component, we considered a general extension of this system in the form of a stochastic differential equation. The extension studied in this report considers the presence of various states and transition probabilities, and it is given herein in its most general form. Two general terms are considered in the system: a deterministic and a probabilistic term. The probabilistic component is governed by a Brownian motion, and various deterministic parameters are incorporated in the mathematical model. For practical purposes, a particular form of that system was considered in this report, and we proposed a nonstandard technique to solve the stochastic model. For comparison purposes, various standard discretizations were studied, namely, a stochastic Euler, a stochastic Runge–Kutta and the Euler–Nurayama methodologies. The deterministic form of the nonstandard approach was fully analyzed for the existence and uniqueness of positive and bounded solutions. We showed that the nonstandard scheme preserves the constant solutions, and the stability of the methodology was thoroughly proved. Some simulations were provided for illustration and comparison purposes. Our numerical results confirm the theoretical findings on the nonstandard approach.

Before closing this manuscript, it is important to mention that the approach employed in the present work may be extended to a wide variety of problems in the sciences and engineering. As one of the reviewers pointed out, our current approach may be applied to conservative systems with non-classical solutions to strongly coupled problems [36], systems of partial differential equations associated to fuzzy differential equations [37], strongly coupled systems with moving density constraints in traffic flow [38], the solution of Volterra integro-differential equations in structural analysis [39], problems of moving bottlenecks in traffic flow [40], the solution of fluid dynamic models for supply chains [41], problems associated to the cardiovascular system [42], supply chains for the production of wine and other products [43], [44], among many other problems [45], [46]. Also, we must point out that the present work imposes various assumptions on the dynamics of each of the sub-populations involved in our continuous model. As example, we are assuming that the rate of birth is equal to the rate of death. Other studies have imposed such kind of hypotheses [10] for simplification purposes. Indeed, the theoretical analysis of the continuous and discrete models turns out to be relatively easier. Obviously, these assumptions represent a limitation in our theoretical study. However, the numerical results provide strong evidence on the advantages in employing a nonstandard approach over the traditional methodologies. Moreover, notice that the extension of the nonstandard technique to account for more general assumptions is relatively straightforward, in spite of the fact that the analysis of the method could be more difficult.

Funding

The corresponding author (J.E.M.-D.) wishes to acknowledge the financial support from the National Council for Science and Technology of Mexico (CONACYT) through grant A1-S-45928.

Ethics approval

This manuscript does not contain any studies with human participants performed by any of the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author, J.E.M.-D., upon reasonable request.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.