Next Article in Journal
An Iteration Function Having Optimal Eighth-Order of Convergence for Multiple Roots and Local Convergence
Next Article in Special Issue
The Convergence Analysis of a Numerical Method for a Structured Consumer-Resource Model with Delay in the Resource Evolution Rate
Previous Article in Journal
On Fundamental Solution for Autonomous Linear Retarded Functional Differential Equations
Previous Article in Special Issue
Lp-Solution to the Random Linear Delay Differential Equation with a Stochastic Forcing Term
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mean Square Convergent Non-Standard Numerical Schemes for Linear Random Differential Equations with Delay

1
Instituto Universitario de Matemática Multidisciplinar, Building 8G, Access C, 2nd Floor, Universitat Politècnica de València, Camino de Vera s/n, 46022 Valencia, Spain
2
Department of Applied Mathematics, University of Alicante, Apdo. 99, 03080 Alicante, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(9), 1417; https://doi.org/10.3390/math8091417
Submission received: 28 July 2020 / Revised: 21 August 2020 / Accepted: 21 August 2020 / Published: 24 August 2020
(This article belongs to the Special Issue Models of Delay Differential Equations)

Abstract

:
In this paper, we are concerned with the construction of numerical schemes for linear random differential equations with discrete delay. For the linear deterministic differential equation with discrete delay, a recent contribution proposed a family of non-standard finite difference (NSFD) methods from an exact numerical scheme on the whole domain. The family of NSFD schemes had increasing order of accuracy, was dynamically consistent, and possessed simple computational properties compared to the exact scheme. In the random setting, when the two equation coefficients are bounded random variables and the initial condition is a regular stochastic process, we prove that the randomized NSFD schemes converge in the mean square (m.s.) sense. M.s. convergence allows for approximating the expectation and the variance of the solution stochastic process. In practice, the NSFD scheme is applied with symbolic inputs, and afterward the statistics are explicitly computed by using the linearity of the expectation. This procedure permits retaining the increasing order of accuracy of the deterministic counterpart. Some numerical examples illustrate the approach. The theoretical m.s. convergence rate is supported numerically, even when the two equation coefficients are unbounded random variables. M.s. dynamic consistency is assessed numerically. A comparison with Euler’s method is performed. Finally, an example dealing with the time evolution of a photosynthetic bacterial population is presented.

1. Introduction

Modeling physical systems for which the future state depends on history due to hereditary characteristics, such as aftereffects or time lags, usually requires the use of delay differential models. The delay may be discrete or continuous, depending on whether a specific or complete past information is used. The inclusion of a delay requires specific techniques for the theoretical study of the differential model [1,2,3,4]. In practice, delay differential models play a key role in different scientific and technical fields [5,6,7,8,9,10].
In the context of delay differential equations, the construction of non-standard finite difference (NSFD) numerical schemes has not been much explored. Historically, NSFD schemes were developed by Mickens in the years 1994 and 2000 [11,12], together with a later edited book in 2005 [13]. Mickens observed that traditional standard finite difference schemes may be modified, on the basis of exact numerical schemes for basic ordinary differential equations, so that the essential properties of the governing continuous model are mimicked [14]. Until relatively recently, NSFD schemes were successfully designed and applied for ordinary, partial and fractional differential equations [15]. However, delay differential equations have not been addressed in detail.
Recently, [16] proposed a NSFD scheme for the general linear delay problem
x ( t ) = α x ( t ) + β x ( t τ ) , t > 0 , x ( t ) = f ( t ) , τ t 0 ,
( τ > 0 ) from an exact scheme on the whole domain, providing high order of accuracy and consistent dynamical behavior with simple computational properties. Such approach was extended to the non-scalar case in [17].
In modeling, the variability of data, due to limited knowledge and fluctuation of the process under study, lack of information, bad calibration machines, etc., gives rise to variability in the model coefficients. Therefore, for a more realistic description of the process, coefficients should be regarded as random quantities on an abstract probability space. When the coefficients are random variables and regular stochastic processes, the solution to the model becomes a differentiable stochastic process, whose realizable trajectories solve the deterministic version of the model. A common treatment of random differential models uses mean square (m.s.) calculus [18,19,20,21,22,23,24]. Of special importance is the computation of the mean and the variance of the solution stochastic process, or even its probability density function.
We are interested on delay random differential equations. Specifically, the randomization of (1) as
x ( t , ω ) = α ( ω ) x ( t , ω ) + β ( ω ) x ( t τ , ω ) , t > 0 , ω Ω , x ( t , ω ) = f ( t , ω ) , τ t 0 , ω Ω .
Here α and β are random variables and f is a stochastic process on a complete probability space ( Ω , F , P ) , where Ω is the sample space formed by the outcomes ω Ω , F is the σ -algebra of events, and P : F [ 0 , 1 ] is the probability measure. The solution x is a differentiable stochastic process.
Only recently, a theoretical study on delay random differential equations was started. General delay random differential equations were analyzed in the m.s. sense in [25], with the goal of extending some of the existing results on random differential equations with no delay from the book [18]. Problem (2) was solved in the m.s. sense in [26], and later generalized to equations with a random forcing term in [27]. On the other hand, in [28] the authors studied (2), but considered the solution in the sample-path sense and computed its probability density function via the random variable transformation technique, for certain forms of the initial condition process.
In this paper, we are concerned with computational aspects of delay random differential equations. Standard finite difference methods have already been applied to random ordinary, partial and fractional differential equations, by establishing the m.s. convergence, and even the convergence of densities, of the numerical discretizations towards the stochastic process solution [29,30,31,32,33,34]. Here we aim at extending the NSFD method from [16] to (2), by assessing the m.s. convergence of the discretizations. This permits approximating the expectation and the variance of the solution with high accuracy, whenever computationally feasible.
The organization of this paper is the following. In Section 2, the main results on m.s. calculus are exposed. The material for this section is essentially taken from [18]. In Section 3, the NSFD numerical scheme from [16] is presented. The randomization of the scheme, its m.s. convergence and its usefulness for approximating moments are discussed in Section 4. Illustration of the theory with numerical examples is conducted in Section 5. Finally, Section 6 draws the main conclusions.

2. M.s. Calculus

