Mean square calculus and random linear fractional differential equations: Theory and applications

It is well-known that the following linear differential equation with initial condition is a proper pattern to model a great number of physical problems

y(i)(t)−λy(t)=c,t>0,y(j)(0)=bj,0≤j<i,i=1,2,j∈ℕ. $$ \begin{equation}\label{IVP_deterministico} y^{(i)}(t)-\lambda y(t)=c,\quad t\gt0,\qquad y^{(j)}(0)=b_j, \qquad 0 \leq j \lt i, \qquad i=1,2,\quad j \in \mathbb{N}. \end{equation} $$(1)

Here ℕ $\mathbb{N}$ denotes the set of positive integers and λ,c,b0∈ℝ $\lambda,c,b_0 \in \mathbb{R}$ . In practice these model parameters need to be fixed using measurements, hence they contain measurement errors. This fact motivates that it is more realistic to treat them as random variables (RVs) instead of deterministic values. Apart from measurement errors, randomness can be attributed to the inherent complexity often encountered in physic phenomena. This justifies the consideration of variability in the forcing term and in the initial condition (c, b₀, b₁). It must be pointed out that randomness may be included also in the diffusion coefficient, however in this contribution, for simplicity, we will assume that λ is a constant.

On the other hand, over the last few decades fractional differential equations are playing a significant role in modelling phenomena with microscopic complex behaviour whose macroscopic dynamics (like viscoelasticity, fluid behaviour, etc.) cannot be properly described using classical derivatives [4]. Throughout this paper, randomness will be considered using the Mean Square Calculus that involves the so-called mean square convergence. In the context of ordinary and fractional random differential equations, this random Calculus has been successfully used to conduct a number of significant problems. Some examples can be checked in [3, 5, 6, 7]. The application mean square calculus to deal with the analysis of random fractional differential equations is adequate because it permits to take advantage of properties of mean square convergence to compute approximations of the expectation and variance of the solution stochastic process. Our contribution is based upon a number of results belonging to the mean square calculus that have been established in the extant literature.

In this contribution we deal with the following generalization of the initial value problem (IVP) (1)

