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

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.


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

M.s. Calculus
We are interested in second order real random variables y : Ω → R, satisfying 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: 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[xy] ([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[ 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 n 0 < t n 1 < . . . < t n r n = b}, such that for any choice of points s n i ∈ [t n i−1 , t n i ], i = 1, . . . , r n , the limit lim n→∞ ∑ r n i=1 z(s n i )(t n i − t n i−1 ) exists in the m.s. sense, and it is denoted as Finally, we mention that the essential supremum norm is defined as 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 < ∞.

NSFD Methods for Linear Deterministic Differential Equations with Delay
Based on the explicit solution to (1), 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.
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.
As detailed in ( [16], Remark 1), the method (8) has the characteristics of a NSFD method: 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.

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 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 We translate Theorem 2 into m.s. convergence. Proof. For n = MN, we have t n = nh = 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 We first consider j = 0, . . . , N − 1 and m − 1 = M. By (12), where C 1 is a constant independent of m and h, . . .
where C 2 is a constant that only depends on τ, α ∞ and β ∞ . Thus, 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: By solving this first order recursive inequality, we derive where C 3 is a constant that only depends on τ, α ∞ and β ∞ . Therefore, In general, we proceed by induction. Suppose that (12) and by induction hypothesis, By solving this first order recursive inequality, we derive (for h < 1) 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) = ax(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 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) By triangular, Jensen's and Cauchy-Schwarz inequalities, 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.

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.   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 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).  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 + βhx 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 ) + α Consequently, e n+1 2 ≤ e n 2 + α ∞ h (λh + e n 2 ) + β ∞ h (λh + e n−N 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: For N + 1 ≤ n ≤ 2N, based on similar calculations, For 2N + 1 ≤ n ≤ 3N, 3N + 1 ≤ n ≤ 4N, 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 CO 2 ) 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.

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