We are interested in second order real random variables y : Ω R , satisfying
E [ y 2 ] = Ω y ( ω ) 2 d P ( ω ) < .
We refer the reader to ([18], Ch. 4), [35]. The set of these random variables is a Hilbert space, denoted as L 2 ( Ω ) and endowed with the inner product y 1 , y 2 = E [ y 1 y 2 ] . This inner product gives rise to the norm y 2 = ( E [ y 2 ] ) 1 / 2 . By Cauchy–Schwarz inequality ([18], p. 19) E [ | y 1 y 2 | ] y 1 2 y 2 2 . Random variables in L 2 ( Ω ) are characterized by having finite variance:
V [ y ] = E [ ( y E [ y ] ) 2 ] = E [ y 2 ] ( E [ y ] ) 2 < .
This is one of the principal reasons for working with second order random variables, since the main statistical information for uncertainty quantification, namely the average value and the dispersion, are well-defined.
Given a stochastic process { z ( t ) : t I R } , it is of second order if the random variable z ( t ) is of second order, for all t I . By Cauchy–Schwarz inequality, it is straightforward to check that a second order stochastic process possesses a correlation function, E [ z ( t 1 ) z ( t 2 ) ] .
Convergence in L 2 ( Ω ) is defined through its norm · 2 : a sequence of random variables { y n } n = 1 converges to y in L 2 ( Ω ) if lim n y n y 2 = 0 . This is referred to as m.s. convergence.
M.s. convergence preserves the convergence of the expectation and the variance. This is a key fact. In general, if { x n } n = 1 and { y n } n = 1 are two sequences of second order random variables such that x n x and y n y as n in the m.s. sense, then E [ x n y n ] E [ x y ] ([18], p. 88).
In the particular case that { x n } n = 1 is a sequence of second order random variables such that its mean and its variance tend to zero, i.e., E [ x n ] 0 and V [ x n ] 0 as n , then { x n } n = 1 is m.s. convergent to zero, since ( x n 2 ) 2 = E [ ( x n ) 2 ] = V [ x n ] + ( E [ x n ] ) 2 0 as n . The converse is also true.
M.s. convergence gives rise to m.s. calculus, where continuity, differentiability and Riemann integrability of a stochastic process are naturally defined by taking m.s. limits in the classical definitions. A stochastic process { z ( t ) : t I R } is m.s. continuous at t 0 I if z ( t ) z ( t 0 ) as t t 0 in the m.s. sense. It is m.s. differentiable at t 0 I if lim h 0 z ( t 0 + h ) z ( t 0 ) h exists in the m.s. sense, which is denoted as z ( t 0 ) . Finally, z ( t ) is m.s. Riemann integrable on an interval [ a , b ] I if there exists a sequence of partitions { P n } n = 1 with mesh tending to 0, P n = { a = t 0 n < t 1 n < < t r n n = b } , such that for any choice of points s i n [ t i 1 n , t i n ] , i = 1 , , r n , the limit lim n i = 1 r n z ( s i n ) ( t i n t i 1 n ) exists in the m.s. sense, and it is denoted as a b z ( t ) d t .
The following important properties will be used: m.s. continuity on an interval implies m.s. Riemann integrability ([18] Section 4.5.1 (1)), a b z ( t ) d t 2 a b z ( t ) 2 d t for any m.s. Riemann integrable process z ( t ) ([18] Section 4.5.1 (3)), and the fundamental theorem of m.s. calculus ([18] Section 4.5.1 (5), (6)).
Finally, we mention that the essential supremum norm is defined as
y = inf { C 0 : | y | C almost surely } .
The set of random variables satisfying y < gives rise to the Banach space L ( Ω ) . Obviously, for two random variables y 1 L ( Ω ) and y 2 L 2 ( Ω ) , it holds y 1 y 2 2 y 1 y 2 2 < .

3. NSFD Methods for Linear Deterministic Differential Equations with Delay

Based on the explicit solution to (1),
x ( t ) = f ( 0 ) k = 0 m 1 β k ( t k τ ) k k ! e α ( t k τ ) + k = 0 m 2 β k + 1 k ! τ 0 ( t ( k + 1 ) τ s ) k e α ( t ( k + 1 ) τ s ) f ( s ) d s + β m ( m 1 ) ! τ t m τ ( t m τ s ) m 1 e α ( t m τ s ) f ( s ) d s ,
for ( m 1 ) τ < t m τ , m 1 , an exact numerical difference scheme for (1) is obtained in [16,17]. It is detailed in the following theorem.
Theorem 1. 
([16], Theorem 2), ([17], Theorem 2) Consider a size mesh h > 0 such that N h = τ , for some integer N 1 . Write t n = n h and x n = x ( t n ) , for n N . Then the numerical solution given by x n = f ( t n ) , for N n 0 , and by
x n + 1 = e α h k = 0 m 1 β k h k k ! x n k N + β m ( m 1 ) ! t n m τ t n m τ + h ( t n m τ + h s ) m 1 e α ( t n m τ + h s ) f ( s ) d s ,
where ( m 1 ) τ n h < m τ and m 1 , defines an exact numerical scheme for (1).
Having an exact numerical scheme is ideal, since it reproduces the exact values of the solution at the points of the mesh. However, a drawback of (7) is that definite integrals need to be numerically computed. The number of definite integrals increases with increasing times. Thus, a NSFD method is proposed to maintain sufficient accuracy and adequate dynamical properties, but reduce the complexity by avoiding definite integrals. In the first M intervals [ 0 , τ ] , , [ ( M 1 ) τ , M τ ] , the exact solution (7) is used (or any other numerical method with sufficiently high accuracy), but afterward the integral part from (7) is discarded. The precision of the method increases with M.
Theorem 2. 
([16], Theorem 3), ([17], Theorem 3) Fix M 1 , and compute the numerical solution to (1) in the intervals ( m 1 ) τ n h m τ , for 0 m M with the exact method (7) or with any other numerical method of global error at most O ( h M ) . Then, for m M + 1 and ( m 1 ) τ n h < m τ , the expression
x n + 1 = e α h k = 0 M β k h k k ! x n k N
defines a NSFD scheme of global error O ( h M ) .
As detailed in ([16], Remark 1), the method (8) has the characteristics of a NSFD method:
x n + 1 x n ( e α h 1 ) / α = α x n + α e α h e α h 1 k = 1 M β k h k k ! x n k N .
Furthermore, in the rest of [16], it is proved and illustrated that the method from Theorem 2 is dynamically consistent with (1), for asymptotic stability, positive preserving properties, and oscillation behavior.