CD0+α<!-- a -->Y(t)−<!-- - -->λ<!-- ? -->Y(t)=c,t><!-- > -->0,0<<!-- < -->α<!-- a -->≤<!-- = -->2,Y(j)(0)=bj,0≤<!-- = -->j<<!-- < -->[α<!-- a -->]−<!-- - -->1, $$ \begin{equation} \left\{\begin{array}{ccl} \left(^CD_{0^{+}}^\alpha Y\right)(t) -\lambda Y(t)& = & c, \quad t\gt0,\,\,\, \hspace{0.2cm} 0 \lt \alpha \leq 2,\\ Y^{(j)}(0) &=& b_j, \quad 0 \leq j \lt[\alpha]-1, \end{array} \right. \end{equation} $$(2)

being [·] the integer part function. To conduct our analysis, we will assume that λ∈ℝ $\lambda\in \mathbb{R}$ and b_j and c are independent second-order RVs (2-RVs), that is, we assume that b_j and c have finite variance. The symbol (CD0+αY)(t) $\left(^CD_{0^{+}}^\alpha Y\right)(t)$ denotes the random Caputo fractional derivative of α order of the stochastic process (SP), Y(t). Although there exist several definitions of derivatives of non-integer order (fractional derivatives), for convenience, throughout our analysis we will consider the Caputo derivative. For the sake of completeness, we recall the definition of this random derivative.

In [8, ch.4], one can find the definition of mean square differentiability to 2-SPs. The random fractional integral (J^n−αf)(t) is a natural extension of its deterministic counterpart (see [4]) using the properties of the mean square Riemann integral of 2-SPs (see [8, ch.4]).

We recall that the correlation function, Γ_f(u₁, u₂), of a 2-SP, f(t), is defined by

Γf(u1,u2):=E[f(u1)f(u2)],u1,u2∈ℝ, $$ \begin{equation*} \Gamma_f(u_1,u_2):=\mathbb{E}[f(u_1)f(u_2)],\quad u_1,u_2\in \mathbb{R}, \end{equation*} $$

where E[⋅] $\mathbb{E}[\cdot]$ denotes the expectation operator. It can be proved, by applying the Schwarz inequality, that Γ_f(u₁, u₂) always is well-defined. An important property of the correlation function of a 2-SP, f(t) and its n-th mean square derivative, f⁽ⁿ⁾(t), is the following

Γf(n)(u1,u2)=∂2nΓf(u1,u2)∂u1n∂u2n,u1,u2∈ℝ. $$ \begin{equation} \Gamma_{f^{(n)}}(u_1,u_2)= \dfrac{\partial^{2n} \Gamma_f(u_1,u_2)} {\partial u_1^n \, \partial u_2^n},\quad u_1,u_2\in \mathbb{R}. \end{equation} $$(4)

Using (4) and the characterization for the existence of the mean square Riemann integral of 2-SP f(t) (see [8, Th. 4.5.1]), one can establish following necessary and sufficient condition for the existence of the random Caputo fractional mean square derivative.

Once we have defined the Caputo derivative of a 2-SP, in Sections 2 and 3 we will construct the solution of random IVP (2) in two steps. First, considering the case when 0<α≤1 $0 \lt \alpha \leq 1$ (Case I) and, secondly, when 1<α≤2 $1 \lt \alpha \leq 2$ (Case II). As we shall see throughout analysis, this latter case will be based upon the study carried out in the former.

In this section we will deal with the construction of the solution SP to the IVP (2) in the case that 0<α≤1 by means of a random generalized power series. It will be done combining the Fröbenious method together with a mean square chain rule for differentiating 2-SPs. This result has been recently established by some of the authors [1].

The solution to the random IVP (2) will be seek in the following form

Y(t)=∑m≥0Xmtαm,t≥0,0<α≤1. $$ \begin{equation} Y(t)=\sum_{m\geq 0} X_m t^{\alpha m},\quad t\geq 0, \qquad 0\lt \alpha \leq 1. \end{equation} $$(6)

If we define

X(t)=∑m≥0Xmtm, $$ \begin{equation} X(t) = \sum_{m\geq 0} X_m t^{m}, \end{equation} $$(7)

then Y(t) = X(t^α), and in this manner we can take advantage of mean square calculus for standard random power series. The random Caputo mean square derivative (3) writes

(CDa+αY)(t)=1Γ(1−α)∫0t(t−u)−αZ(u)du, $$ \begin{equation} \left(^CD_{a+}^\alpha Y\right)(t) = \frac{1}{\Gamma(1-\alpha)}\int_0^t (t-u)^{-\alpha} Z(u) \mathrm{d}u, \end{equation} $$(8)

where Z(t) : (X(t^α))′ is the mean square derivative of the SP X(t) compounded with the deterministic function t^α.

If we want to apply the Fröbenius method, firstly we need to obtain the Caputo derivative of the generalized power series given by (6). This is just done by applying Theorem 2 assuming the following hypotheses:

C1

X(v), defined in expression (7), is differentiable in mean square sense at v = t^α. Also, X′(tα)=∑m≥1mXmtα(m−1). $$ \begin{equation} X'(t^{\alpha}) = \sum_{m\geq 1}m X_m t^{\alpha (m-1)} . \end{equation} $$(9)

As 0 < α ≤ 1, one gets that V=[0,Tα] $\mathcal{V}=[0,T^{\alpha}]$ contains the range of g(t) = t^α, that is g([a1,a2])=g([0,T])=[0,Tα]=V $g([a_1,a_2])=g([0,T])=[0,T^{\alpha}]=\mathcal{V}$ . By Theorem 2 we can check that X(g (t)) is mean square differentiable at t and its mean square derivative has the following expression

Z(t):=Y′(t)=(X(tα))′=αtα−1X′(tα). $$ \begin{equation} Z(t):=Y'(t)=(X(t^{\alpha}))'= \alpha t^{\alpha-1} X'(t^{\alpha}). \end{equation} $$(10)

Therefore

(CDa+αY)(t)=1Γ(1−α)∫0t(t−u)−α∑m=1∞αmXmuαm−1du. $$ \begin{equation} \left(^CD_{a+}^\alpha Y\right)(t) = \frac{1}{\Gamma(1-\alpha)}\int_0^t (t-u)^{-\alpha} \sum_{m=1}^\infty \alpha m X_m u^{\alpha m-1}\mathrm{d}u. \end{equation} $$(11)

Furthermore, assuming that C3

The random power series ∑m=1∞mXmtαm−1 $\sum_{m=1}^\infty m X_m t^{\alpha m-1}$ is mean square convergent on the interval 0 ≤ t ≤ T,

the infinite sum and the integral in (11) can be commuted and then it can be proved that the resulting expression writes as

(CD0+αY)(t)=∑m≥0(Xm+1Γ(α(m+1)+1)Γ(αm+1)tαm). $$ \begin{equation} \begin{split} \left(^CD_{0^{+}}^\alpha Y \right)(t) & = \sum_{m\geq 0} \left( X_{m+1} \dfrac{\Gamma(\alpha (m+1)+1)}{\Gamma(\alpha m+1)} t^{\alpha m} \right). \end{split} \end{equation} $$(12)

In order to calculate the coefficients X_m of Y(t), we need to take into account that Y (0) = X₀ = b₀. Then, using the Fröbenius method one gets

Xm+1=λm+1b0+λmcΓ(α(m+1)+1),m≥0. $$ \begin{equation*} X_{m+1} = \dfrac{\lambda^{m+1} b_0 + \lambda^{m} c}{\Gamma(\alpha (m+1) +1)}, \quad m\geq 0. \end{equation*} $$

A candidate solution to the random IVP (2) when 0 < α ≤ is given by

Y(t)=X(tα),X(v)=∑m≥0Xm,1vm+∑m≥1Xm,2vm, $$ \begin{equation} Y(t)=X(t^{\alpha}),\qquad X(v)=\sum_{m\geq 0} X_{m,1} v^m+ \sum_{m\geq 1} X_{m,2} v^m, \end{equation} $$(13)

where

Xm,1=λmb0Γ(αm+1),Xm,2=λm−1cΓ(αm+1). $$ \begin{equation} X_{m,1} = \dfrac{\lambda^m b_0}{\Gamma(\alpha m +1)},\qquad X_{m,2} = \dfrac{\lambda^{m-1} c}{\Gamma(\alpha m +1)}. \end{equation} $$(14)

This solution must verify conditions C1–C3 in order to legitimate the assumptions that allows us to define the random Caputo mean square derivative of the 2-SP Y(t). In order to check that C1 is satisfied, we will need the following result:

Now, we check the hypotheses of Proposition 3 for the two series defined in (13)–(14). First, let us define

Xm,1(t):=Xm,1tm=λmb0Γ(αm+1)tm, $$ \begin{equation*} X_{m,1}(t):= X_{m,1}t^m=\dfrac{\lambda^m b_0}{\Gamma(\alpha m + 1)} t^m, \end{equation*} $$

and observe that

0<‖Xm,1(tα+h)−Xm,1(tα)h−X′m,1(tα)‖2≤|λ|m‖b0‖2Γ(αm+1)|(tα+h)m−tαmh−mtα(m−1)|→h→00, $$ \begin{equation} \begin{split} 0\lt \left\| \dfrac{X_{m,1}(t^{\alpha}+h)-X_{m,1}(t^{\alpha})}{h}-X_{m,1}'(t^{\alpha}) \right\|_2 % %&\ \leq \dfrac{|\lambda|^m \left\| b_0 \right\|_2}{\Gamma(\alpha m +1)} \left| \dfrac{(t^{\alpha}+h)^{m}-t^{\alpha m}}{h} - m t^{\alpha (m-1)} \right| \xrightarrow[h \to 0]{} 0, \end{split} \end{equation} $$(15)

where in the limit we have used that the deterministic function v^m is differentiable at v = t^α and that ‖b0‖2<+∞$\left\| b_0\right\|_2\lt+\infty$ since b₀ is a 2-RV. As a consequence, X_m,1(v) is m.s. differentiable at v = t^α being X′m,1(tα)=mλmb0Γ(αm+1)tα(m−1) $X_{m,1}'(t^\alpha)=\frac{m \lambda^m b_0}{\Gamma(\alpha m +1)} t^{\alpha (m-1)}$ its m.s. derivative.

Secondly, using a similar reasoning it can be checked that X′_m,1(v) = mλ^mb₀v^m−1/Γ(α m + 1) is mean square continuous at v = t^α.

Thirdly, let us show that the random power series ∑m≥0Xm,1(v)=∑m≥0Xm,1vm$\sum_{m\geq 0} X_{m,1}(v)=\sum_{m\geq 0} X_{m,1}v^m$ is mean square convergent for each v: 0 < v ≤ T^α. To this end, we first majorize this series

‖Xm,1‖2tm=‖λmb0Γ(αm+1)‖2tm≤|λ|m‖b0‖2tmΓ(αm+1):=δm(t). $$ \begin{equation} \begin{split} \left\| X_{m,1} \right\|_2 t^m = \left\| \dfrac{\lambda^m b_0}{\Gamma(\alpha m +1)} \right\|_2 t^m\leq |\lambda|^m \left\| b_0 \right\|_2 \dfrac{t^m}{\Gamma(\alpha m +1)}:=\delta_m(t). \end{split} \end{equation} $$(16)

Then, we check that series with general term δ_m(t) is convergent by applying the D’Alembert test together with the generalized Stirling formula, Γ(x+1)≈xxe−x2πx $\Gamma(x+1) \approx x^x e^{-x} \sqrt{2 \pi x}$ , x → + ∞. This leads

limm→+∞δm+1(t)δm(t)=|λ|(limm→+∞1(α(m+1))αmm+1)t=0. $$ \begin{equation} \begin{split} \lim_{m\to +\infty}\dfrac{\delta_{m+1}(t)}{\delta_m(t)} &= |\lambda| \left( \lim_{m\to +\infty} \dfrac{1}{(\alpha (m+1))^{\alpha}} \sqrt{\dfrac{m}{m+1}} \right) t=0. \end{split} \end{equation} $$(17)

Finally, the mean square uniform convergence of ∑m=0∞Xm,1(t) $\sum_{m=0}^\infty X_{m,1} (t)$ can be established using similar arguments to the ones shown earlier. In this way it is justified that the random power series ∑m≥0Xm,1vm$\sum_{m\geq 0}X_{m,1}v^m$ is mean square differentiable at v = t^α. Analogously, it can be proved that the second power series ∑m≥1Xm,2tm$\sum_{m\geq 1}X_{m,2}t^m$ in (13) is mean square differentiable at v = t^α. As a consequence, the random power series X(t) is mean square differentiable at v = t^α. Then, condition C1 is satisfied. Based upon similar arguments, it can be shown that X(v) also satisfies conditions C2 and C3.

In this section we will construct the solution of the random IVP (2) when 1<α≤2 $1 \lt \alpha \leq 2$ applying similar arguments to ones exhibited in previous section. Thus we seek the solution SP in the following form

Y(t)=Y1(t)+Y2(t), $$ \begin{equation*} Y(t)=Y_1(t)+Y_2(t), \end{equation*} $$

where

Y1(t)=∑m≥0Xmtαm,Y2(t)=∑m≥0Ymtαm+1. $$ \begin{equation} Y_1(t)=\sum_{m \geq 0}X_m t^{\alpha m}, \qquad Y_2(t)= \sum_{m \geq 0}Y_m t^{\alpha m + 1}. \end{equation} $$(18)

In order to obtain the expression of the Caputo derivative of order α to Y₁(t), we define Y^1(t)=∑m≥0Xmtm, $\hat{Y}_1(t)= \sum_{m \geq 0}X_m t^{m}$ , hence Y1(t)=Y^1(tα) $Y_1(t)=\hat{Y}_1(t^\alpha)$ . According to (3), the random mean square Caputo derivative is given by

(CD0αY1)(t)=(CD0αY^1)(tα)=(J2−αZ)(t),1<α≤2, $$ \begin{equation*} \left( ^C D_0^\alpha Y_1 \right)(t)=\left( ^C D_0^\alpha \hat{Y}_1 \right)(t^\alpha)= \left( J^{2-\alpha} Z \right)(t),\qquad 1 \lt \alpha \leq 2, \end{equation*} $$

where Z(t)=Y^″1'(tα) $Z(t)=\hat{Y}_1^{''}(t^\alpha)$ . To compute Z(t) we will apply twice Theorem 2, assuming that hypotheses i)-iii) of Theorem 2 hold for both Y^1(t) $\hat{Y}_1(t)$ and Y^′1(t) $\hat{Y}_1^{'}(t)$ . In that case

