Applying the random variable transformation method to solve a class of random linear differential equation with discrete delay

Applying the Random Variable Transformation method to solve a class of random linear differential equation with discrete delay Tom´ as Caraballoa, J.-C. Cort´ esb, A. Navarro-Quilesc,d,∗ aDpto. Ecuaciones Diferenciales y An´alisis Num´erico, Universidad de Sevilla, c/Tarfia s/n, 41012, Sevilla, Spain bInstituto Universitario de Matem´atica Multidisciplinar, Universitat Polit`ecnica de Val`encia, c/Camino de Vera s/n, 46022, Valencia, Spain cDeustoTech, University of Deusto, 48007 Bilbao, Basque Country, Spain dFacultad de Ingenieria, Universidad de Deusto, Avda.Universidades, 24, 48007, Bilbao, Basque Country, Spain. Abstract We randomize the following class of linear differential equations with delay, x0 τ(t)=axτ(t)+ bxτ(t−τ), t>0, and initial condition, xτ(t)=g(t), −τ≤t≤0, by assuming that coefficients aand bare random variables and the initial condition g(t) is a stochastic process. We consider two cases, depending on the functional form of the stochastic process g(t), and then we solve, from a probabilistic point of view, both random initial value problems by determining explicit expressions to the first probability density function, f(x,t;τ), of the corresponding solution stochastic processes. Afterwards, we establish sufficient conditions on the involved random input parameters in order to guarantee that f(x,t;τ) converges, as τ→0+, to the first probability density function, say f(x,t), of the corresponding associated random linear problem without delay (τ=0). The paper concludes with several numerical experiments illustrating our theoretical findings. Keywords: Random linear differential equation with delay, Probability density function, Random Variable Transformation technique. 1. Introduction and motivation1 Ordinary differential equations are useful mathematical tools to model phenomena in areas2 like Physics, Engineering, Epidemiology, Economics, etc. In many applications, differential3 equations are formulated using the principle of causality based upon the fact that future state of a4 system under study is independent of its past state and is solely determined by the current state.5 Although many phenomena can be properly described using this tenet, there are other situations6 ∗Corresponding author Email addresses: [email protected] (Tom´ as Caraballo), [email protected] (J.-C. Cort´ es), [email protected] (A. Navarro-Quiles ) Preprint submitted to Applied Mathematics and Computation (Elsevier) March 7, 2019 where it may be more realistic to model the current state of a physical system (understood in7 a wide sense) in terms of past information. In that case, it is more suitable to describe the8 dynamics of the physical system by means of differential equations incorporating the past history9 of the system under analysis. These kinds of differential equations are usually referred to as10 delay differential equations (DDEs). One commonly distinguishes two main classes of DDEs11 depending on the type of delay considered therein. If only a part of the history has a relevant12 influence on the current state, then discrete DDEs are formulated, while continuous DDEs are13 those whose delay is unbounded or infinite. In this latter case, the whole past history is taken into14 account to describe the phenomenon under study. In this paper, we will deal with the following15 class of linear discrete DDEs with initial condition16 (x0 τ(t)=axτ(t)+bxτ(t−τ),t>0, τ > 0, xτ(t)=g(t),−τ≤t≤0,(1) where τ > 0 denotes a prefixed finite delay, ais the coefficient of the non-delay term, xτ(t), b17 is the coefficient of the delay term, xτ(t−τ), and g(t) is an arbitrary function (initial condition)18 defined on the interval [−τ, 0].To avoid confusion with the notation introduced for the initial19 condition and the solution of the initial value problem (IVP) formulated in (1), hereinafter the20 solution will be denoted by xτ(·). This notation is necessary to explicitly indicate the dependence21 of the solution on the τparameter since later we will study the convergence of the solution of the22 IVP (1) as τ→0+to the solution of the corresponding IVP without delay. Now, if we assume23 that g(t) is continuous in [−τ, 0] and differentiable in ]−τ, 0[ then, according to Theorem 1 below,24 there exists an exact expression of the solution of IVP (1).25 Theorem 1 ([1, 2]). Let us consider IVP (1) and assume that g(t)is continuous in [−τ, 0] and26 differentiable in ]−τ, 0[, i.e., g(·)∈ C1([−τ, 0]). Then, IVP (1) has a unique solution xτ(·)∈27 C0([−τ, ∞]) ∩C1([−τ, 0]) ∩C1([0,∞]) given by28 xτ(t)=ea(t+τ)eb1,t τg(−τ)+Z0 −τ ea(t−s)eb1,t−τ−s τ(g0(s)−ag(s)) ds,(2) where b1=e−aτb, and eb1,t τand eb1,t−τ−s τdenote the delayed exponential function, ec,t τ, evaluated29 at (c,t)=(b1,t) and (c,t)=(b1,t−τ−s), respectively. This function appears in a natural way30 in dealing with the linear discrete DDE (1) since its solution is constructed segment by segment31 (see for example [2]). For the sake of completeness, below we recall its definition.32 Definition 1. ([2]) Let c be a real number and τ > 0, then the function33 ec,t τ=                                            0,−∞ <t<−τ, 1,−τ≤t<0, 1+ct 1!,0≤t< τ, 1+ct 1! +c2(t−τ)2 2! , τ ≤t<2τ, . . .. . . n X k=0 ck(t−(k−1)τ)k k!,(n−1)τ≤t<nτ, (3) is called the delayed exponential function, where n =bt/τc+1, being bxcthe greatest integer less34 than or equal to x.35 2 So far we have revised the main definitions and results involving the solution of deterministic36 IVP (1) for the linear discrete DDE. When this class of equations is applied to model the dynam-37 ics of real phenomena, its input parameters, i.e. the coefficients aand b, and the initial condition,38 g(t), must be fixed from experimental data which often involve uncertainties because they are39 obtained after measurements and sampling. This fact allows us to treat the input coefficients,40 aand b, as random variables (RVs), and the initial condition, g(t), as a stochastic process (SP)41 rather than deterministic constants and a classical function, respectively. This leads to the full42 randomization of IVP (1)43 (x0 τ(t;ω)=a(ω)xτ(t;ω)+b(ω)xτ(t−τ;ω),t>0, τ > 0, xτ(t;ω)=g(t;ω),−τ≤t≤0,(4) where a(ω) and b(ω) are assumed to be absolutely continuous RVs and g(t;ω) is a SP, being44 all of them defined on a common complete probability space (Ω,F,P). In accordance with45 Theorem 1, and using the so-called “sample random calculus” for SPs [3, App. I], hereinafter it46 will be assumed that the initial condition SP, g(t;ω), satisfies the following condition,47 g(·;ω)∈ C1([−τ, 0] ×Ω),a.e.,(5) where, as usual, a.e. stands for ‘almost everywhere’, so that the existence of a unique sample48 solution SP, xτ(t;ω), to random IVP (4) can be guaranteed.49 The study of differential equations with delay involving uncertainties has been studied from50 different approaches. In [4], authors study a class of stochastic impulsive differential equations51 involving Bernoulli distribution where trial lengths vary randomly. In [5], the complete controlla-52 bility property of a class of nonlinear stochastic differential equation with delay, in the fractional53 sense, is investigated assuming that delays are described by Poisson jumps. Stochastic differen-54 tial equations with delay have been proposed to model interesting real problems. For example, in55 [6] the nonlinear delay differential neoclassical growth model is analysed assuming that stochas-56 tic perturbations of the white noise type. In [7], authors provide sufficient conditions for stability57 in probability of the equilibrium point of a social obesity epidemic model with distributed delay58 and stochastic perturbations. In dealing with delay differential equations it is usual to study the59 behaviour of the solution when the delay tends to zero, i.e., to investigate conditions under which60 there is convergence of the solution of the IVP with delay to the corresponding solution of the61 IVP when the delay vanishes. In [8] it is shown that the solution of a mixed stochastic delay62 differential equation depends continuously on the coefficients and the initial data. In addition,63 authors prove the convergence of the solutions to equations with vanishing delay to the solution64 of corresponding equations without delay. In [9] one deals with the mean square convergence65 and mean square exponential stability of an Euler scheme for a linear impulsive stochastic delay66 differential equation.67 Solving a random (ordinary/partial/fractional/delay/etc.) differential equation means not just68 to obtain an exact/approximate solution SP, but also its main statistical properties, like the mean69 and the variance functions. These equations are said to be solved, from a probabilistic standpoint,70 when the first probability density function (1-PDF) of the solution SP is exact/approximately71 obtained since from this deterministic function one can completely characterize the probabilistic72 behaviour of the solution SP at every time instant. As a consequence, mean, variance, skewness,73 etc., as well as any one-dimensional moment of the solution SP can be derived from the 1-PDF,74 provided these moments exist. To be specific, if f(x,t;τ) denotes the 1-PDF of the solution SP75 xτ(t;ω) to random IVP (1), then the mean, E[xτ(t;ω)], and the variance, V[xτ(t;ω)], functions76 3 can be obtained by77 Eh(xτ(t;ω))ki=ZR xkf(x,t;τ) dx,k=1,2, . . . as78 E[xτ(t;ω)] =ZR x f (x,t;τ) dx,V[xτ(t;ω)] =ZR x2f(x,t;τ) dx−(E[xτ(t;ω)])2,(6) respectively. Fixed a time instant, say ˆ t, the computation of the probability that the solution,79 xτ(ˆ t;ω), lies within an interval of specific interest can also be computed just by integrating the80 1-PDF81 Pω∈Ω:a≤xτ(ˆ t;ω)≤b=Zb a f(x,ˆ t;τ) dx. At this point is important to stress that, to the best of our knowledge, this approach has already82 been dealt with some classes of random fractional, ordinary and partial differential equations and83 of random difference equations as well (see for instance, [10], [11, 12, 13, 14, 15], [16, 17, 18,84 19] and [20, 21], respectively), but the corresponding analysis for random DDEs has not been85 addressed yet.86 In the spirit of these previous contributions, the main objective of this paper is solving, from87 a probabilistic point of view, random IVP (4) by obtaining the 1-PDF, f(x,t;τ), of its solution88 SP, xτ(t;ω), which, according to the deterministic solution formulated in (2), for each t∈[(n−89 1)τ, nτ[, with τ > 0 fixed and n=1,2, . . ., is given by90 xτ(t;ω)=ea(ω)(t+τ)eb1(ω),t τg(−τ;ω)+Z0 −τ ea(ω)(t−s)eb1(ω),t−τ−s τg0(s;ω)−a(ω)g(s;ω)ds,(7) where b1(ω)=e−a(ω)τb(ω). The key tool that will be applied to achieve this goal is the Random91 Variable Transformation (RVT) method. This technique allows us to obtain the PDF of a random92 vector, which results from mapping another random vector whose PDF is known. The following93 result provides the RVT technique in its multidimensional version.94 Theorem 2 (Multidimensional RVT method, [3]). Let x(ω)=(x1(ω),...,xm(ω)) and y(ω)= (y1(ω),...,ym(ω)) be two m-dimensional absolutely continuous random vectors defined on a complete probability space (Ω,F,P). Let r:Rm→Rmbe a one-to-one deterministic transformation of x(ω)into y(ω), i.e., y(ω)=r(x(ω)). Assume that ris continuous in xand has continuous partial derivatives with respect to x. Then, if fx(x)denotes the joint probability density function of the random vector x(ω), and s=r−1=(s1(y1,...,ym),...,sm(y1,...,ym)) represents the inverse mapping of r=(r1(x1,...,xm),...,rm(x1,...,xm)), the joint probability density function of the random vector y(ω)is given by gy(y)=fx(s(y)) |Jm|, where |Jm|, which is assumed to be different from zero, denotes the absolute value of the Jacobian defined by the determinant Jm=det                      ∂s1(y1,...,ym) ∂y1··· ∂sm(y1,...,ym) ∂y1 . . ..... . . ∂s1(y1,...,ym) ∂ym··· ∂sm(y1,...,ym) ∂ym                      . 4 Once we have obtained the 1-PDF, f(x,t;τ), of the solution SP to random IVP (4), the next95 objective of this paper is to study the relationship between f(x,t;τ) and the 1-PDF, say f(x,t),96 of the corresponding random IVP without delay, i.e.,97 (x0(t;ω)=(a(ω)+b(ω))x(t;ω),t>0, x(0; ω)=g(0; ω)=g0(ω),(8) where g0(ω) is an absolutely continuous RV. To be specific, we will establish conditions in order98 to guarantee that99 lim τ→0+f(x,t;τ)=f(x,t),for each (x,t)∈ D(xτ(t;ω)) ∩D(x(t;ω)) ×[(n−1)τ, nτ[ fixed,(9) being n=1,2, . . . and where D(xτ(t;ω)) and D(x(t;ω)) denote the codomains of SPs xτ(t;ω)100 and x(t;ω), respectively. This analysis will focus on the following choices of the initial condition101 g(t;ω) for random IVP (4)102 •Case I: g(t;ω)=ea(ω)t+c(ω).103 •Case II: g(t;ω)= m X j=0 cj(ω)tj,m≥0.104 The choice in Case I has been made because, as it will be seen later, it allows us to deal with105 ascenario that illustrates adequately the main ideas of our approach. While Case II involves106 the general case in which the initial condition is a random polynomial of degree m, which has107 interest by itself since, similarly to what happens in the deterministic context, many important108 random SPs can be approximated by appropriate random polynomials.109 The rest of this paper is organized as follows. Section 2 is split into two subsections, in110 the first one, we determine an explicit expression to the 1-PDF, f(x,t;τ), of the solution SP to111 random IVP (1) in Case I, i.e., when g(t;ω)=ea(ω)t+c(ω). Subsection 2.1 is addressed to establish112 sufficient conditions upon the involved random input parameters a(ω), b(ω) and c(ω), in order113 to guarantee that f(x,t;τ) converges, as τ→0+, to the 1-PDF, f(x,t), of the corresponding114 associated random linear problem without delay (τ=0). For the sake of clarity, Section 3 is115 organized analogously as Section 2 to conduct the corresponding study in Case II, i.e., when116 the initial condition to random IVP (1) is a random polynomial g(t;ω)=Pm j=0cj(ω)tj,m≥117 0. Several numerical examples, corresponding to Cases I and II are exhibited in Section 4.118 Conclusions are drawn in Section 5.119 2. Case I: Computing the 1-PDF of the solution SP and study of the convergence120 Let us consider random IVP (4) where the initial condition is given by the SP g(t;ω)=121 ea(ω)t+c(ω). Observe that this situation corresponds to the case where g(t;ω) is the solution of the122 following random IVP123 (g0(t;ω)=a(ω)g(t;ω),t≥ −τ, τ > 0, g(−τ;ω)=e−a(ω)τ+c(ω).(10) In other words, we are then implicitly considering a stochastic control problem defined by (4)124 and (10). In agreement with (7), the solution SP of this control problem is given by125 xτ(t;ω)=ea(ω)t+c(ω)eb1(ω),t τ,t∈[(n−1)τ, nτ[,n∈Nfixed,(11) 5 where b1(ω)=e−a(ω)τb(ω). Henceforth, we will assume that the random inputs a(ω), b(ω) and126 c(ω) are absolutely continuous dependent RVs with a joint PDF denoted by fc,a,b(c,a,b), which127 is assumed to be known. Observe that g(t;ω) satisfies condition (5).128 2.1. Computing the PDF129 First at all, notice that if ˆ tis such that xτ(ˆ t;ω)=0, then clearly the 1-PDF is given by130 f(x,ˆ t;τ)=δ(x), −∞ <x<∞, being δ(·) the Dirac delta function. Thus, in order to determine the131 1-PDF of the solution SP, xτ(t;ω), given by (11), we will only consider time instants tsuch that132 xτ(t;ω),0 a.e. As a consequence, as xτ(t;ω)=ea(ω)t+c(ω)eb1(ω),t τ, one derives that eb1(ω),t τ,0133 a.e., at every time instant twhere the 1-PDF f(x,t;τ) is going to be determined.134 Let t∈[(n−1)τ, nτ[ be fixed, next we will apply the RVT method (see Theorem 2), in order to135 obtain the PDF, f(x,t;τ), of the solution SP, xτ(t;ω), given by (11). This PDF will be expressed136 in terms of the joint PDF fc,a,b(c,a,b). To this end, consider the following deterministic mapping137 r:R3−→ R3 138 x1=r1(c,a,b)=eat+ceb1,t τ, x2=r2(c,a,b)=a, x3=r3(c,a,b)=b, where b1=e−aτb. The inverse mapping, s:R3−→ R3, of ris given by139 c=s1(x1,x2,x3)=ln x1e−x2t eb1,t τ!, a=s2(x1,x2,x3)=x2, b=s3(x1,x2,x3)=x3, where b1=e−x2τx3. Observe that the Jacobian of sis |J3|=1/|x1|,0. Moreover, |J3|is well-140 defined since xτ(t;ω)=x1,0 a.s. for each instant time t. Therefore, applying the RVT method141 one obtains the PDF of the random vector (x1(ω),x2(ω),x3(ω)),142 f(x1,x2,x3)=fc,a,b ln x1e−x2t eb1,t τ!,x2,x3!1 |x1|,where b1=e−x2τx3. Then, marginalizing with respect to the random vector (x2(ω),x3(ω)) =(a(ω),b(ω)) and taking143 t∈[(n−1)τ, nτ[ arbitrary, one gets the 1-PDF of the solution SP, xτ(t;ω),144 f(x,t;τ)=ZR2 fc,a,b ln xe−at eb1,t τ!,a,b!1 |x|dadb,where b1=e−aτb.(12) 2.2. Convergence145 As it has been indicated previously, this subsection is devoted to investigate the relationship146 between the 1-PDF, f(x,t;τ), given in (12), as τ→0+, and the 1-PDF, f(x,t), of the solution147 SP to random IPV (8). In order for the corresponding analysis makes sense (put τ→0+in the148 initial condition of (10)), we will take as initial condition in this latter IVP g0(ω)=ec(ω), so the149 solution SP to random IPV (8) is given by150 x(t;ω)=ec(ω)e(a(ω)+b(ω)) t. 6 By applying the RVT method, it can be checked that the 1-PDF of x(t;ω) is given by151 f(x,t)=ZR2 fc,a,bln xe−(a+b)t,a,b1 |x|dadb.(13) In order to find out sufficient conditions that guarantee the convergence given in (9), where152 f(x,t;τ) and f(x,t) are given by (12) and (13), respectively, we first establish the following result153 that permits relating the delayed exponential function to the (classical) exponential function as154 the delay tends to zero. The proof of this result is based upon the ideas exhibited in [22, Theorem155 A.3.].156 Theorem 3 (Convergence of the delayed exponential function). Let c ∈R,τ0>0,α=1+157 eτ0|c|>1. Then, for any τ∈]0, τ0[,158 |ect −ec,t τ| ≤ τ|c|eα|c|T+e|c|(T−τ),t∈[0,T]. If τ=0then ect =ec,t τ.159 Proof Let τ∈]0, τ0[. We will apply mathematical induction in order to prove that for any n∈N160 |ect −ec,t τ| ≤ τ|c|eα|c|nτ+e|c|(n−1)τ,for t∈[(n−1)τ, nτ].(14) •If n=1, then t∈[0, τ] and, by the definition of the delayed exponential function given in161 (3), ec,t τ=1+ct. Therefore, we must prove162 |ect −(1 +ct)| ≤ τ|c|eα|c|τ+1. Applying the Fundamental Calculus Theorem (FCT) and the Mean Value Theorem (MVT)163 for integration, one derives164 |ect −(1 +ct)|≤|ect −1|+|ct|FCT =Zt 0"d dx(ecx)#dx +|c|t≤Zt 0 d dx(ecx)dx+|c|t =|c|Zt 0 ecx dx+|c|tt≤τ ≤ |c|Zτ 0 ecx dx+|c|τMVT ≤ |c|τecδ+|c|τ(I) ≤ |c|τeα|c|τ+1. Now, we justify Step (I) previously applied.165 Step (I): By the MVT, δ∈[0, τ]. Let α > 1 , then166 cδ≤ |c|δ≤ |c|τ≤α|c|τ, and, as the exponential is an increasing function, ecδ≤eα|c|τ.167 •Now, assuming that claim (14) is satisfied for n≥1 (induction hypothesis), we will apply168 the FCT to prove inequality (14) for n+1. Let t∈[nτ, (n+1)τ],169 7 |ec,t τ−ect |(II) ≤ |c|τeα|c|nτ+e|c|(n−1)τ+Z(n+1)τ nτ d ds(ec,s τ−ecs)ds =|c|τeα|c|nτ+e|c|(n−1)τ+|c|Z(n+1)τ nτec,s−τ τ−ecsds ≤ |c|τeα|c|nτ+e|c|(n−1)τ+|c|Z(n+1)τ nτec,s−τ τ−ec(s−τ)ds +|c|Z(n+1)τ nτec(s−τ)−ecsds (III) ≤ |c|τeα|c|nτ+e|c|(n−1)τ+|c|Z(n+1)τ nτ|c|τeα|c|nτ+e|c|(n−1)τds +|c|Z(n+1)τ nτZs s−τ d dσ(ecσ)dσds =|c|τeα|c|nτ+e|c|(n−1)τ+(|c|τ)2eα|c|nτ+e|c|(n−1)τ +|c|2Z(n+1)τ nτZs s−τ ecσdσds (IV) ≤ |c|τeα|c|nτ+e|c|(n−1)τ+(|c|τ)2eα|c|nτ+e|c|(n−1)τ +|c|2Z(n+1)τ nτZs s−τ eα|c|nτ+|c|τdσds =|c|τeα|c|nτ+e|c|(n−1)τ+(|c|τ)2eα|c|nτ+e|c|(n−1)τ+|c|2τ2eα|c|nτ+|c|τ =|c|τeα|c|nτ1+|c|τ+|c|τe|c|τ+e|c|(n−1)τ(1 +|c|τ) (V) ≤ |c|τeα|c|(n+1)τ+e|c|nτ. Now, we justify Steps (II)–(V) previously applied.170 Step (II): By the FCT171 Zt nτ d ds(ec,s τ−ecs)ds=ec,t τ−ect −ec,nτ τ+ecnτ. 8 Then by the induction hypothesis and taking into account that t≤(n+1)τ, we have172 |ec,t τ−ect |≤|ec,nτ τ−ecnτ|+Zt nτ d ds(ec,s τ−ecs)ds ≤ |c|τeα|c|nτ+e|c|(n−1)τ+Z(n+1)τ nτ d ds(ec,s τ−ecs)ds. Step (III): In this step we apply both the FCT and the hypothesis of induction (14).173 Step (IV): Following the same argument of Step (I), being α > 1174 cσ≤ |c|σ≤ |c|s≤ |c|(n+1)τ≤α|c|nτ+|c|τ, so, ecσ≤eα|c|nτ+|c|τ.175 Step (V): Applying that ep≥1+p,∀p≥0 and taking α=1+e|c|τ0, as 0 < τ ≤τ0, then176 1+|c|τ+|c|τe|c|τ=1+|c|τ(1 +e|c|τ)≤1+|c|τ(1 +e|c|τ0)=1+α|c|τ≤eα|c|τ. And analogously, 1 +|c|τ≤e|c|τ.177 This concludes the proof.178 179 Before showing the convergence (9), we need to establish some technical results that will be180 required in the subsequent analysis.181 Let c∈Rbe arbitrary and let us take τ0→0+(thus τ→0+, too) in Theorem 3. Then182 lim τ→0+ec,t τ=ect,t∈[0,T],(15) where, according to Definition 3, t=(n−1)τ,n→+∞and τ→0+. Furthermore, as the delayed183 exponential function is continuous, one derives that184 if lim τ→0+f(τ)=ˆ f,then lim τ→0+ef(τ),t τ=eˆ f t ,0.(16) In particular, if we take f(τ)=e−aτbwith a,b∈Rfixed, as e−aτbτ→0+ −→ bin (16), then185 lim τ→0+ee−aτb,t τ=ebt ,0. Therefore, there exists τ1>0 such that186 eb1,t τ,0,b1=e−aτb,∀τ∈[0, τ1]. As a consequence of the continuity of the logarithm function and of the delayed exponential187 function (15), one further obtains188 lim τ→0+ln xe−at eb1,t τ!=ln xe−(a+b)t. 9 which is just the right-hand side of (VII).258 Step (VIII): In this step we will prove that (A) and (D) are bounded, and (B) and (C) tend to zero259 as τ→0+.260 •Let us see that expression (B) tends to zero as τ→0+.261 According to (20), for each ω∈Ω, given a(ω)=aand b(ω)=b, it is known that262 ϕ0 τ(t,a,b)=ea(t+τ)eb1,t τ−aR0 −τea(t−s)eb1,t−τ−s τds, then by the MVT and the convergence of263 the delayed exponential function to the exponential one, Theorem 3, we have264 lim τ→0+ϕ0 τ(t,a,b)=e(a+b)t. By continuity, it is clear that 1/ϕ0 τ(t,a,b)τ→0+ −→ 1/e(a+b)t, then for all 0>0 there exists τ1 265 such that for every τ∈]0, τ1[∩]0, τ0[ (being τ0defined in Theorem 3),266  1 ϕ0 τ(t,a,b)−1 e(a+b)t < 0.(25) •Let us see that expression (A) is bounded.267 Let F0=fc0(0) >0, then by the Lipschitz condition, hypothesis ˆ H1,268 fc0        x−Pm j=1ϕj τ(t,a,b)cj ϕ0 τ(t,a,b)       ≤Lfc0 1 |ϕ0(t,a,b;τ)||x|+ m X j=1|ϕj τ(t,a,b)||cj|+F0 ≤Lfc00+e|a|Te|b|T       |x|+ m X j=1 j|cj|       +F0, where last inequality is justified by formula (25),269  1 ϕ0 τ(t,a,b)≤ 1 ϕ0 τ(t,a,b)−1 e(a+b)t +1 e(a+b)t< 0+e−(a+b)t≤0+e|a|Te|b|T,t∈[0,T]. On the other hand, by (21), ϕj τ(t,a,b)=ea(t+τ)eb1,t τ(−τ)j+R0 −τea(t−s)eb1,t−τ−s τ(j−as)sj−1ds,270 1≤j≤m, then by the MVT, the convergence of the delayed exponential function to the271 exponential function and, Theorem 3, we have272 lim τ→0+ϕj τ(t,a,b)=0,t∈[0,T]. Thus, for each j=1,...,mand j>0 arbitrary, there exists τj+1>0 such that for every273 τ∈]0, τj+1[∩]0, τ0[274 ϕj τ(t,a,b)< j.(26) •Let us see that expression (C) tends to zero as τ→0+.275 16 By the Lipschitz condition, hypothesis ˆ H1, and expressions (25)–(26), one derives276  fc0        x−Pm j=1ϕj τ(t,a,b)cj ϕ0 τ(t,a,b)       −fc0xe−(a+b)t≤Lfc0 x−Pm j=1ϕj τ(t,a,b)cj ϕ0 τ(t,a,b)−xe−(a+b)t ≤Lfc0       |x| 1 ϕ0 τ(t,a,b)−e−(a+b)t +1 |ϕ0 τ(t,a,b;τ)| m X j=1|ϕj τ(t,a,b)||cj|        <Lfc0       |x|0+0+e|a|Te|b|Tm X j=1 j|cj|       , with τ∈ ∩m+1 j=0]0, τj[.277 •Let us see that expression (D) is bounded. Indeed, it is clear that278 e−(a+b)t≤e|a|Te|b|T,t∈[0,T]. Then, the right hand-side of the inequality of Step (VIII) is obtained.279 Finally, we use the definition of the expectation and the independence between the RVs a(ω),280 b(ω) and cj(ω), j=1,2,...,m, obtaining the last expression.281 Now, we assume the following hypothesis in order to prove the convergence282 ˆ H3 :Ee|a(ω)|T2=ka<+∞,Ee|b(ω)|T2=kb<+∞,Eh|cj|i=kj<+∞. Notice that if Eh(y(ω))2i=ky<∞, then by the Cauchy-Schwarz inequality Ey(ω)≤283 Eh(y(ω))2i1/2=k1/2 y<∞.284 Therefore, for τ∈]0, τ∗[, τ∗=min{τj: 0 ≤j≤m}, (x,t)∈ D(xτ(t;ω)) ∩D(x(t;ω)) ×[(n−285 1)τ, nτ[⊂R×[0,T] all of them being fixed286 |f(x,t;τ)−f(x,t)| ≤ 0       Lfc00+Ehe|a(ω)|TiEhe|b(ω)|Ti       |x|+ m X j=1 jEh|cj(ω)|i       +F0        +Lfc0|x|0Ehe|a(ω)|TiEhe|b(ω)|Ti +Lfc00Ehe|a(ω)|TiEhe|b(ω)|Ti+Ehe2|a(ω)|TiEhe2|b(ω)|Ti m X j=1 jEh|cj(ω)|i ≤0       Lfc00+k1/2 ak1/2 b       |x|+ m X j=1 jkj       +F0        +Lfc0       |x|0k1/2 ak1/2 b+0k1/2 ak1/2 b+kakbm X j=1 jkj       . 17 Summarizing, the following result has been established.287 Theorem 5. Consider the random discrete delay differential equation (4) with g(t;ω)given288 by (18) and whose solution SP, xτ(t;ω), is given by (19)–(21). Let us assume that w(ω)=289 (c0(ω),c1(ω),...,cm(ω),a(ω),b(ω)), is an absolutely continuous random vector being fw(w)290 their joint PDF. Then, the 1-PDF, f (x,t;τ), of the solution SP xτ(t;ω)is given by (22). Fur-291 ther assume that hypotheses ˆ H1–ˆ H3 hold, then f (x,t;τ)converges to the 1-PDF, f (x,t), of the292 solution SP to the random linear differential equation (23), according to (9).293 4. Numerical examples294 This section is devoted to illustrate our theoretical findings by means of several numerical295 experiments where the 1-PDF of the solution SP to random IVP (4) is computed in the two cases296 previously investigated. For the sake the clarity in the presentation, the explicit expressions297 of the 1-PDFs, in each one of the numerical examples, are reported in Appendix B since their298 mathematical representation are somewhat cumbersome. It is important to point out that we have299 chosen a wide variety of probability distributions for the input parameters a(ω), b(ω) and c(ω).300 Example 1. Let us consider random IVP (4) whose initial condition has the form g(t;ω)=301 ea(ω)t+c(ω), which corresponds to Case I. We will assume that a(ω), b(ω)and c(ω)are independent302 continuous absolutely RVs (hence hypothesis H1 is fulfilled) with the following distributions:303 •a(ω)is a Gaussian RV with zero mean and standard deviation 0.1, i.e., a(ω)∼N(0; 0.1).304 •b(ω)is a Beta RV with parameters 2and 3, i.e., b(ω)∼Be(2; 3).305 •c(ω)is an Exponential RV with mean 1/20, i.e., c(ω)∼Exp(20).306 Since the first derivative of the PDF of RV c(ω)is bounded, hypothesis H2 also holds. Therefore,307 assumptions of Theorem 4 are fulfilled. Now, we check numerically that the 1-PDF, f (x,t;τ),308 of the solution SP of random IVP (4) with g(t;ω)=ea(ω)t+c(ω)converges to the 1-PDF, f (x,t),309 of the solution SP of random IVP (8). To this end, in Figure 1 we have plotted f (x,t)together310 with f (x,t;τ)with different delays τ∈ {0.01,0.05,0.1,0.5,2}at different time instants t =0.1,311 t=0.5and t =1. From this graphical representation we can observe that f (x,t;τ)converges312 rapidly to f (x,t)as τ→0+. To numerically assess this convergence, in Table 1 we show the313 error between f (x,t)and f (x,t;τ)for the values of the delays and the time instants previously314 indicated, according to the following error measure315 ePDF τ(t)=ZR|f(x,t;τ)−f(x,t)|dx.(27) We observe that for t fixed, the error ePDF τ(t)decreases as τ→0+, as expected. We also observe316 that the velocity of the convergence decreases as t increases.317 Example 2. In this second example we consider that the initial condition in random IVP (4) is318 a polynomial of degree m, g(t;ω)=Pm j=0cj(ω)tj,m≥0which corresponds with Case II studied319 before. We will consider two problems, when the initial condition is a constant RV (i.e., random320 polynomial of degree m =0) and when is a random polynomial of degree m =1.321 Problem 1: g(t;ω)=c0(ω).322 323 Let T =0.5and assume that a(ω), b(ω)and c0(ω)are independent RVs (hence hypothesis324 ˆ H1 is fulfilled) with the following distributions:325 18 t=0.1 f1(x,t) f1(x,t;2) f1(x,t;0.5) f1(x,t;0.1) f1(x,t;0.05) f1(x,t;0.01) 1.055 1.060 1.065 1.070 x 8.6 8.8 9.0 9.2 1.0 1.2 1.4 1.6 1.8 2.0 2.2 x 2 4 6 8 10 t=0.5 f1(x,t) f1(x,t;2) f1(x,t;0.5) f1(x,t;0.1) f1(x,t;0.05) f1(x,t;0.01) 1.20 1.22 1.24 1.26 1.28 x 2.45 2.50 2.55 2.60 2.65 1.0 1.5 2.0 2.5 3.0 x 0.5 1.0 1.5 2.0 2.5 3.0 t=1 f1(x,t) f1(x,t;2) f1(x,t;0.5) f1(x,t;0.1) f1(x,t;0.05) f1(x,t;0.01) 1.42 1.44 1.46 1.48 1.50 x 1.10 1.12 1.14 1.16 1.18 1 2 3 4 x 0.5 1.0 1.5 Figure 1: Graphical representation of the PDFs, f(x,t) and f(x,t;τ), with different delays τ∈ {0.01,0.05,0.1,0.5,2}, at different time instants: t=0.1 (top), t=0.5 (center) and t=1 (bottom) in the context of Example 1. 19 ePDF τ(t)τ=2τ=0.5τ=0.1τ=0.05 τ=0.01 t=0.10 0.040793 0.027599 0.021251 0.015753 0.003949 t=0.50 0.215938 0.172538 0.061891 0.032553 0.006153 t=1.00 0.382194 0.252025 0.069224 0.035766 0.007786 Table 1: Error measure ePDF τ(t), defined by (27), with different delays τ∈ {0.01,0.05,0.1,0.5,2}, at different time instants, t∈ {0.1,0.5,1}, in the context of Example 1. •a(ω)follows a truncated Gaussian distribution on the interval I =[−1,1] with zero mean326 and standard deviation 0.1, i.e., a(ω)∼NI(0; 0.1).327 •b(ω)is an Exponential RV with mean 1/50, i.e., b(ω)∼Exp(50).328 •c0(ω)is a Beta RV with parameters 2and 3, i.e., c0(ω)∼Be(2; 3).329 Since f 0 c0(c0)is bounded, hypothesis ˆ H2 also holds. Finally, we compute the three values330 given in hypothesis ˆ H3 with T =0.5, obtaining331 ka=1.08507,kb=1.02041,k0=0.4, which are all finite. Therefore, assumptions of Theorem 5 hold. Now, as in the previous example,332 in order to see the convergence of f (x,t;τ)to f (x,t), in Figure 2 we have plotted f (x,t)together333 with f (x,t;τ)with different delays τ∈ {0.1,0.5,1,2,3}at the time instants t =0.1and t =0.5.334 In addition, in Table 2 we have calculated the error given in formula (27). From both we can see335 graphical and numerical the convergence as τ→0+.336 t=0.1 f(x,t) f(x,t;3) f(x,t;2) f(x,t;1) f(x,t;0.5) f(x,t;0.1) 0.702 0.704 0.706 0.708 0.710 x 0.66 0.68 0.70 0.72 0.74 0.76 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x 0.5 1.0 1.5 t=0.5 f(x,t) f(x,t;3) f(x,t;2) f(x,t;1) f(x,t;0.5) f(x,t;0.1) 0.702 0.704 0.706 0.708 0.710 x 0.66 0.68 0.70 0.72 0.74 0.76 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x 0.5 1.0 1.5 Figure 2: Graphical representation of the PDFs f(x,t) and f(x,t;τ), with different delays τ∈ {0.1,0.5,1,2,3}, at different time instants t=0.1 (left) and t=0.5 (right) in the context of Problem 1 in Example 2. In addition, in Figure 3 the mean and the variance of x(t;ω)and xτ(t;ω)for different τ∈337 {0.1,0.5,1,2,3}in the time interval [0,0.5] have been represented. To calculate the mean and338 the variance of xτ(t;ω)we have used the expressions of E[xτ(t;ω)] and V[xτ(t;ω)], respectively,339 given in (6), where f (x,t;τ)is defined by (22). In a similar way we have computed the mean and340 the variance of x(t;ω)but using (24). We can observe the convergence as τ→0+. For sake of341 clarity, the error defined by (28) has been calculated in Table 3.342 20 ePDF τ(t)τ=3τ=2τ=1τ=0.5τ=0.1 t=0.10 0.425626 0.261867 0.098525 0.028427 0.001875 t=0.50 0.449224 0.287860 0.115015 0.043369 0.005402 Table 2: Error measure ePDF τ(t), defined by (27), with different delays τ∈ {0.1,0.5,1,2,3}, at different time instants, t∈ {0.1,0.5}, in the context of Problem 1 of Example 2. eE τ=ZT 0|E[xτ(t;ω)] −E[x(t;ω)]|dt,eV N=ZT 0|V[xτ(t;ω)] −V[x(t;ω)]|dt.(28) [x(t,ω)] [x3(t,ω)] [x2(t,ω)] [x1(t,ω)] [x0.5 (t,ω)] [x0.1 (t,ω)] 0.152 0.154 0.156 0.158 0.160 t 0.402 0.403 0.404 0.405 0.406 0.0 0.1 0.2 0.3 0.4 0.5 t0.35 0.40 0.45 0.50 0.55 [x(t,ω)] [x3(t,ω)] [x2(t,ω)] [x1(t,ω)] [x0.5 (t,ω)] [x0.1 (t,ω)] 0.152 0.154 0.156 0.158 0.160 t 0.042 0.044 0.046 0.048 0.050 0.0 0.1 0.2 0.3 0.4 0.5 t0.0 0.1 0.2 0.3 0.4 0.5 Figure 3: Graphical representation of the mean (left) and the variance (right) of the solutions SP x(t;ω) and xτ(t;ω), with different delays τ∈ {0.01,0.05,0.1,0.5,2}, at the time interval [0,1] in the context of Problem 1 in the Example 2. τ=3τ=2τ=1τ=0.5τ=0.1 eE τ0.0218545 0.0101997 0.00301312 0.000995992 0.00011704 eV τ0.049849 0.0214272 0.00574011 0.00172316 0.000165323 Table 3: Error measures eE τand eV τ, defined by (28), with different delays τ∈ {0.1,0.5,1,2,3}, in the context of Problem 1 in Example 2. Problem 2: g(t;ω)=c0(ω)+c1(ω)t.343 344 We take T =0.5and we will assume that a(ω), b(ω), c0(ω)and c1(ω)are independent RVs345 (hence hypothesis ˆ H1 is fulfilled) with the following distributions:346 •a(ω)is a Beta RV with parameters 2and 3, i.e., a(ω)∼Be(2; 3).347 •b(ω)follows a Gaussian distribution with mean −1and standard deviation 0.2, i.e., b(ω)∼348 N(−1; 0.2).349 •c0(ω)is a Uniform RV in the interval [−1,1], i.e., c0(ω)∼U([−1,1]).350 21 •c1(ω)follows an Exponential distribution with mean 1/20, i.e., c1(ω)∼Exp(20).351 Since c0(ω)follows a Uniform distribution, then its PDF is a constant and therefore the352 derivative of the PDF is 0, then hypothesis ˆ H2 also holds. Finally, we compute the three values353 given in hypothesis ˆ H3 with T =0.5, obtaining354 ka=1.52247,kb=2.77319,k0=0.5,k1=0.05, which are all finite. Therefore, assumptions of Theorem 5 are satisfied. Now, as in the previous355 example, in order to see the convergence of f (x,t;τ)to f (x,t), in Figure 4 we have plotted f (x,t)356 together with f (x,t;τ)with different delays τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}at time instants t =357 0.1and t =0.5. In addition, in Table 4 we have calculated the error given in formula (27). From358 both we can see graphical and numerical the convergence as τ→0+.359 t=0.1 f(x,t) f(x,t;2.5) f(x,t;2) f(x,t;1.5) f(x,t;1) f(x,t;0.5) f(x,t;0.1) f(x,t;0.05) -0.04 -0.02 0.02 0.04 x 0.535 0.540 0.545 0.550 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 x 0.5 1.0 1.5 t=0.5 f(x,t) f(x,t;2.5) f(x,t;2) f(x,t;1.5) f(x,t;1) f(x,t;0.5) f(x,t;0.1) f(x,t;0.05) -0.04 -0.02 0.02 0.04 x 0.685 0.690 0.695 0.700 0.705 0.710 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 x 0.5 1.0 1.5 Figure 4: Graphical representation of the PDFs f(x,t) and f(x,t;τ), with different delays τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}, at different time instants t=0.1 (left) and t=0.5 (right) in the context of Problem 2 in Example 2. ePDF τ(t)τ=2.5τ=2τ=1.5τ=1 t=0.10 0.853492 0.666295 0.454429 0.236062 t=0.50 0.808335 0.639961 0.440874 0.248945 ePDF τ(t)τ=0.5τ=0.1τ=0.05 t=0.10 0.0631106 0.00712427 0.00581714 t=0.50 0.182441 0.0520957 0.0178104 Table 4: Error measure ePDF τ(t), defined by (27), with different delays τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}, at different time instants: t∈ {0.10,0.50}, in the context of Problem 2 in Example 2. In Figure 5 the mean and the variance of f (x,t)and f1(x,t;τ)for different τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}360 in the time interval [0,0.5] have been represented. We can observe the convergence as τ→0+.361 As in the previous Problem 1, in Table 5 we show the error for the mean and the variance defined362 by (28). We observe that the difference of these errors becomes smaller as τgoes to zero.363 22 [x(t,ω)] [x2.5 (t,ω)] [x2(t,ω)] [x1.5 (t,ω)] [x1(t,ω)] [x0.5 (t,ω)] [x0.1 (t,ω)] [x0.05 (t,ω)] 0.1520.1540.1560.1580.160 t 0.00005 0.00010 0.00015 0.00020 0.00025 0.1 0.2 0.3 0.4 0.5 t0.00 0.05 0.10 0.15 0.20 0.25 0.30 [x(t,ω)] [x2.5 (t,ω)] [x2(t,ω)] [x1.5 (t,ω)] [x1(t,ω)] [x0.5 (t,ω)] [x0.1 (t,ω)] [x0.05 (t,ω)] 0.152 0.154 0.156 0.158 0.160 t 0.280 0.285 0.290 0.295 0.300 0.305 0.0 0.1 0.2 0.3 0.4 0.5 t0.0 0.2 0.4 0.6 0.8 1.0 Figure 5: Graphical representation of the mean (left) and the variance (right) of the solutions SP x(t;ω) and xτ(t;ω), with different delays τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}, at the time interval [0,0.5] in the context of Problem 2 of Example 2. τ=2.5τ=2τ=1.5τ=1 eE τ0.0642834 0.0434024 0.0258283 0.0121818 eV τ0.11049 0.0967728 0.0742367 0.0479234 τ=0.5τ=0.1τ=0.05 eE τ0.0033024 0.000101439 0.000011333 eV τ0.0282562 0.0208731 0.0194024 Table 5: Error measures eE τand eV τ, defined by (28), with different delays τ∈ {0.05,0.1,0.5,1,1.5,2,2.5}, in the context of Problem 2 in Example 2. 23 5. Conclusions364 The main goal of the paper has been to compute the first probability density function (1-PDF)365 of the solution stochastic process of an important class of random linear differential equations366 with discrete delay. The study has been split into two cases depending on the functional form of367 the initial condition. It is important to point out that the success of our approach relies upon the368 knowledge of an explicit expression of the solution. Otherwise, RVT method would need to be369 combined with numerical strategies. In the examples, we have shown how to compute the mean370 and the variance of the solution stochastic process from the 1-PDF. It is important to point out371 that the analysis performed in this contribution is also useful to deal with other types of initial372 conditions different from the ones presented here. To the best of our knowledge, this is the first373 time that this class of random differential equations with delay are studied following the proposed374 approach. In a future contribution, we plan to go a step forward by studying the case which the375 coefficients of the differential equation are stochastic processes rather than random variables.376 Acknowledgements377 This work has been partially supported by the Ministerio de Econom´ ıa y Competitividad378 grant MTM2017–89664–P and MTM2015-63723-P and Junta de Andaluc´ ıa under Proyecto de379 Excelencia P12-FQM-1492. Ana Navarro Quiles acknowledges the Fundaci´ o Ferran Sunyer i380 Balaguer and the Instituto de Estudios Catalanes for its contribution through the Borsa Ferran381 Sunyer i Balaguer. Ana Navarro Quiles acknowledges the postdoctoral contract financed by Dy-382 Con project funding from the European Research Council (ERC) under the European Unions383 Horizon 2020 research and innovation programme (grant agreement No 694126-DYCON).384 The authors express their deepest thanks and respect to the editors and reviewers for their valu-385 able comments.386 Conflict of Interest Statement387 The authors declare that there is no conflict of interests regarding the publication of this388 article.389 390 Appendix A. Proof of Theorem 1: Existence and uniqueness of the solution391 First, we will show how to construct the solution given in expression (2). To this end, first it392 is convenient to construct the solution of the following IVP393 (x0 τ(t)=bxτ(t−τ),t>0, τ > 0, xτ(t)=1,−τ≤t≤0,(A.1) which corresponds to the IVP (1) with a=0 and g(t)=1. We will see that its solution is the394 delayed exponential function, xτ(t)=eb,t τ. This fact will check in the subintervals [−τ, 0], [0, τ],395 [τ, 2τ], . . .:396 •If t∈[−τ, 0], it is clear, taken into account the definition of the delayed exponential func-397 tion (see (3)), that xτ(t)=1=eb,t τ,−τ≤t≤0,satisfies the IVP (A.1).398 24 •If t∈[0, τ], we integrate the DDE in (A.1) on the interval [0,t]399 Zt 0 x0 τ(s) ds=Zt 0 bxτ(s−τ) ds. As s∈[0,t] and t∈[0, τ], s−τ∈[−τ, 0], then by the definition of the delayed exponential400 function xτ(s−τ)=1. Therefore, applying the FCT401 xτ(t)−xτ(0) =Zt 0 bdsxτ(0)=1 =⇒xτ(t)=1+bt =eb,t τ,0≤t≤τ. •If t∈[τ, 2τ], we integrate the DDE in (A.1) on the interval [τ, t]402 Zt τ x0 τ(s) ds=Zt τ bxτ(s−τ) ds. Analogously, s∈[τ, t] and t∈[τ, 2τ] implies s−τ∈[0, τ], then, by the previous step,403 xτ(s−τ)=1+b(s−τ). Therefore, applying the FCT and integrating one gets404 xτ(t)−xτ(τ)=Zt τ b(1+b(s−τ)) dsxτ(τ)=1+bτ =⇒xτ(t)=1+bt+b2(t−τ)2 2=eb,t τ, τ ≤t≤2τ. Following this reasoning, it is straightforward we conclude by induction that the solution of IVP405 (A.1) is given by xτ(t)=eb,t τ, for each t∈[(n−1)τ, nτ], n=0,1,2, . . ..406 On the other hand, suppose that the initial condition in (A.1) is a time-dependent function xτ(t)=407 g(t), with −τ≤t≤0. By the method of variation of constants [], the solution can be written as408 xτ(t)=x∗ τ(t)c+Z0 −τ x∗ τ(t−τ−s)φ0(s) ds,(A.2) being x∗ τ(t)=eb,t τthe solution of (A.1), can unknown constant and φ(s) is an unknown continu-409 ously differentiable function which must be determined. If −τ≤t≤0, then xτ(t)=g(t). Thus,410 substituting in (A.2)411 g(t)=x∗ τ(t)c+Z0 −τ x∗ τ(t−τ−s)φ0(s) ds,−τ≤t≤0. We split this integral into the two following integrals412 g(t)=x∗ τ(t)c+Zt −τ x∗ τ(t−τ−s)φ0(s) ds+Z0 t x∗ τ(t−τ−s)φ0(s) ds. Now, by the properties of the delayed exponential function x∗ τ(t)=eb,t τone gets413 • −τ≤t≤0 implies x∗ τ(t)=1.414 •When s∈[−τ, t], then t−τ−s∈[−τ, 0]. Therefore, x∗ τ(t−τ−s)=1.415 •When s∈[t,0], then t−τ−s∈[−2τ, −τ]. Therefore, x∗ τ(t−τ−s)=0 if s,tand416 x∗ τ(t−τ−s)=1 if s=t.417 25