4. NSFD Methods for Linear Random Differential Equations with Delay: Approximations of Moments

When α and β are random variables and f is a stochastic process, problem (1) is randomized. These inputs depend on each outcome ω Ω and (2) is obtained. The numerical schemes from Section 3 are also randomized. The exact scheme (7) becomes
x n + 1 ( ω ) = e α ( ω ) h k = 0 m 1 β ( ω ) k h k k ! x n k N ( ω ) + β ( ω ) m ( m 1 ) ! t n m τ t n m τ + h ( t n m τ + h s ) m 1 e α ( ω ) ( t n m τ + h s ) f ( s , ω ) d s ,
where the integral is considered in the m.s. sense (it is assumed that f is a m.s. integrable stochastic process), while (8) becomes
x n + 1 ( ω ) = e α ( ω ) h k = 0 M β ( ω ) k h k k ! x n k N ( ω ) .
We translate Theorem 2 into m.s. convergence.
Theorem 3.
Suppose that α and β are bounded random variables, and that f is a m.s. continuous stochastic process on [ τ , 0 ] . Fix M 1 , and compute the numerical stochastic solution to (2) in the intervals ( m 1 ) τ n h m τ , for 0 m M with the exact method (10). Then, for m M + 1 and ( m 1 ) τ n h < m τ , the expression (11) defines a random NSFD scheme of m.s. global error O ( h M ) .
Proof. 
For n = M N , we have t n = n h = M τ and t n + 1 = ( n + 1 ) h > M τ . For t k , k n , the exact scheme (10) is used, so that x ( t k ) x k 2 = 0 . By (10) and (11), one gets
x ( t n + j + 1 ) x n + j + 1 2 e α h k = 0 M β k h k k ! x ( t n + j k N ) x n + j k N 2 + β m ( m 1 ) ! t n + j m τ t n + j m τ + h ( t n + j m τ + h s ) m 1 e α ( t n + j m τ + h s ) f ( s ) 2 d s ,
for ( m 1 ) τ t n + j < m τ , m 1 M , j 0 .
Let M 0 = β , and M 1 , M 2 > 0 such that e α s < M 1 and f ( s ) 2 < M 2 for s [ 0 , h ] . We first consider j = 0 , , N 1 and m 1 = M . By (12),
x ( t n + 1 ) x n + 1 2 β m ( m 1 ) ! t n m τ t n m τ + h ( t n m τ + h s ) m 1 e α ( t n m τ + h s ) f ( s ) 2 d s M 0 m M 1 M 2 ( m 1 ) ! t n m τ t n m τ + h ( t n m τ + h s ) m 1 d s = M 0 m M 1 M 2 m ! h m C 1 h m = C 1 h M + 1 ,
where C 1 is a constant independent of m and h,
x ( t n + 2 ) x n + 2 2 e α h x ( t n + 1 ) x n + 1 2 + C 1 h M + 1 C 1 e α h + 1 h M + 1 ,
x ( t n + 3 ) x n + 3 2 e α h x ( t n + 2 ) x n + 2 2 + C 1 h M + 1 C 1 e 2 α h + e α h + 1 h M + 1 ,
x ( t n + N ) x n + N 2 C 1 j = 0 N 1 e j α h h M + 1 = C 1 e N α h 1 e α h 1 h M + 1 = C 1 e α τ 1 e α h 1 h M + 1 = C 1 e α τ 1 α h e α h 1 1 α h M C 2 h M ,
where C 2 is a constant that only depends on τ , α and β . Thus,
max n + 1 j n + N x ( t j ) x j 2 = O ( h M ) .
We continue by evaluating x ( t n + N + j ) x n + N + j 2 , 1 j N , by starting from (12) and by employing the bounds already obtained:
x ( t n + N + j ) x n + N + j 2 e α h x ( t n + N + j 1 ) x n + N + j 1 2 + e α h β h C 2 h M + C 1 h M + 1 .
By solving this first order recursive inequality, we derive
x ( t n + N + j ) x n + N + j 2 e j α h x ( t n + N ) x n + N 2 + k = 1 j h M + 1 C 1 + C 2 β e α h e α h ( j k ) C 2 h M e N α h + h M + 1 C 1 + C 2 β e α h e α h j 1 e α h 1 C 2 h M e α τ + h M C 1 + C 2 β e α h e α τ 1 α h e α h 1 1 α C 3 h M ,
where C 3 is a constant that only depends on τ , α and β . Therefore,
max n + N + 1 j n + 2 N x ( t j ) x j 2 = O ( h M ) .
In general, we proceed by induction. Suppose that x ( t j ) x j 2 C h M for n + 1 j n + l N , where C = C ( τ , α , β ) > 0 is constant and l 1 . We prove that max n + l N + 1 j n + ( l + 1 ) N x ( t j ) x j 2 = O ( h M ) . Fix 1 j N . From (12) and by induction hypothesis,
x ( t n + l N + j ) x n + l N + j 2 e α h x ( t n + l N + j 1 ) x n + l N + j 1 2 + e α h k = 1 M β k h k k ! C h M + C h M + 1 .
By solving this first order recursive inequality, we derive (for h < 1 )
x ( t n + l N + j ) x n + l N + j 2 e j α h x ( t n + l N ) x n + l N 2 + k = 1 j C h M + 1 e α h r = 1 M β r h r 1 r ! + 1 e α h ( j k ) C h M e α τ + C h M + 1 e α + β + 1 e α h j 1 e α h 1 C h M e α τ + e α + β + 1 e α τ 1 h e α h 1 C ˜ h M ,
where C ˜ = C ˜ ( τ , α , β ) > 0 is constant. This concludes the proof by induction. □
Remark 1.
As shall be seen in the numerical computations from Section 5, the boundedness of α and β from Theorem 3 is sufficient, but not necessary. Nonetheless, for unbounded α and/or β, if one wants to ensure the m.s. convergence of the NSFD scheme a priori, it is possible to properly truncate the support of α and β. Indeed, since lim m P [ α ( m , m ) ] = 1 , one may take a sufficiently big interval ( m * , m * ) in such a way that P [ α ( m * , m * ) ] 1 , and truncate the support of α to ( m * , m * ) (analogously for β). In fact, by applying the generalized Markov’s inequality, it may be demonstrated that any second order random variable can be truncated to the interval [mean±10×deviation], so that this interval contains 99% of the probability mass irrespective of the probability distribution. In the theory of m.s. calculus, the boundedness of the random input coefficient must be usually imposed: as proved in ([36], Example p. 541), in order for an autonomous and homogeneous first order linear random differential equation (i.e., x ( t ) = a x ( t ) , where a is a random variable) to possess a m.s. solution for every m.s. integrable initial condition x ( 0 ) = x 0 , the coefficient a must be bounded.
M.s. convergence guarantees convergence of the expectation and the variance. In the computer, x n ( ω ) is explicitly and symbolically expressed in terms of α ( ω ) , β ( ω ) and f ( · , ω ) , by employing either the exact scheme (10) along all the integration region or, at cheaper cost, the NSFD scheme from Theorem 3. By using the linearity of the expectation E , one can explicitly compute E [ x n ] . If the NSFD scheme is being used, one is approximating the true expectation of the solution to (2). Complexity is severely increased for large times t, large N (small h), and moderate or high dimension of the random space. The variance may also be approximated by symbolically expressing x n ( ω ) 2 and explicitly computing V [ x n ] = E [ x n 2 ] ( E [ x n ] ) 2 , although the complexity becomes significantly affected because the symbolic expression handled is larger.
Notice that working with these symbolic expressions for x n ( ω ) seems to be necessary. Indeed, if one applies the expectation operator directly in (11), for instance, then E [ x n + 1 ] = k = 0 M E [ e α h β k h k k ! x n k N ] . Since each x n k N depends on α and β , the expectation E [ e α h β k h k k ! x n k N ] cannot be split as E [ e α h ] E [ β k ] h k k ! E [ x n k N ] , unless both α and β are nonrandom. So there does not seem to exist a recursive relation formula for { E [ x n ] } n .
Notice that by Jensen’s and Cauchy–Schwarz inequalities,
| E [ x n ] E [ x ( t n ) ] | = | E [ x n x ( t n ) ] | E [ | x n x ( t n ) | ] x n x ( t n ) 2 .
By triangular, Jensen’s and Cauchy–Schwarz inequalities,
| V [ x n ] V [ x ( t n ) ] | = | E [ ( x n ) 2 ] ( E [ x n ] ) 2 E [ ( x ( t n ) ) 2 ] + ( E [ x ( t n ) ] ) 2 | E [ | ( x n ) 2 ( x ( t n ) ) 2 | ] + | ( E [ x n ] ) 2 ( E [ x ( t n ) ] ) 2 | = E [ | x n x ( t n ) | | x n + x ( t n ) | ] + | E [ x n ] E [ x ( t n ) ] | | E [ x n ] + E [ x ( t n ) ] | x n x ( t n ) 2 x n + x ( t n ) 2 + | E [ x n ] E [ x ( t n ) ] | ( | E [ x n ] | + | E [ x ( t n ) ] | ) x n x ( t n ) 2 ( x n 2 + x ( t n ) 2 ) + | E [ x n ] E [ x ( t n ) ] | ( | E [ x n ] | + | E [ x ( t n ) ] | ) .
So the approximations of the expectation and the variance inherit the rate of convergence corresponding to the m.s. norm, which is O ( h M ) when the exact numerical scheme (10) is used for the first M intervals of length τ and (11) is used for the subsequent intervals (Theorem 3).
In the following section, the m.s. convergence of the NSFD scheme is illustrated with some numerical computations. We point out that the boundedness of α and β from Theorem 3 is sufficient, but not necessary. It may be possible that a random coefficient is unbounded and the NSFD scheme converges in the m.s. sense.