Z(t)=[(Y^1(tα))′]′=[αtα−1Y^1′(v)|v=tα]′=α(α−1)tα−2Y^(v)|v=tα+αtα−1αtα−1Y^1′​′(v)|v=tα=α(α−1)tα−2Y^1(v)|v=tα+α2t2α−2Y^1′​′(v)|v=tα=α(α−1)∑m=0∞(m+1)Xm+1tα(m+1)−2+α2∑m=0∞(m+2)(m+1)Xm+2tα(m+2)−2. $$ \begin{equation*} \begin{aligned} Z(t)&=\left [\left ( \hat{Y}_1(t^\alpha) \right)' \right]' = \left [ \alpha t^{\alpha-1} {\hat{Y}_1^{'}(v)}\Bigl\lvert_{v=t^\alpha}\right]'\\ & = \alpha (\alpha-1) t^{\alpha-2} {\hat{Y}_1^{'}(v)}\Bigl\lvert_{v=t^\alpha} + \alpha t^{\alpha-1} \alpha t^{\alpha-1} {\hat{Y}_1^{''}(v)}\Bigl\lvert_{v=t^\alpha} \\ & = \alpha (\alpha-1) t^{\alpha-2} {\hat{Y}_1^{'}(v)}\Bigl\lvert_{v=t^\alpha}+ \alpha^2 t^{2 \alpha-2} {\hat{Y}_1^{''}(v)}\Bigl\lvert_{v=t^\alpha}\\ & = \alpha (\alpha-1) \sum_{m=0}^\infty (m+1) X_{m+1} t^{\alpha (m+1) -2} + \alpha^2 \sum_{m=0}^\infty (m+2) (m+1) X_{m+2} t^{\alpha (m+2)-2}. \end{aligned} \end{equation*} $$

Observe that, we have applied Property (4.126) of [8, p.96] to compute the mean square derivative of the product of a deterministic function (α t^α−1) and a mean square 2-SP (Y^1′(tα)) $(\hat{Y}_1^{'}(t^{\alpha}))$ . Moreover, assuming that the random power series ∑m=0∞(m+1)Xm+1tα(m+1)−2 $\sum_{m=0}^\infty (m+1) X_{m+1} t^{\alpha (m+1) -2}$ and ∑m=0∞(m+2)(m+1)Xm+2tα(m+2)−2 $\sum_{m=0}^\infty (m+2)(m+1) X_{m+2} t^{\alpha (m+2)-2}$ are mean square convergent, we can obtain the random mean square Caputo derivative as follows

(CD0∞Y1)(t)=(J2−αZ)(t)=J2−α(α(α−1)∑m=0∞(m+1)Xm+1tα(m+1)−2+α2∑m=0∞(m+2)(m+1)Xm+2tα(m+2)−2)=α(α−1)∑m=0∞(m+1)Xm+1J2−α(tα(m+1)−2)+α2∑m=0∞(m+2)(m+1)Xm+2J2−α(tα(m+2)−2)=∑m≥0Γ(α(m+1)+1)Γ(am+1)Xm+1tαm,$$\begin{equation*} \begin{aligned} \left( ^C D_0^\alpha Y_1 \right)(t)&= \left( J^{2-\alpha} Z \right)(t)\\ &= J^{2-\alpha} \left( \alpha (\alpha-1) \sum_{m=0}^\infty (m+1) X_{m+1} t^{\alpha (m+1) -2} + \alpha^2 \sum_{m=0}^\infty (m+2)(m+1) X_{m+2} t^{\alpha (m+2)-2} \right)\\ &= \alpha (\alpha-1) \sum_{m=0}^\infty (m+1) X_{m+1} \, J^{2-\alpha}\left( t^{\alpha (m+1) -2} \right) + \alpha^2 \sum_{m=0}^\infty (m+2)(m+1) X_{m+2} \, J^{2-\alpha} \left (t^{\alpha (m+2)-2} \right)\\ &= \sum_{m \geq 0} \frac{\Gamma(\alpha (m+1)+1)}{\Gamma(\alpha m +1)} X_{m+1}t^{ \alpha m}, \end{aligned} \end{equation*}$$

where in the last inequality we have simplified and used the reproductive property of gamma function, Γ(γ+1) = γ Γ(γ), γ > 0.

The next step is to compute the random mean square Caputo derivative of Y₂(t). Note that by (3) one gets

(CD0αY2)(t)=(J2−αY2″)(t)=(J2−α(Y2′)′)(t)=(CD0α−1Y2′)(t). $$ \begin{equation*} \left( ^C D_0^\alpha Y_2 \right)(t)= \left( J^{2-\alpha} Y_2'' \right) (t) = \left( J^{2-\alpha} (Y_2')' \right) (t) = \left( ^C D_0^{\alpha-1} Y_2' \right)(t). \end{equation*} $$

As 1<α≤2 $1 \lt \alpha \leq 2$ , and Y2′(t)=∑m≥0(am+1)Ymtam$Y_2'(t)= \sum_{m \geq 0} (\alpha m +1) Y_m t^{\alpha m}$, we can recast α^=α−1∈(0,1] $\hat{\alpha}=\alpha-1\in (0,1]$ , Y^m=(αm+1)Ym $\hat{Y}_m = (\alpha m +1) Y_m$ and to compute the random mean square Caputo derivative of order α^ $\hat{\alpha}$ of ∑m≥0Y^mtam $\sum_{m \geq 0} \hat{Y}_mt^{\alpha m}$ using the same method used in (18)–(12). This leads

(CD0αY2)(t)=∑m≥0Ym+1Γ(α(m+1)+2)Γ(αm+2)tαm+1. $$ \begin{equation*} \left( ^C D_0^\alpha Y_2 \right)(t)=\sum_{m \geq 0} Y_{m+1} \frac{\Gamma(\alpha (m+1)+2)}{\Gamma(\alpha m +2)} t^{\alpha m+1}. \end{equation*} $$

Once we have obtained the mean square Caputo derivative of both series in (18), we need to compute their coefficients X_m and Y_m taking into account the initial conditions Y (0) = X₀ = b₀ and Y′(0) = Y₀ = b₁. After handling the corresponding recurrence relationships, one gets

Xm=λmb0+λm−1cΓ(αm+1),Ym=λmb1Γ(αm+2),m≥1. $$ \begin{equation*} X_m = \frac{\lambda ^m b_0+\lambda^{m-1} c}{\Gamma(\alpha m +1)}, \qquad Y_m= \frac{\lambda^m b_1}{\Gamma(\alpha m +2)}, \qquad m \geq 1. \end{equation*} $$

Therefore, a candidate solution SP to the random IVP (2) with 1<α≤2 $1 \lt \alpha \leq 2$ is given by

Y(t)=∑m≥0Xm,1tαm+∑m≥1Xm,2tαm+∑m≥0Ymtαm+1, $$ \begin{equation*} Y(t)= \sum_{m \geq 0}X_{m,1} t^{\alpha m} + \sum_{m \geq 1}X_{m,2} t^{\alpha m} + \sum_{m \geq 0}Y_m t^{\alpha m + 1}, \end{equation*} $$

where

Xm,1=λmb0Γ(αm+1),Xm,2=λm−1cΓ(αm+1),Ym=λmb1Γ(αm+2). $$ \begin{equation*} X_{m,1}= \frac{\lambda ^m b_0}{\Gamma(\alpha m +1)}, \qquad X_{m,2}= \frac{\lambda ^{m-1} c}{\Gamma(\alpha m +1)}, \qquad Y_m = \frac{\lambda^m b_1}{\Gamma(\alpha m +2)}. \end{equation*} $$

An important feature when solving random differential equations is that the main goal is not only to compute the solution SP but its statistical functions. This section is addressed to compute the mean, the variance, the covariance and the cross-covariance functions of the solution SP to random IVP (2) when 0 ≤ α < 1. The method to compute these statistical functions in the case that 0 < α ≤ 2. is analogous and we then skip in the subsequent development. To this end, the following property that will play a key role later. At this point, it is important to stress that that crucial property holds for mean square convergence being not true when another kind of stochastic convergences are considered.

In order to construct the approximations of the mean and the variance, it is convenient to introduce the truncation of order M of the solution SP, that is,

YM(t)=∑m=0MXm,1tαm+∑m=1MXm,2tαm=∑m=0Mλmb0Γ(αm+1)tαm+∑m=1Mλm−1cΓ(αm+1)tαm. $$ \begin{equation} \begin{aligned} Y_M(t) = \sum_{m= 0}^M X_{m,1} t^{\alpha m} + \sum_{m= 1}^M X_{m,2} t^{\alpha m} = \sum_{m= 0}^M \dfrac{\lambda^m b_0}{\Gamma(\alpha m +1)} t^{\alpha m} + \sum_{m= 1}^M \dfrac{\lambda^{m-1} c}{\Gamma(\alpha m +1)} t^{\alpha m}. \end{aligned} \end{equation} $$(20)

By applying the the expectation operator on the previous expressions, one gets

E[YM(t)]=E[b0]∑m=0MλmΓ(αm+1)tαm+E[c]∑m=1Mλm−1Γ(αm+1)tαm. $$ \begin{equation} \begin{split} \mathbb{E}[Y_M(t)] = \,\, \mathbb{E}[b_0]\sum_{m= 0}^M \dfrac{\lambda^m}{\Gamma(\alpha m +1)} t^{\alpha m} + \mathbb{E}[ c]\sum_{m= 1}^M \dfrac{\lambda^{m-1}}{\Gamma(\alpha m +1)} t^{\alpha m}. \end{split} \end{equation} $$(21)

In order to obtain the expression of the approximation of the variance function, V[YM(t)] $\mathbb{V}[Y_M(t)]$ , let us first recall that

V[YM(t)]=E[(YM(t))2]−(E[YM(t)])2, $$ \begin{equation} \mathbb{V}[Y_M(t)]=\mathbb{E}[(Y_M(t))^2]-\left( \mathbb{E}[Y_M(t)]\right)^2, \end{equation} $$(22)

then it is enough to compute E[(YM(t))2] $\mathbb{E}[(Y_M(t))^2]$ . To do that, we will take advantage of statistical independence of RVs b₀ and c. Henceforth,

E[YM(t)2]=E[(∑m=0Mλmb0Γ(am+1)tαm+∑m=1Mλm−1cΓ(αm+1)tαm)2]=E[(b0)2]∑m=0Mλ2m(Γ(αm+1))2t2αm+2E[(b0)2]∑m=1M∑n=0m−1λm+nΓ(αm+1)Γ(αn+1)tα(m+n)E[c2]∑m=1Mλ2(m−1)(Γ(αm+1))2t2αm+2E[c2]∑m=2M∑n−1m−1λm+n−2Γ(αm+1)Γ(αn+1)tα(m+n)+2E[b0]E[c]∑m=0M∑n=1Mλm+n−1Γ(αm+1)Γ(αn+1)tα(m+n). $$ \begin{equation}\label{momento_orden_2_Y} \begin{split} \mathbb{E}\left[(Y_M(t))^2\right] =& \mathbb{E}\left[ \left( \displaystyle \sum_{m= 0}^M \dfrac{\lambda^m b_0}{\Gamma(\alpha m +1)} t^{\alpha m} + \displaystyle \sum_{m= 1}^M \dfrac{\lambda^{m-1} c}{\Gamma(\alpha m +1)} t^{\alpha m} \right)^2 \right]\\ = & \displaystyle\mathbb{E}\left[(b_0)^2\right] \sum_{m= 0}^M \dfrac{\lambda^{2m} }{(\Gamma(\alpha m +1))^2} t^{2\alpha m}\\ & + \displaystyle 2\mathbb{E}\left[(b_0)^2\right]\sum_{m= 1}^M \sum_{n= 0}^{m-1} \dfrac{\lambda^{m+n} }{\Gamma(\alpha m +1)\Gamma(\alpha n +1)} t^{\alpha (m+n)} \\ &+ \displaystyle \mathbb{E}\left[c^2\right] \sum_{m= 1}^M \dfrac{\lambda^{2(m-1)} }{(\Gamma(\alpha m +1))^2} t^{2\alpha m}\\ & + \displaystyle 2 \mathbb{E}\left[c^2\right] \sum_{m= 2}^M \sum_{n= 1}^{m-1} \dfrac{\lambda^{m+n-2} }{\Gamma(\alpha m +1)\Gamma(\alpha n +1)} t^{\alpha (m+n)} \\ & +\displaystyle2\mathbb{E}[b_0]\mathbb{E}[c] \sum_{m= 0}^M \sum_{n= 1}^M \dfrac{\lambda^{m+n-1}}{\Gamma(\alpha m +1)\Gamma(\alpha n +1)} t^{\alpha (m+n)}. \end{split} \end{equation} $$(23)

In order to compute an approximation of the cross-covariance function of the solution SP, let us consider N, M ≥ 1 and t, s > 0. Using the properties of covariance, this approximation can be represented by the following expression

where

Observe that ℂYM,YN(t,t) $\mathbb{C}_{Y_M,Y_N}(t,t)$ corresponds to the covariance of the approximations of order M and N of the solution SP at the time instant t while ℂYM,YM(t,s) $\mathbb{C}_{Y_M,Y_M}(t,s)$ gives the covariance of the approximation of order M at the two time instants t and s of the solution SP.

In this section we illustrate the theoretical results established by means of several numerical examples. Let us consider the random fractional IVP (2), where 0<α≤1 $0 \lt \alpha \leq 1$ , λ=34 $\lambda=\frac{3}{4}$ , b₀ and c are 2-RVs such that

E[b0]=E[c]=12,V[b0]=V[c]=14. $$ \begin{equation} \mathbb{E}[b_0]=\mathbb{E}[c]=\frac{1}{2}, \quad \mathbb{V}[b_0]=\mathbb{V}[c]=\frac{1}{4}. \end{equation} $$(24)

Figure 1 shows the approximations of the mean and the standard deviation of the solution SP to the random IVP (2) with α = 0.7. In both plots we have carried out computations taking different orders of truncations M = 6, 7, 8, 9, 10.

Figure 1

To illustrate the role of λ parameter in the velocity of the convergence of approximations, let us now consider that λ parameter is greater than 1, for example λ=54 $\lambda=\frac{5}{4}$ . In that case, the approximations of the mean and standard deviation converge slower than in the previous case. In Figure 2, we show the approximations considering higher truncation order than in the previous case, namely M = 10, 12, 14, 16, 18.

Figure 2

On the one hand, if we assume that the mean of b₀ and c are negative, the mean of the solution SP decreases as t increases, while the standard deviation of the solution SP increases. In Figure 3, the mean and standard deviation have been represented taking E[b0]=E[c]=−1 $\mathbb{E}[b_0]=\mathbb{E}[c]=-1$ .

Figure 3

In Figures 1–3 we observe that approximations of the mean and the standard deviation improve as M increases although the error is greater as we depart from the origin (t = 0) where the IVP is stated. These results are consistent with Proposition 4.

On the other hand is important to illustrate the behaviour of the solution SP when the fractional order, α, of the derivative changes. In Figure 4, we have fixed the order of truncation, M = 20, and we show the trajectories of the expectation and standard deviation to the solution SP for different orders of the fractional derivative. Note that the case α = 0.99 is very close to the classical derivative (α = 1).

Figure 4

In this contribution we have provided a probabilistic study of the randomized deterministic fractional linear differential equation. The study has been carried out using some results belonging to random mean square calculus that have been previously established in the extant literature. We have constructed a mean square convergent power generalized series solution stochastic process by means of a random Fröbenius method. Furthermore, we have constructed approximations for both the mean and the variance of the solution stochastic process. Our numerical experiments are in agreement with the theoretical results. This paper can stimulate the extension of well-known results of deterministic fractional differential equations to the random framework using mean square calculus.