5. Numerical Examples

The theoretical discussion is illustrated with some numerical computations. We consider specific probability distributions for α , β and/or f, and a fixed delay τ > 0 . We denote by x n the discretization of the random NSFD method from Theorem 3. The exact solution x ( t n ) is computed with the exact scheme (10). Both random variables are explicitly and symbolically expressed in terms of α ( ω ) , β ( ω ) and f ( · , ω ) . These expansions are employed to compute the expectation and the variance, by using the linearity of the expectation. To check the accuracy, the absolute errors in the approximations of the mean value, ϵ N , M = | E [ x ( t n ) ] E [ x n ] | , and the variance, δ N , M = | V [ x ( t n ) ] V [ x n ] | , are calculated for different values of N and M. According to the theoretical discussion, the errors should decay as O ( h M ) when h 0 , which entails accuracy up to a significant number of digits. We remark that such level of accuracy cannot be achieved by Monte Carlo simulation, since its error decreases as the reciprocal of the square root of the number of realizations.
The implementations and computations are performed with Mathematica® (Wolfram Research, Inc, Mathematica, Version 12.0, Champaign, IL, USA, 2019), owing to its capability to handle both symbolic and numeric computations.
Example 1.
Let τ = 0.35 . Consider f ( t ) = 1 and α = 1 , while β is a random variable, uniformly distributed on the interval [ 0.1 , 0.2 ] .
Figure 1 plots absolute errors ϵ N , M of the approximation of the expectation. First, N = 10 is fixed and M { 1 , 2 , 3 , 4 } varies. Second, M = 1 is fixed and N { 5 , 7 , 10 } varies. In addition, third, these errors are divided by h to show, because of the overlapping, the decrease O ( h M ) as h 0 . Observe that the error is exactly 0 on the first M intervals of length τ. In the fourth panel, E [ x n ] is plotted, where x n is the output of the NSFD scheme of Theorem 3; the estimated expected values are validated by Monte Carlo simulation. Figure 2 is analogous with the variance and the absolute error of its approximation, δ N , M . According to ([16], Lemma 4, Theorem 7), since α + β < 0 and α β almost surely, the NSFD scheme converges to 0 almost surely as t (it is asymptotically stable almost surely). In both figures, observe that E [ x n ] and V [ x n ] tend to 0 as t , which means that the NSFD scheme is asymptotically stable in the m.s. sense. Finally, the numerical solution is always positive because β > 0 almost surely and f ( t ) > 0 ([16], Theorem 8).
Example 2.
Let τ = 0.35 . Consider f ( t ) = 1 , α = 0 , and β random with Gaussian distribution, of zero mean and 0.3 standard deviation. Notice that the support of β is unbounded; however, we will see that m.s. convergence of the NSFD scheme described in Theorem 3 holds.
Figure 3 reports absolute errors ϵ N , M of the approximation of the expectation, where N and M take on the same values as in Example 1. The decay O ( h M ) as h 0 is captured again. The last panel of the figure plots the expectation of the numerical solution from the NSFD scheme of Theorem 3, E [ x n ] , together with Monte Carlo simulation. Figure 4 is analogous with the variance and the absolute error of its approximation, δ N , M . This example is interesting from the dynamics viewpoint. By ([16], Lemma 4), the probability that the zero solution to the realizations of (2) is asymptotically stable is the probability that β < 0 and τ < τ * = 1 / | β | , i.e., 1 / τ < β < 0 . Taking into account the Gaussian distribution of β, this probability is 0.5 up to 12 decimals, i.e., approximately half of the time a realizable NSFD scheme tends to 0 as t , and half of the time it does not. The m.s. treatment mixes these two behaviors. In the figures, both E [ x n ] and V [ x n ] seem to increase as t advances, which means that the NSFD scheme is unstable in the m.s. sense. Finally, notice that β has one half of probability of being negative and the mean of the solution is positive ([16], Theorem 8).
Example 3.
Let τ = 0.35 . Consider f ( t ) = γ , where γ is a random variable. It is assumed that α, β and γ are independent random quantities, uniformly distributed on the interval [ 0.1 , 0.2 ] .
As Example 1, Figure 5 reports absolute errors ϵ N , M of the approximation of the expectation. First, for fixed N = 10 and M { 1 , 2 , 3 , 4 } . Second, for fixed M = 2 and N { 5 , 7 , 10 } . In addition, third, these errors are divided by h 2 to highlight, because of the overlapping, the decrease O ( h M ) as h 0 . The fourth panel plots E [ x n ] with the discretization x n computed via the NSFD scheme of Theorem 3, which is validated by Monte Carlo simulation. For the variance, computations become more expensive, due to the dimension three of the random space. In particular, the symbolic expression of the exact scheme (10) becomes unfeasible, so the exact error δ N , M of the variance approximation cannot be reported. In Figure 6, we plot V [ x n ] with the discretization x n from the NSFD scheme of Theorem 3. Comparison is performed with Monte Carlo simulation, showing agreement of the estimates. Based on ([16], Lemma 4, Theorem 7), the condition α + β > 0 almost surely entails that the NSFD scheme does not approach 0 as t (almost sure instability). This fact agrees with the plots of E [ x n ] and V [ x n ] , which seem to increase as t grows; this behavior entails that the NSFD scheme is unstable in the m.s. sense. Finally, β > 0 and γ > 0 almost surely implies the positivity of the numerical solution ([16], Theorem 8).
We would like to remark that even when M = 1 and the global error of the NSFD scheme is O ( h ) , its error is lower than Euler’s method, given by x n + 1 = ( 1 + α h ) x n + β h x n N . Euler’s method has already been employed for random differential equations (ordinary and fractional) in the m.s. sense [29,30,34]. In this Example 3, Figure 7 plots errors ϵ N for approximations of the mean, for comparing Euler’s method and the NSFD scheme with M = 1 fixed. Observe that in log scale, both errors are located in parallel, but the error corresponding to the NSFD scheme is lower; this may be due to the non-standard nature of the method and being error-free on [ 0 , τ ] . Although the proposed NSFD method is restricted to the linear random differential equation with delay, it may provide the foundation for designing new non-standard numerical methods for delay nonlinear equations.
Remark 2.  
In the recent literature, the m.s. convergence of Euler’s method has not been formally proved for delay random differential equations. If randomness is not incorporated into the system by the coefficients, but by a Wiener noise instead (Itô stochastic delay differential equation, which gives rise to non-differentiable solutions), then Euler’s method (Euler-Maruyama’s method, which considers discrete increments of the driving Wiener process) was rigorously studied and its m.s. convergence was proved in [37]. In the context of delay random differential equations (those with randomness manifested in coefficients and no Wiener noise), we focus on the linear case (2) studied in this paper assuming, as in Theorem 3, that α and β are bounded random variables, so that α and β are finite. The exact m.s. solution to (2) satisfies x ( t n + 1 ) = x ( t n ) + α t n t n + h x ( s ) d s + β t n t n + h x ( s τ ) d s , where the integrals are m.s. Riemann. If e n = x ( t n ) x n denotes the difference between the exact solution and the discretization at the mesh point t n , then
e n + 1 = e n + α t n t n + h ( x ( s ) x n ) d s + β t n t n + h ( x ( s τ ) x n N ) d s ,
by a simple subtraction. Given any interval [ τ , T ] , T > 0 , the m.s. Lipschitz condition x ( s 1 ) x ( s 2 ) 2 = s 2 s 1 x ( s ) d s 2 | s 2 s 1 x ( s ) 2 d s | λ | s 1 s 2 | holds, λ = λ ( T ) = max [ τ , T ] x ( s ) 2 > 0 . Then x ( s ) x m 2 x ( s ) x ( t m ) 2 + e m 2 λ h + e m 2 , m N , s [ t m , t m + h ] , t m + h T . Consequently,
e n + 1 2 e n 2 + α h λ h + e n 2 + β h λ h + e n N 2 = 1 + α h e n 2 + β h e n N 2 + α + β h 2 λ .
For 0 n N (notice that e n N = 0 ), by solving the first order recursive inequality for e n 2 , we derive the following bounds:
e n 2 i = 1 n α + β h 2 λ 1 + α h n i = h 2 λ α + β 1 + α h n 1 α h h λ α + β α 1 + α τ / N N 1 h λ α + β α e α τ 1 = C 1 h ,
where C 1 = C 1 ( λ , τ , α , β ) is constant. Then
max 0 n N e n 2 = O ( h ) .
For N + 1 n 2 N , based on similar calculations,
e n 2 1 + α h n N e N 2 + i = N + 1 n β e i N 2 + α + β λ h h 1 + α h n i C 1 h e α τ + C 1 β + α + β λ h 2 1 + α h n N 1 α h C 2 h ,
where C 2 = C 2 ( λ , τ , α , β ) is constant. Then
max N + 1 n 2 N e n 2 = O ( h ) .
For 2 N + 1 n 3 N , 3 N + 1 n 4 N , etc. one proceeds similarly. This proves that the random Euler’s method has m.s. global error O ( h ) . This proof corresponds to the linear case, although it may be extendible to delay random differential equations satisfying a m.s. Lipschitz condition.
Example 4.
In this example, the linear model (2) is considered for fitting the time evolution of a photosynthetic bacterial population, Rhodobacter capsulatus (R. capsulatus) [38], under infrared lighting conditions. Direct cell counts were made for the first 7 days, every two to three days, during which the population grew with no effect of competition for resources (light and/or CO2) that would yield logistic nonlinearities. For days 0, 2, 4 and 7, the population sizes, measured in cells/mL scaled by one million, were 0.583, 0.635, 1.08 and 3.20, respectively. For delay τ = 1 and initial function f ( t ) = 0.583 on [ τ , 0 ] , the least-squares estimates for α and β are 1.20426 and 1.18024 , respectively. The effect of small random displacements on the coefficients is studied here. Let us suppose 0.5 % displacements of α and β with respect to their least-squares estimates, with zero mean values. According to the maximum entropy principle [39], α and β follow truncated exponential distributions, with rates 1 / 1.20426 and 1 / 1.18024 respectively. The expectation and the variance of the output are approximated with the random NSFD scheme. Figure 8 plots the results for M = 2 and distinct N (mean values in solid line, and mean ± 2 × standard deviation in dashed lines), together with the least-squares fitting. Observe that a small uncertainty of 0.5 % for parameters may cause significant changes in the final solution, up to 30% variation for the seventh day compared to an idealized situation containing no uncertainty. Observe also that as N increases, the approximations from the NSFD scheme tend to overlap, thus indicating convergence.

6. Conclusions

In this paper, we have extended a NSFD numerical scheme recently proposed for deterministic linear differential equations with delay to the random framework. Incorporating randomness into models is important to account for measurement errors in data. M.s. convergence of the numerical discretizations has been established when the two equation coefficients are bounded random variables and the initial condition is a regular stochastic process, with rate of convergence given by O ( h M ) , where h is the step size and M is the number of intervals of length τ where the exact scheme is applied. M.s. convergence allows for approximating the expectation and the variance of the solution at inherited rate O ( h M ) , by symbolically expanding the discretizations in terms of the random inputs. The numerical examples have illustrated and assessed the proposed approach. The convergence rate O ( h M ) has been supported numerically, even when the two equation coefficients are unbounded random variables. The asymptotic behavior of the expectation and the variance as the time t grows has been evaluated, graphically and taking into account theoretical results on deterministic stability and instability of the zero solution. A comparison with Euler’s method has been performed when M = 1 ; although both methods have global errors O ( h ) , the error of the NSFD scheme is lower, possibly due to the non-standard nature of the scheme and being error-free on [ 0 , τ ] . Also, we have considered an example dealing with actual experimental data for a bacterial specie growing under infrared lighting conditions, and have calculated numerical solutions after randomizing the input parameters according to the maximum entropy principle.
The advantage of the random NSFD scheme is the high accuracy to approximate some statistics, which cannot be achieved with Monte Carlo methods. In addition, the procedure is simple: one only symbolically expands the discretizations in terms of the random inputs and afterward applies the corresponding statistical operator. However, this strategy possesses some limitations. Obviously, the necessity of symbolically expressing the discretizations restricts the applicability of the NSFD scheme to moderate step size h and time variable t, as well as small dimension of the random space. Although the calculation of the expectation of the discretization seems to be quite feasible in the computer, the calculation of the variance may become a big issue with this approach, let alone other statistics of order greater than two. We ask ourselves about the possibility of accurately calculating statistics with the random NSFD scheme without relying on symbolic expansions. We admit that for the moment, Monte Carlo simulation seems the best option for large time variable t, small step size h or large dimension of the random space, where each realization of the governing delay model is numerically solved by employing a NSFD scheme. For estimating densities, the symbolic expression is too complex and kernel methods are the preferable.
Mickens’ methodology on NSFD schemes has shown fruitful applications along the years on ordinary, partial and fractional deterministic differential equations. However, only very recently, the application of NSFD numerical schemes to delay deterministic differential equations has been explored, in the context of linear models. Thus, this is the first contribution that proposes the use of NSFD schemes for quantifying uncertainty for delay random differential equations. Further study of NSFD methods for delay deterministic and random differential equations needs to be conducted, especially for nonlinear equations, for applications to modeling of real-life systems with aftereffects or time lags.
We propose specific lines of research for possible future developments:
  • In the deterministic setting, the extension of the NSFD method to nonlinear delay differential equations. We believe that this extension may be done by linearization or by applying the empirical rules proposed by Mickens.
  • Randomization of the NSFD method for delay random differential equations and applications without relying on symbolic expansions. Symbolic computations are the main drawback of the method proposed in the present paper.
  • A theoretical analysis of m.s. dynamic consistency.

Author Contributions

Investigation, J.C. and M.J.; Methodology, J.C. and M.J.; Software, M.J.; Supervision, J.C.C. and F.R.; Validation, all authors; Visualization, all authors; Writing—original draft, J.C. and M.J.; Writing—review & editing, all authors. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been supported by the Spanish Ministerio de Economía, Industria y Competitividad (MINECO), the Agencia Estatal de Investigación (AEI) and Fondo Europeo de Desarrollo Regional (FEDER UE) grant MTM2017–89664–P.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this article.

References

  1. Smith, H. An Introduction to Delay Differential Equations with Applications to the Life Sciences; Texts in Applied Mathematics; Springer: New York, NY, USA, 2011. [Google Scholar]
  2. Driver, Y. Ordinary and Delay Differential Equations; Applied Mathematical Science Series; Springer: New York, NY, USA, 1977. [Google Scholar]
  3. Diekmann, O.; van Cils, S.A.; Verduyn Lunel, S.M.; Walther, H.-O. Delay Equations; Springer: New York, NY, USA, 1995. [Google Scholar]
  4. Saaty, T.L. Modern Nonlinear Equations; Dover Publications: New York, NY, USA, 1981. [Google Scholar]
  5. Bocharov, G.A.; Rihan, F.A. Numerical modelling in biosciences using delay differential equations. J. Comput. Appl. Math. 2000, 125, 183–199. [Google Scholar] [CrossRef] [Green Version]
  6. Jackson, M.; Chen-Charpentier, B.M. Modeling plant virus propagation with delays. J. Comput. Appl. Math. 2017, 309, 611–621. [Google Scholar] [CrossRef]
  7. Chen-Charpentier, B.M.; Diakite, I. A mathematical model of bone remodeling with delays. J. Comput. Appl. Math. 2016, 291, 76–84. [Google Scholar] [CrossRef]
  8. Erneux, T. Applied Delay Differential Equations; Surveys and Tutorials in the Applied Mathematical Sciences Series; Springer: New York, NY, USA, 2009. [Google Scholar]
  9. Kyrychko, Y.N.; Hogan, S.J. On the use of delay equations in engineering applications. J. Vib. Control 2017, 16, 943–960. [Google Scholar] [CrossRef]
  10. Harding, L.; Neamtu, M. A dynamic model of unemployment with migration and delayed policy intervention. Comput. Econ. 2018, 51, 427–462. [Google Scholar] [CrossRef] [Green Version]
  11. Mickens, R.E. Nonstandard Finite Difference Models of Differential Equations; World Scientific: Singapore, 1994. [Google Scholar]
  12. Mickens, R.E. Applications of Nonstandard Finite Difference Schemes; World Scientific: Singapore, 2000. [Google Scholar]
  13. Mickens, R.E. Advances on Applications of Nonstandard Finite Difference Schemes; World Scientific: Singapore, 2005. [Google Scholar]
  14. Mickens, R.E. Dynamic consistency: A fundamental principle for constructing nonstandard finite difference schemes for differential equations. J. Differ. Equ. Appl. 2005, 11, 645–653. [Google Scholar] [CrossRef]
  15. Patidar, K.C. Nonstandard finite difference methods: Recent trends and further developments. J. Differ. Equ. Appl. 2016, 22, 817–849. [Google Scholar] [CrossRef]
  16. García, M.A.; Castro, M.A.; Martín, J.A.; Rodríguez, F. Exact and nonstandard numerical schemes for linear delay differential models. Appl. Math. Comput. 2018, 338, 337–345. [Google Scholar] [CrossRef]
  17. Castro, M.Á.; García, M.A.; Martín, J.A.; Rodríguez, F. Exact and Nonstandard Finite Difference Schemes for Coupled Linear Delay Differential Systems. Mathematics 2019, 7, 1038. [Google Scholar] [CrossRef] [Green Version]
  18. Soong, T.T. Random Differential Equations in Science and Engineering; Academic Press: New York, NY, USA, 1973. [Google Scholar]
  19. Villafuerte, L.; Braumann, C.A.; Cortés, J.-C.; Jódar, L. Random differential operational calculus: Theory and applications. Comput. Math. Appl. 2010, 59, 115–125. [Google Scholar] [CrossRef] [Green Version]
  20. Cortés, J.-C.; Jódar, L.; Roselló, M.-D.; Villafuerte, L. Solving initial and two-point boundary value linear random differential equations: A mean square approach. Appl. Math. Comput. 2012, 219, 2204–2211. [Google Scholar] [CrossRef]
  21. Calatayud, J.; Cortés, J.-C.; Jornet, M.; Villafuerte, L. Random non-autonomous second order linear differential equations: Mean square analytic solutions and their statistical properties. Adv. Differ. Equ. 2018, 392, 1–29. [Google Scholar] [CrossRef]
  22. Calatayud, J.; Cortés, J.-C.; Jornet, M. Improving the approximation of the first- and second-order statistics of the response stochastic process to the random Legendre differential equation. Mediterr. J. Math. 2019, 16, 68. [Google Scholar] [CrossRef]
  23. Licea, J.A.; Villafuerte, L.; Chen-Charpentier, B.M. Analytic and numerical solutions of a Riccati differential equation with random coefficients. J. Comput. Appl. Math. 2013, 239, 208–219. [Google Scholar] [CrossRef]
  24. Burgos, C.; Calatayud, J.; Cortés, J.-C.; Villafuerte, L. Solving a class of random non-autonomous linear fractional differential equations by means of a generalized mean square convergent power series. Appl. Math. Lett. 2018, 78, 95–104. [Google Scholar] [CrossRef] [Green Version]
  25. Calatayud, J.; Cortés, J.-C.; Jornet, M. Random differential equations with discrete delay. Stoch. Anal. Appl. 2019, 37, 699–707. [Google Scholar] [CrossRef]
  26. Calatayud, J.; Cortés, J.-C.; Jornet, M. Lp-calculus approach to the random autonomous linear differential equation with discrete delay. Mediterr. J. Math. 2019, 16, 85. [Google Scholar] [CrossRef]
  27. Cortés, J.-C.; Jornet, M. Lp-solution to the random linear delay differential equation with a stochastic forcing term. Mathematics 2020, 8, 1013. [Google Scholar] [CrossRef]
  28. Caraballo, T.; Cortés, J.-C.; Navarro-Quiles, A. Applying the Random Variable Transformation method to solve a class of random linear differential equation with discrete delay. Appl. Math. Comput. 2019, 356, 198–218. [Google Scholar] [CrossRef]
  29. Cortés, J.-C.; Jódar, L.; Villafuerte, L. Numerical solution of random differential equations: A mean square approach. Math. Comput. Model. 2007, 45, 757–765. [Google Scholar] [CrossRef]
  30. Cortés, J.-C.; Jódar, L.; Villafuerte, L. Mean square numerical solution of random differential equations: Facts and possibilities. Comput. Math. Appl. 2007, 53, 1098–1106. [Google Scholar] [CrossRef] [Green Version]
  31. El-Tawil, M.A. The approximate solutions of some stochastic differential equations using transformations. Appl. Math. Comput. 2005, 164, 167–178. [Google Scholar] [CrossRef]
  32. Calatayud, J.; Cortés, J.-C.; Díaz, J.A.; Jornet, M. Density function of random differential equations via finite difference schemes: A theoretical analysis of a random diffusion-reaction Poisson-type problem. Stochastics 2020, 92, 627–641. [Google Scholar] [CrossRef]
  33. Calatayud, J.; Cortés, J.-C.; Díaz, J.A.; Jornet, M. Constructing reliable approximations of the probability density function to the random heat PDE via a finite difference scheme. Appl. Numer. Math. 2020, 151, 413–424. [Google Scholar] [CrossRef]
  34. Burgos, C.; Cortés, J.-C.; Villafuerte, L.; Villanueva, R.-J. Mean square convergent numerical solutions of random fractional differential equations: Approximations of moments and density. J. Comput. Appl. Math. 2020, 378, 112925. [Google Scholar] [CrossRef]
  35. Loève, M. Probability Theory; Springer: New York, NY, USA, 1977; Volume I–II. [Google Scholar]
  36. Strand, J.L. Random ordinary differential equations. J. Differ. Equ. 1970, 7, 538–553. [Google Scholar] [CrossRef] [Green Version]
  37. Buckwar, E. Introduction to the numerical analysis of stochastic delay differential equations. J. Comp. Appl. Math. 2000, 125, 297–307. [Google Scholar] [CrossRef] [Green Version]
  38. Stanescu, D.; Chen-Charpentier, B.; Jensen, B.J.; Colberg, P.J. Random coefficient differential models of growth of anaerobic photosynthetic bacteria. Electron. Trans. Numer. Anal. 2009, 34, 44–58. [Google Scholar]
  39. Dorini, F.A.; Sampaio, R. Some results on the random wear coefficient of the Archard model. J. Appl. Mech. 2012, 79, 051008–051014. [Google Scholar] [CrossRef]
Figure 1. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 1.
Figure 1. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 1.
Mathematics 08 01417 g001
Figure 2. Upper left panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with N = 5 , 7 , 10 , 20 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the variance with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 1.
Figure 2. Upper left panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with N = 5 , 7 , 10 , 20 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the variance with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 1.
Mathematics 08 01417 g002
Figure 3. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 2.
Figure 3. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 2.
Mathematics 08 01417 g003
Figure 4. Upper left panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the variance with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 2.
Figure 4. Upper left panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right panel: Absolute errors (log-scale) in the approximation of the variance with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 1 . Lower left panel: Errors from the upper right panel divided by h. Lower right panel: Approximation of the variance with the NSFD scheme, with M = 1 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 2.
Mathematics 08 01417 g004
Figure 5. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 2 . Lower left panel: Errors from the upper right panel divided by h 2 . Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 2 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 3.
Figure 5. Upper left panel: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with M = 1 , 2 , 3 , 4 and N = 10 ( h = τ / N ). Upper right: Absolute errors (log-scale) in the approximation of the mean value with the NSFD scheme, with N = 5 , 7 , 10 ( h = τ / N ) and M = 2 . Lower left panel: Errors from the upper right panel divided by h 2 . Lower right panel: Approximation of the expectation with the NSFD scheme, with M = 2 and N = 7 , and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 3.
Mathematics 08 01417 g005
Figure 6. Approximation of the variance with the NSFD scheme, with M = 1 and N = 5 ( h = τ / N ), and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 3.
Figure 6. Approximation of the variance with the NSFD scheme, with M = 1 and N = 5 ( h = τ / N ), and comparison with Monte Carlo simulation (circles) using 10,000 realizations. This figure corresponds to Example 3.
Mathematics 08 01417 g006
Figure 7. Absolute errors (log-scale) in the approximation of the mean value, with Euler’s method (dashed lines) and with the NSFD scheme M = 1 (solid lines). For step size h = τ / N , the left panel corresponds to N = 7 and the right panel to N = 10 . This figure corresponds to Example 3.
Figure 7. Absolute errors (log-scale) in the approximation of the mean value, with Euler’s method (dashed lines) and with the NSFD scheme M = 1 (solid lines). For step size h = τ / N , the left panel corresponds to N = 7 and the right panel to N = 10 . This figure corresponds to Example 3.
Mathematics 08 01417 g007
Figure 8. Application of the random NSFD scheme for the model of the R. capsulatus bacterial population, for M = 2 and different values of N (mean values in solid line, and mean ± 2 × standard deviation in dashed lines). The least-squares fitting is also plotted. A zoom for a particular region is included for a better appreciation of convergence as N grows. This figure corresponds to Example 4.
Figure 8. Application of the random NSFD scheme for the model of the R. capsulatus bacterial population, for M = 2 and different values of N (mean values in solid line, and mean ± 2 × standard deviation in dashed lines). The least-squares fitting is also plotted. A zoom for a particular region is included for a better appreciation of convergence as N grows. This figure corresponds to Example 4.
Mathematics 08 01417 g008

Share and Cite

MDPI and ACS Style

Calatayud, J.; Cortés, J.C.; Jornet, M.; Rodríguez, F. Mean Square Convergent Non-Standard Numerical Schemes for Linear Random Differential Equations with Delay. Mathematics 2020, 8, 1417. https://doi.org/10.3390/math8091417

AMA Style

Calatayud J, Cortés JC, Jornet M, Rodríguez F. Mean Square Convergent Non-Standard Numerical Schemes for Linear Random Differential Equations with Delay. Mathematics. 2020; 8(9):1417. https://doi.org/10.3390/math8091417

Chicago/Turabian Style

Calatayud, Julia, Juan Carlos Cortés, Marc Jornet, and Francisco Rodríguez. 2020. "Mean Square Convergent Non-Standard Numerical Schemes for Linear Random Differential Equations with Delay" Mathematics 8, no. 9: 1417. https://doi.org/10.3390/math8091417

APA Style

Calatayud, J., Cortés, J. C., Jornet, M., & Rodríguez, F. (2020). Mean Square Convergent Non-Standard Numerical Schemes for Linear Random Differential Equations with Delay. Mathematics, 8(9), 1417. https://doi.org/10.3390/math8091417

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop