Solving Second-Order Linear Differential Equations with Random Analytic Coefficients about Regular-Singular Points

Abstract

In this contribution, we construct approximations for the density associated with the solution of second-order linear differential equations whose coefficients are analytic stochastic processes about regular-singular points. Our analysis is based on the combination of a random Fröbenius technique together with the random variable transformation technique assuming mild probabilistic conditions on the initial conditions and coefficients. The new results complete the ones recently established by the authors for the same class of stochastic differential equations, but about regular points. In this way, this new contribution allows us to study, for example, the important randomized Bessel differential equation.

1. Introduction and Motivation

As a main difference with respect to deterministic (or classical) differential equations, solving a random differential equation (RDE) does not only consist of determining, exactly or approximately, its solution stochastic process (SP), say $X\left(t\right)$, but also of computing its main statistical properties such as the mean and the variance. Even more, the approximation of the first probability density function (1-PDF), say $f_1(x,t)$, of the solution SP is a more ambitious and useful goal since from it, one can calculate, via its integration, the mean and the variance of $X\left(t\right)$,

$${\displaystyle {\mu }_X\left(t\right)=\mathbb{E}\left[X\left(t\right)\right]={\int }_{-\infty }^{\infty }x{f}_1(x,t)\mathrm{d}x,\mathbb{V}\left[X\left(t\right)\right]={\int }_{-\infty }^{\infty }{(x-{\mu }_X\left(t\right))}^2{f}_1(x,t)\mathrm{d}x,}$$

(1)

as well as its higher one-dimensional statistical moments (such as asymmetry, kurtosis, etc.),

$$\mathbb{E}\left[{\left(X\left(t\right)\right)}^k\right]={\int }_{\mathbb{R}}{x}^k{f}_1(x,t)\mathrm{d}x,k=1,2,3,\dots ,$$

provided they exist. Moreover, given a fixed time instant $\widehat{t}$, the 1-PDF permits calculating key information as the probability that the solution varies in any specific set of interest, say $[a,b]$, just via its integration,

$$\mathbb{P}\left[\left\{\omega \in \mathsf{\Omega }:a\le x(\widehat{t},\omega )\le b\right\}\right]={\int }_a^b{f}_1(x,\widehat{t})\mathrm{d}x,a,b\in \mathbb{R}.$$

The computation of the 1-PDF in the setting of RDEs and their applications is currently a cutting-edge topic for which a number of interesting advances have been achieved in the recent literature. Here, we highlight the following contributions that are strongly related to the goal of the present paper [1,2,3,4]. The main aim of this paper is to contribute to advance the field of RDEs by extending important classical results to the stochastic context. Specifically, we address the problem of constructing reliable approximations for the 1-PDF of the solution SP to second-order linear differential equations whose coefficients are analytic SPs depending on a random variable (RV) about a regular-singular point and whose initial conditions (ICs) are also RVs. Therefore, we are dealing with what is often called RDEs having a finite degree of randomness [5] (p. 37). In this manner, we complete the analysis performed in our previous contribution [6], where we dealt with the ordinary point case. It is important to emphasize that the analysis of second-order linear random differential equations has been performed in previous contributions using the random mean square calculus [5] (Ch. 4) (see [7,8,9,10,11,12] for the Airy, Hermite, Legendre, Laguerre, Chebyshev, and Bessel equations, respectively) and using other approaches like the random homotopy method [13,14], the Adomian method [15], the differential transformation method [16], the variational iteration method [17], etc. However, in all these contributions, only approximations for the two first statistical moments (mean and variance) of the solution were calculated. In contrast, in [6] and in the present paper, we deal with the general form of the aforementioned second-order linear RDEs, and furthermore, we provide additional key information of the solution via the computation of approximations for the 1-PDF that, as has been previously indicated, permits calculating not only the mean and the variance, but also higher moments of the solution, as well as further relevant information as the probability that the solution lies in specific sets of interest. To the best of our knowledge, these results for the 1-PDF of second-order linear RDEs about regular-singular points are new, and then, they contribute to advance the setting of this important class of RDEs. At this point, is important to underscore the main differences between our previous contribution [6] and the present paper. In [6], we provided a comprehensive study of second-order linear differential equations with random analytic coefficients about ordinary points. Now, we propose a natural continuation of [6] by extending the analysis for the same class of differential equations, but about singular-regular points, whose mathematical nature is completely distinct. Our aim is to complete the stochastic analysis for this important type of differential equation inspired by the well-known extension of the deterministic setting, the first one dealing with ordinary points and secondly with singular-regular points, by applying the Fröbenius theorem. It is important to point out that apparently, the problems have a strong similarity, but they are applied to differential equations of a different nature. In this regard, this new contribution allows us to study, for example, the important randomized Bessel differential equation, which does not fall within the mathematical setting of the previous contribution [6].

For the sake of completeness, below, we summarize the main results about the deterministic theory of second-order linear differential equations about ordinary and regular-singular points. These results will be very useful for the subsequent development. Let us then consider the second-order linear differential equation:

$$a_0\left(t\right)\ddot{x}\left(t\right)+a_1\left(t\right)\dot{x}\left(t\right)+a_2\left(t\right)x\left(t\right)=0,$$

(2)

where coefficients $a_i\left(t\right)$, $i=0,1,2$ are analytic functions at a certain point, say $t_0$, i.e., they admit convergent Taylor series expansions about $t_0$ (in practice, these coefficients are often polynomials, which are analytic everywhere). To study when this equation admits a power series solution centered at the point $t_0$, say $x\left(t\right)=x_0+x_1(t-t_0)+x_2{(t-t_0)}^2+\cdots$ (where the coefficients $x_i$, $i=0,1,2,\dots$ must be determined so that this series satisfies the differential equation), it is convenient to recast Equation (2) in its standard or canonical form:

$$\ddot{x}\left(t\right)+p\left(t\right)\dot{x}\left(t\right)+q\left(t\right)x\left(t\right)=0,$$

(3)

where $p\left(t\right)=a_1\left(t\right)/a_0\left(t\right)$ and $q\left(t\right)=a_2\left(t\right)/a_0\left(t\right)$. The key question is how must we pick the center of the previous expansion, $t_0$, because this choice fully determines the region of convergence of the power series. To this end, $t_0$ is classified as an ordinary or a singular point. Specifically, $t_0$ is called an ordinary point of the differential Equation (3) (or equivalently (2)) if coefficients $p\left(t\right)$ and $q\left(t\right)$ are both analytic at $t_0$, otherwise $t_0$ is termed a singular point. Recall that the quotient of analytic functions is also an analytic function provided the denominator is distinct from zero. Therefore, if $a_0\left(t_0\right)\ne 0$ in (2) (and using that their coefficients are analytic functions about $t_0$), then $t_0$ is an ordinary point. In the particular case that $a_i\left(t\right)$, $i=0,1,2$ are polynomials without common factors, then $t_0$ is an ordinary point of (2) if and only if $a_0\left(t_0\right)\ne 0$. It is well known that when $t_0$ is an ordinary point of differential Equation (2), then its general solution can be expressed as $x\left(t\right)={\sum }_{n\ge 0}{x}_n{(t-t_0)}^n={\eta }_1{x}_1\left(t\right)+{\eta }_2{x}_2\left(t\right)$, where ${\eta }_1,{\eta }_2$ are free constants and $x_1\left(t\right)$ and $x_2\left(t\right)$ are two linearly independent series solutions of (2) centered at the point $t_0$. These series are also analytic about $t_0$, i.e., they converge in a certain common interval centered at $t_0$: $0\le |t-t_0|<\rho$ ($\rho =min({\rho }_1,{\rho }_2)$ being ${\rho }_i$ the radius of convergence of $x_i\left(t\right)$, $i=1,2$, respectively). In the important case that coefficients $a_i\left(t\right)$, $i=0,1,2$ are polynomials, the radius of convergence $\rho >0$ of both series is at least as great as the distance from $t_0$ to the nearest root of $a_0\left(t\right)=0$. As a consequence, if the leading coefficient $a_0\left(t\right)$ is constant, then a power series solution expanded about any point can be found, and this power series will converge on the whole real line. In [6], we solved, in the probabilistic sense previously explained, the randomization of Equation (3) by assuming that coefficients $p\left(t\right)=p(t;A)$ and $q\left(t\right)=q(t;A)$ depend on a common RV, denoted by A, together with two random ICs fixed at the ordinary point $t_0$, namely $X\left(t_0\right)=Y_0$ and $\dot{X}\left(t_0\right)=Y_1$.

The study of Equation (3) (or equivalently, (2)) about a singular point requires further distinguishing between regular-singular points and irregular-singular points. In this paper, we shall deal with the former case, which happens when $p\left(t\right)$ approaches infinity no more rapidly than $1/(t-t_0)$ and $q\left(t\right)$ no more rapidly than $1/{(t-t_0)}^2$, as $t\to t_0$. In other words, $p\left(t\right)$ and $q\left(t\right)$ have only weak singularities at $t_0$, i.e., writing Equation (3) in the form:

$$\ddot{x}\left(t\right)+\frac{\widehat{p}\left(t\right)}{t-t_0}\dot{x}\left(t\right)+\frac{\widehat{q}\left(t\right)}{{(t-t_0)}^2}x\left(t\right)=0,$$

(4)

where:

$$\widehat{p}\left(t\right)=(t-t_0)p\left(t\right)\mathrm{and}\widehat{q}\left(t\right)={(t-t_0)}^{2}q\left(t\right),$$

(5)

then $\widehat{p}\left(t\right)$ and $\widehat{q}\left(t\right)$ are analytic about $t_0$. Otherwise, the point $t_0$ is called an irregular-singular point. In the case that $\widehat{p}\left(t\right)$ and/or $\widehat{q}\left(t\right)$ defined in (4) become indeterminate forms at $t_0$, the situation is determined by the limits:

$${\widehat{p}}_0=\widehat{p}\left(t_0\right)=\underset{t\to t_0}{lim}\widehat{p}\left(t\right)=\underset{t\to t_0}{lim}(t-t_0)p\left(t\right)\mathrm{and}{\widehat{q}}_0=\widehat{q}\left(t_0\right)=\underset{t\to t_0}{lim}\widehat{q}\left(t\right)=\underset{t\to t_0}{lim}{(t-t_0)}^{2}q\left(t\right).$$

(6)

If ${\widehat{p}}_0=0={\widehat{q}}_0$, then $t=t_0$ may be an ordinary point of the differential equation ${(t-t_0)}^2\ddot{x}\left(t\right)+(t-t_0)\widehat{p}\left(t\right)\dot{x}\left(t\right)+\widehat{q}\left(t\right)x\left(t\right)=0$ (or equivalently, dividing it by $t-t_0$ and taking into account (5), of the differential Equation (3), i.e., $\ddot{x}\left(t\right)+p\left(t\right)\dot{x}\left(t\right)+q\left(t\right)x\left(t\right)=0$). Otherwise, if both limits in (6) exist and are finite (and distinct form zero), then $t=t_0$ is a regular-singular point, while if either limit fails to exist or is infinite, then $t=t_0$ is an irregular-singular point. As has been previously underlined, the most common case in applications, in dealing with differential equations of the form (4), is when $\widehat{p}\left(t\right)$ and $\widehat{q}\left(t\right)$ are both polynomials. In such a case, ${\widehat{p}}_0$ and ${\widehat{q}}_0$ are simply the coefficients of the terms ${(t-t_0)}^0$ of these polynomials, if they are expressed in powers of $t-t_0$, so $t_0$ is a regular-singular point. In dealing with the case that $t_0$ is a regular-singular point, once the differential Equation (2) is written in the form (4), we are formally working in a deleted (or punctured) neighborhood of $t_0$, say $0<|t-t_0|<\rho$, $\rho >0$. The solution of (4) (equivalently of (2)) is then sought via generalized power series (often called Fröbenius series) centered at the point $t_0$, $x\left(t\right)=|t-t_0{|}^r{\sum }_{n\ge 0}{x}_n{(t-t_0)}^n$, $0<|t-t_0|<\rho$, being $x_0\ne 0$ and where r is a (real or complex) value to be determined, as well as coefficients $x_n$, $n=0,1,2,\dots ,$ by imposing that this series satisfies the differential Equation (4) (equivalently of (2)). This leads to the fact that parameter r must satisfy the following quadratic equation:

$$r(r-1)+{\widehat{p}}_{0}r+{\widehat{q}}_0=0.$$

(7)

This equation is called the indicial equation of differential Equation (8), and its two roots, $r_1$ and $r_2$ (possibly equal), are called the exponents of the differential equation at the regular-singular point $t_0$,

$${(t-t_0)}^2\ddot{x}\left(t\right)+(t-t_0)\widehat{p}\left(t\right)\dot{x}\left(t\right)+\widehat{q}\left(t\right)x\left(t\right)=0,$$

(8)

which is obtained after multiplying (4) by ${(t-t_0)}^2$. The full solution of differential Equation (8) (or equivalently, (2)) is given in the following theorem in terms of the nature of $r_1$ and $r_2$.

Theorem 1.

(Fröbenius method) [18] (p. 240). Let us consider the differential Equation (2) whose coefficients are analytic about $t_0$, and assume that $\widehat{p}\left(t\right)$ and $\widehat{q}\left(t\right)$ defined in (4) and (5) are analytic about $t_0$. Let us assume that $t_0$ is a regular-singular point of the differential Equation (2). Let $r_1$ and $r_2$ be the roots of the indicial equation associated with (8) at the point $t_0$:

$$r(r-1)+{\widehat{p}}_{0}r+{\widehat{q}}_0=0,$$

(9)

where ${\widehat{p}}_0$ and ${\widehat{q}}_0$ are defined in (6). Without loss of generality, we assume that $\mathrm{Re}\left(r_1\right)\ge \mathrm{Re}\left(r_2\right)$, where $\mathrm{Re}(·)$ stands for the real part.

Then, the differential Equation (2) has two linearly independent solutions, $x_1\left(t\right)$ and $x_2\left(t\right)$, in a certain deleted neighborhood centered at $t_0$, of the following form:

(a)

If $r_1-r_2$ is not a non-negative integer (i.e., $r_1-r_2\ne 0,1,2,\cdots$),

$$x_1\left(t\right)=|t-t_0{|}^{r_1}\sum _{n=0}^{\infty }{x}_{1,n}{(t-t_0)}^{n}and{x}_2\left(t\right)={|t-t_0|}^{r_2}\sum _{n=0}^{\infty }{x}_{2,n}{(t-t_0)}^n,$$

where $x_{1,0}\ne 0$ and $x_{2,0}\ne 0$.

(b)

If $r_1-r_2$ is a positive integer (i.e., $r_1-r_2=1,2,\dots$),

$$x_1\left(t\right)=|t-t_0{|}^{r_1}\sum _{n=0}^{\infty }{x}_{1,n}{(t-t_0)}^{n}and{x}_2\left(t\right)=|t-t_0{|}^{r_2}\sum _{n=0}^{\infty }{x}_{2,n}{(t-t_0)}^n+c{x}_1\left(t\right)log|t-t_0|,$$

where $x_{1,0}\ne 0$, $x_{2,0}\ne 0$, and c is a constant, which may or may not be distinct from zero.

(c)

If $r_1=r_2$,

$$x_1\left(t\right)=|t-t_0{|}^{r_1}\sum _{n=0}^{\infty }{x}_{1,n}{(t-t_0)}^n,and{x}_2\left(t\right)=|t-t_0{|}^{r_1+1}\sum _{n=0}^{\infty }{x}_{2,n}{(t-t_0)}^n+x_1\left(t\right)log|t-t_0|,$$

where $x_{1,0}\ne 0$.

As was previously indicated, in this contribution, the objective is to extend the analysis performed in [6] for the aforementioned randomization of differential Equation (3) (equivalently of (2)) about an ordinary point $t_0$, to the case that $t_0$ is a regular-singular point. Specifically, we will consider the following second-order linear RDE:

$$\ddot{X}\left(t\right)+p(t;A)\dot{X}\left(t\right)+q(t;A)X\left(t\right)=0,$$

(10)

together with the following random ICs:

$$X\left(t_1\right)=Y_0,\dot{X}\left(t_1\right)=Y_1,$$

(11)

fixed at an initial instant $t_1$ belonging to a certain deleted interval centered at $t_0$, ${\widehat{T}}_0=\{t:0<|t-t_0|<\rho \}$, $\rho >0$. Then, we will construct approximations of the 1-PDF, $f_1(x,t)$, of the SP, $X\left(t\right)$, to this random initial value problem (IVP) via a random Fröbenius series centered at the regular-singular point $t_0$. In (10), $p(t;A)$ and $q(t;A)$ are SPs satisfying certain hypotheses, which will be stated later, which depend on the common RV denoted by A. Here, we assume that we are working on a probability space $(\mathsf{\Omega },\mathcal{F},\mathbb{P})$, which is complete, to which the real RVs A, $Y_0$, and $Y_1$ belong. Furthermore, we assume that these three RVs have a joint density, i.e., the random vector $(A,Y_0,Y_1)$ is absolutely continuous.

As usual, observe that in (10) and (11), we are writing RVs with capital letters, for instance, $Z:\mathsf{\Omega }⟶{\mathcal{D}}_Z\subset \mathbb{R}$ (here, ${\mathcal{D}}_Z$ is referred to as the range of Z). Moreover, each realization of an RV, say Z, will be denoted by $z\left(\omega \right)$, $\omega \in \mathsf{\Omega }$ or simply $z=z\left(\omega \right)\in {\mathcal{D}}_Z$. To be as general as possible, in the subsequent analysis, A, $Y_0$, and $Y_1$ are assumed to be dependent absolutely continuous RVs, and $f_{A,Y_0,Y_1}(a,y_0,y_1)$ will denote their joint PDF. Notice that, if A, $Y_0$, and $Y_1$ are independent RVs, then their PDF can be factorized as the product of their respective marginal PDFs, i.e.,

$$f_{A,Y_0,Y_1}(a,y_0,y_1)=f_A\left(a\right)f_{Y_0}\left(y_0\right)f_{Y_1}\left(y_1\right).$$

(12)

Based on our previous contribution [6], together with Theorem 1 and further reasons that will be apparent later, hereinafter, we will assume the following hypotheses:

H0: A is a bounded RV, i.e., $\exists M_A>0:|a\left(\omega \right)|<M_A,\forall \omega \in \mathsf{\Omega }.$

H1: $f_{A,Y_0,Y_1}(a,y_0,y_1)$ is continuous in the second component and bounded, i.e.,

$$\begin{array}{c}\exists M_f>0:\left|f_{A,Y_0,Y_1}(a,y_0,y_1)\right|<M_f,\forall (a,y_0,y_1)\in {\mathcal{D}}_{A,Y_0,Y_1}.\end{array}$$

In addition, we will assume that the SPs $\widehat{p}(t;A)=(t-t_0)p(t;A)$ and $\widehat{q}(t;A)={(t-t_0)}^{2}q(t;A)$ are analytic about $(t_0;a_0\left(\omega \right))$ for every $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$, i.e.,

H2: There exists a common neighborhood ${\mathcal{N}}_{\widehat{p},\widehat{q}}(t_0;a_0\left(\omega \right))\subset {\widehat{T}}_0\times {\mathcal{D}}_A$ where:

$$\begin{array}{c}\widehat{p}(t;A)=(t-t_0)p(t;A),\widehat{q}(t;A)={(t-t_0)}^{2}q(t;A)\mathrm{are}\mathrm{analytic}\forall a_0\left(\omega \right)\in {\mathcal{D}}_A,\omega \in \mathsf{\Omega }.\end{array}$$

Here, ${\widehat{T}}_0$ denotes a deleted interval centered at $t_0$ that has been previously defined. Recall that here, the SP $\widehat{p}(t;A)=(t-t_0)p(t;A)$ is analytic about a point $(t_0,a_0\left(\omega \right))$, $a_0\left(\omega \right)\in {\mathcal{D}}_A$, for all $\omega \in \mathsf{\Omega }$, if the deterministic function $(t-t_0)p(t;a\left(\omega \right))$ is analytic about $(t_0,a_0\left(\omega \right))$ (equivalently to the SP $\widehat{q}(t;A)={(t-t_0)}^{2}q(t;A)$). To simplify notation, we shall assume that $\widehat{p}(t;A)=(t-t_0)p(t;A)$ and $\widehat{q}(t;A)={(t-t_0)}^{2}q(t;A)$ are analytic in a common neighborhood that will be denoted by ${\mathcal{N}}_{\widehat{p},\widehat{q}}(t_0;a_0\left(\omega \right))$. In practice, this neighborhood is determined intersecting the domains of analyticity of $\widehat{p}(t;A)=(t-t_0)p(t;A)$ and $\widehat{q}(t,A)={(t-t_0)}^{2}q(t;A)$. In [5] (Th. 4.4.3), a characterization was stated of analyticity of second-order SPs (those having finite variance) in terms of the analyticity of the correlation function. Moreover, to ensure that the IVP (10)–(11) has a unique solution, both SPs $p(t;A)$ and $q(t;A)$ are assumed to satisfy all the needed conditions. (see [5] (Th. 5.1.2), for instance).

According to the Fröbenius method (see Theorem 1) and under the analyticity condition assumed in Hypothesis H2, the solution SP of RDE (10) about a regular-singular point, $t_0$, can be written as a linear combination of two uniformly convergent independent random series, $X_1(t;A)$ and $X_2(t;A)$,

$$X\left(t\right)=K_0(A,Y_0,Y_1)X_1(t;A)+K_1(A,Y_0,Y_1)X_2(t;A),$$

where $t\in {\widehat{T}}_0$, and the coefficients $K_0(A,Y_0,Y_1)$ and $K_1(A,Y_0,Y_1)$ can be obtained in terms of the random ICs given in (11), which are established at the time instant $t_1\in {\widehat{T}}_0$. In the subsequent development, we will take advantage of the above representation along with the application of the random variable transformation (RVT) technique [5] (p. 25), [6] (Th. 1), to construct the 1-PDF, $f_1^N(x,t)$, corresponding to approximations, $X^N\left(t\right)$, of the solution SP, $X\left(t\right)$, about the regular-singular point $t_0$. We shall provide sufficient conditions so that $f_1^N(x,t)$ approximate the 1-PDF, $f_1(x,t)$, of $X\left(t\right)$. The RVT method has been successfully applied to obtain exact or approximate representations of the 1-PDF of the solution SP to random ordinary/partial differential and difference equations and systems [1,2,3,4]. In particular, the RVT technique has also been applied to conduct the study for the random IVP (10) and (11) in the case that $t_0$ is an ordinary point [6]. Thus, the present contribution can be considered as a natural continuation of our previous paper [6]. This justifies that in our subsequent development, we directly refer to [6] when applying a number of technical results already stated in the above-mentioned contribution. In this manner, we particularly avoid repeating the multidimensional random variable transformation technique [6] (Th. 1), Poincare’s expansion [6] (Th. 2), as well as several results related to uniform convergence that can be inferred from classical real Analysis [6] (Prop. 1-4). For the sake of clarity, we point out that we keep identical the notation in both contributions. These results were also extensively applied in [6,19,20,21].

The paper is organized as follows. In Section 2, the 1-PDF, $f_1^N(x,t)$, of the approximate solution SP, $X^N\left(t\right)$, to the random IVP (10) and (11) is formally constructed. This function is obtained by applying the RVT method to $X^N\left(t\right)$, which follows from truncating the random generalized power series solution derived after applying the Fröbenius method stated in Theorem 1. Section 3 is devoted to rigorously proving the convergence of approximations $f_1^N(x,t)$ to the exact 1-PDF, $f_1(x,t)$, associated with the solution SP, $X\left(t\right)$. In Section 4, several illustrative examples are shown to demonstrate the usefulness of the theoretical results established in Section 2 and Section 3. Our main conclusions are drawn in Section 5.

2. Computation of the 1-PDF of the Truncated Solution SP

As was indicated in the Introduction, by the Fröbenius method, Theorem 1, the solution SP of the random IVP (10) and (11) can be written as:

$$X\left(t\right)=K_0(A,Y_0,Y_1)X_1(t;A)+K_1(A,Y_0,Y_1)X_2(t;A),$$

(13)

where $t\in {\widehat{T}}_0$ and $X_1(t;A)$ and $X_2(t;A)$ are determined taking into account the values of the random roots of the associated indicial Equation (9). Since A is assumed to be an absolutely continuous RV, Cases (b) and (c) in Theorem 1 have null probability to occur, because punctual probabilities considering continuous RVs are zero. Thus, by (a), for each $t\in {\widehat{T}}_0$, the solution SP is given by Expression (13) being:

$$X_1(t;A)=|t-t_0{|}^{r_1\left(A\right)}\sum _{n=0}^{\infty }{X}_{1,n}\left(A\right){\left(t-t_0\right)}^n\mathrm{and}{X}_2(t;A)={|t-t_0|}^{r_2\left(A\right)}\sum _{n=0}^{\infty }{X}_{2,n}\left(A\right){\left(t-t_0\right)}^n,$$

(14)

where $r_1\left(A\right)$ and $r_2\left(A\right)$ are the random roots of the random indicial Equation (9). Random coefficients $X_{1,n}\left(A\right)$ and $X_{2,n}\left(A\right)$ are recursively determined in practice. At this point, it is convenient to underline that the factors $|t-t_0{|}^{r_1\left(A\right)}$ and $|t-t_0{|}^{r_2\left(A\right)}$ are the distinctive and major difference in dealing with the analysis of RDEs about regular-singular points with respect to the case of RDEs about ordinary points presented in [6].

By imposing that Expression (13) satisfies the ICs given in (11), we obtain the following random algebraic system:

$$\begin{array}{ccccc}{Y}_0\hfill & =& X\left(t_1\right)\hfill & =& K_0(A,Y_0,Y_1)X_1(t_1;A)+K_1(A,Y_0,Y_1)X_2(t_1;A),\hfill \\ Y_1\hfill & =& \dot{X}\left(t_1\right)\hfill & =& K_0(A,Y_0,Y_1){\dot{X}}_1(t_1;A)+K_1(A,Y_0,Y_1){\dot{X}}_2(t_1;A).\hfill \end{array}$$

Solving for $K_0$ and $K_1$, one gets:

$${\displaystyle K_0(A,Y_0,Y_1)=\frac{Y_0{\dot{X}}_2(t_1;A)-Y_1{X}_2(t_1;A)}{W(t_1;A)},K_1(A,Y_0,Y_1)=\frac{-Y_0{\dot{X}}_1(t_1;A)+Y_1{X}_1(t_1;A)}{W(t_1;A)},}$$

(15)

where $W(t;A)=X_1(t;A){\dot{X}}_2(t;A)-X_2(t;A){\dot{X}}_1(t;A)$ is the Wronskian of the fundamental system $\{X_1\left(t\right),X_2\left(t\right)\}$ given by (14).

Remark 1.

Since A is a continuous RV, then:

$$\mathbb{P}\left[\{\omega \in \mathsf{\Omega }:W(t_1;a\left(\omega \right))\ne 0\}\right]=1.$$

Moreover, there exists a positive constant, $m_W$, such that:

$$0<m_W<|W(t_1;a\left(\omega \right))|,\forall \omega \in \mathsf{\Omega },$$

being this lower bound $m_W$ independent of $\omega \in \mathsf{\Omega }$ since, according to Hypothesis H0, RV A is bounded.

Hence, according to Remark 1, the functions $K_0(A,Y_0,Y_1)$ and $K_1(A,Y_0,Y_1)$, given in (15), are well defined. Summarizing, the solution SP (13) and (14) can be written in the form:

$$X\left(t\right)=Y_0{S}_1(t;A)+Y_1{S}_2(t;A),$$

(16)

where $Y_0$ and $Y_1$ are the random ICs and:

$$\begin{array}{c}{\displaystyle S_1(t;A)=\frac{1}{W(t_1;A)}\left({\dot{X}}_2(t_1;A)X_1(t;A)-{\dot{X}}_1(t_1;A)X_2(t;A)\right),}\hfill \\ {\displaystyle S_2(t;A)=\frac{1}{W(t_1;A)}\left(-X_2(t_1;A)X_1(t;A)+X_1(t_1;A)X_2(t;A)\right).}\hfill \end{array}$$

(17)

From an analytical point of view, it could be infeasible to obtain the explicit expression of both series, $X_1(t;A)$ and $X_2(t;A)$. This fact makes it more advisable to consider the approximate solution SP that comes from truncating both series:

$$X^N\left(t\right)=Y_0{S}_1^N(t;A)+Y_1{S}_2^N(t;A),$$

(18)

where:

$$\begin{array}{c}{\displaystyle S_1^N(t;A)=\frac{1}{W^N(t_1;A)}\left({\dot{X}}_2^N(t_1;A)X_1^N(t;A)-{\dot{X}}_1^N(t_1;A)X_2^N(t;A)\right),}\hfill \\ {\displaystyle S_2^N(t;A)=\frac{1}{W^N(t_1;A)}\left(-X_2^N(t_1;A)X_1^N(t;A)+X_1^N(t_1;A)X_2^N(t;A)\right).}\hfill \end{array}$$

(19)

with $W^N(t_1;A)=X_1^N(t_1;A){\dot{X}}_2^N(t_1;A)-X_2^N(t_1;A){\dot{X}}_1^N(t_1;A)$, and $X_1^N(t;A)$ and $X_2^N(t;A)$ result from the truncation of $X_1(t;A)$ and $X_2(t;A)$, defined in (14), at a common order N, i.e.,

$$X_1^N(t;A)=|t-t_0{|}^{r_1\left(A\right)}{\sum }_{n=0}^N{X}_{1,n}\left(A\right){\left(t-t_0\right)}^n\mathrm{and}{X}_2^N(t;A)={|t-t_0|}^{r_2\left(A\right)}{\sum }_{n=0}^N{X}_{2,n}\left(A\right){\left(t-t_0\right)}^n.$$

(20)

In the context of the RVT method [6] (Th. 1) and following the same methodology as in the contribution [6], we compute the 1-PDF of the truncated solution SP (18) and (19), obtaining:

$$f_1^N(x,t)={\int }_{{\mathcal{D}}_{A,Y_1}}{f}_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{S_1^N(t;a)}\right|\mathrm{d}a\mathrm{d}{y}_1,$$

(21)

for every $t\in {\widehat{T}}_0$ arbitrary.

3. Study of the Convergence

In this section, we give sufficient conditions to guarantee the convergence of the approximation, $f_1^N(x,t)$, given in (21), to the exact function $f_1(x,t)$ when N tends to infinity, i.e., we study when:

$$\underset{N\to +\infty }{lim}{f}_1^N(x,t)=f_1(x,t),\mathrm{for}\mathrm{each}(x,t)\in \mathbb{R}\times {\widehat{T}}_0\cap [t_1,+\infty [\mathrm{fixed},$$

(22)

being:

$$f_1(x,t)={\int }_{{\mathcal{D}}_{A,Y_1}}{f}_{A,Y_0,Y_1}\left(a,{\displaystyle \frac{x-y_1{S}_2(t;a)}{S_1(t;a)}},y_1\right)\left|{\displaystyle \frac{1}{S_1(t;a)}}\right|\mathrm{d}{y}_1\mathrm{d}a,$$

(23)

where $S_1(t;a)$ and $S_2(t;a)$ are the realizations of the random series $S_1(t;A)$ and $S_2(t;A)$ defined in (17).

Now, we include here some remarks about the convergence and boundedness of the series involved in (23), which will play a key role in our subsequent developments.

Remark 2.

Expanding first the terms $|t-t_0{|}^{r_1\left(A\right)}$ and $|t-t_0{|}^{r_2\left(A\right)}$ as their corresponding series and interpreting the RV A as a parameter indexed by $\omega \in \mathsf{\Omega }$, $A=\left\{a\right(\omega ):\omega \in \mathsf{\Omega }\}$, under the hypothesis H2, Poincaré’s theorem, stated in [6](Th. 2), permits representing the solution $X\left(t\right)$ as a series of parameter $a\left(\omega \right)$, $\omega \in \mathsf{\Omega }$. Therefore, $S_1(t;A)$ and $S_2(t;A)$ do. On the other hand, taking into account the uniqueness of the solution of IVP (10)–(11), both series expansions (as powers of $t-t_0$ and as powers of $a\left(\omega \right)-a_0\left(\omega \right)$) match. Henceforth, the series $X_1(t;A)$ and $X_2(t;A)$, given by (14), are convergent in $t_0$-deleted neighborhoods ${\mathcal{N}}_{X_1}(t_0;a_0\left(\omega \right))$ and ${\mathcal{N}}_{X_2}(t_0;a_0\left(\omega \right))$, respectively, for all $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$. In addition, uniform convergence takes place in every closed set contained in ${\mathcal{N}}_{X_1}(t_0;a_0\left(\omega \right))$ and ${\mathcal{N}}_{X_2}(t_0;a_0\left(\omega \right))$. Notice that the domain of convergence, ${\mathcal{N}}_X(t_0;a_0\left(\omega \right))$, of the series solution $X\left(t\right)$ given in (13), satisfies ${\mathcal{N}}_{\widehat{p},\widehat{q}}(t_0;a_0\left(\omega \right))\subseteq {\mathcal{N}}_X(t_0;a_0\left(\omega \right))={\mathcal{N}}_{X_1}(t_0;a_0\left(\omega \right))\cap {\mathcal{N}}_{X_2}(t_0;a_0\left(\omega \right))$ for all $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$, where ${\mathcal{N}}_{\widehat{p},\widehat{q}}(t_0;a_0\left(\omega \right))$ is defined in Hypothesis H1.

Remark 3.

Functions $S_1(t;A)$ and $S_2(t;A)$ are linear combinations of series $X_1(t;A)$ and $X_2(t;A)$ (see (17)), which, by Remark 2, are uniformly convergent in every closed set ${\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$ contained in ${\mathcal{N}}_X(t_0;a_0\left(\omega \right))$. Then, $S_1(t;A)$ and $S_2(t;A)$ also converge uniformly in ${\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$.

Remark 4.

Note that for every $\omega \in \mathsf{\Omega }$:

$$S_1^N(t_1;a\left(\omega \right))=\frac{1}{W^N(t_1;a\left(\omega \right))}\left({\dot{X}}_2^N(t_1;a\left(\omega \right))X_1^N(t_1;a\left(\omega \right))-{\dot{X}}_1^N(t_1;a\left(\omega \right))X_2^N(t_1;a\left(\omega \right))\right)\equiv 1.$$

Then, according to Remark 2, there exist a certain neighborhood, ${\mathcal{N}}_X(t_1;a_0\left(\omega \right))$, and a positive constant, $m_{S_1}$, such that for all $N\ge 0$:

$$0<m_{S_1}\le |S_1^N(t;a\left(\omega \right))|,\forall (t,a\left(\omega \right))\in {\mathcal{N}}_X(t_1;a_0\left(\omega \right)),\forall a_0\left(\omega \right)\in {\mathcal{D}}_A,\omega \in \mathsf{\Omega }.$$

(24)

On the other hand, by Remark 3, $S_i(t;A)$, $i=1,2$, converge uniformly in every closed set ${\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$ contained in ${\mathcal{N}}_X(t_0;a_0\left(\omega \right))$ for all $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$. This guarantees the existence of constants $M_{S_1},M_{S_2}>0$ such that:

$$\left|S_i^N(t;a\left(\omega \right))\right|<M_{S_i},$$

(25)

$$\forall (t,a\left(\omega \right))\in {\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))\subset {\mathcal{N}}_X(t_0;a_0\left(\omega \right)),\forall a_0\left(\omega \right)\in {\mathcal{D}}_A,\omega \in \mathsf{\Omega },\forall N\ge 0,i=1,2.$$

Let $(x,t)\in \mathbb{R}\times {\widehat{T}}_0\cap [t_1,+\infty [$ be fixed; below, we establish conditions to assure the convergence stated in (22). First, let us take limits as $N\to \infty$ in Expression (21):

$$\underset{N\to +\infty }{lim}{f}_1^N(x,t)=\underset{N\to +\infty }{lim}{\int }_{{\mathcal{D}}_{A,Y_1}}{f}_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{S_1^N(t;a)}\right|\mathrm{d}a\mathrm{d}{y}_1.$$

To commute the limit as $N\to +\infty$ and the above double integral, we apply [6] (Prop. 1). We assume the following hypothesis:

H3: ${\mathcal{D}}_{A,Y_1}$ is a Lebesgue measurable set of ${\mathbb{R}}^2$ with finite measure such that:

$$\begin{array}{c}\widehat{M}=sup\left\{|y_1\left(\omega \right)|:(a\left(\omega \right),y_1\left(\omega \right))\in {\mathcal{D}}_{A,Y_1},\omega \in \mathsf{\Omega }\right\}<+\infty .\end{array}$$

(26)

Then, we shall prove that $\left\{k^N(a,y_1):N\ge 0\right\}\subset {\mathcal{L}}^1\left({\mathcal{D}}_{A,Y_1}\right)$ and $k_N(a,y_1)\underset{N\to +\infty }{\overset{\mathrm{uniformly}{\mathcal{D}}_{A,Y_1}}{\to }}k(a,y_1)$, being:

$${\displaystyle k^N(a,y_1)=f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{S_1^N(t;a)}\right|}$$

and:

$${\displaystyle k(a,y_1)=f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},y_1\right)\left|\frac{1}{S_1(t;a)}\right|.}$$

By Remark 4, considering the lower bound given in (24) for $|S_1^N(t;a)|$, one gets:

$$k^N(a,y_1)\le \frac{1}{m_{S_1}}{f}_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right).$$

Therefore, as $f_{A,Y_0,Y_1}$ is a PDF, we conclude that $k^N(a,y_1)$ is Lebesgue integrable in ${\mathcal{D}}_{A,Y_1}$, i.e., $\left\{k^N(a,y_1):N\ge 0\right\}\subset {\mathcal{L}}^1\left({\mathcal{D}}_{A,Y_1}\right)$. To prove that ${\displaystyle \underset{N\to +\infty }{lim}{k}^N(a,y_1)=k(a,y_1)}$ uniformly in ${\mathcal{D}}_{A,Y_1}$, we will apply the same argument as in the previous contribution [6]. We demonstrate that for every $ϵ>0$, there exists $N_0$ such that for all $N\ge N_0$:

$$\left|k^N(a,y_1)-k(a,y_1)\right|<ϵ,\forall (a,y_1)\in {\mathcal{D}}_{A,Y_1}.$$

Adding and subtracting the term $f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{S_1(t,a)}\right|$ and applying the triangular inequality, one gets:

$$\begin{array}{ccc}\left|k^N(a,y_1)-k(a,y_1)\right|\hfill & \le & {\displaystyle f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{|S_1^N(t,a)|}-\frac{1}{|S_1(t,a)|}\right|}\hfill \\ & & {\displaystyle +\frac{1}{|S_1(t;a)|}\left|f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)-f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},y_1\right)\right|.}\hfill \end{array}$$

Let $(x,t)\in \mathbb{R}\times {\widehat{T}}_0\cap [t_1,+\infty [$ be fixed and $N\ge 0$ be an arbitrary non-negative integer, and consider $(t,a)=(t,a\left(\omega \right))\in {\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$, for all $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$. Then, we apply [6] (Prop. 2) to $z_1=t$ fixed, $z_2=a\in {\mathcal{D}}_A$ arbitrary, $D=\{(t,a)\in {\mathbb{R}}^2:a\in {\mathcal{D}}_A\}$, $g_N(t,a)=1$, $g(t,a)=1$, $h_N(t,a)=|S_1^N(t;a)|$, and $h(t,a)=|S_1(t;a)|$. Note that, by Remark 4, $h_N(t,a)\times h(t,a)\ne 0$ for all $(t,a)\in D$. In addition, it is obvious that $g_N(t,a)$ converges uniformly to $g(t,a)$ on D. Regarding the real function $h_N(t,a)$, it converges uniformly to $h(t,a)$, by the arguments shown in Remarks 2 and 3. As $g_N(t,a)=1$, its absolute value is bounded. By Remark 4, there exists constants $m_h=m_{S_1}$ and $M_h=M_{S_1}$, such that $0<m_h<|h_N(t,a)|<M_h$, $\forall a\in {\mathcal{D}}_A$. Taking into account the uniform convergences, the bounds, and the hypothesis H1 (boundedness of the PDF $f_{A,Y_0,Y_1}(a,y_0,y_1)$), for every $ϵ>0$, there exists $N_0$ (which depends on $ϵ$) such that:

$$f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{|S_1^N(t;a)|}-\frac{1}{|S_1(t;a)|}\right|\le M_f\frac{ϵ}{2M_f}=\frac{ϵ}{2},\forall N\ge N_0,$$

independently of the values of $(a,y_1)\in {\mathcal{D}}_{A,Y_1}$.

Now, we obtain an analogous result for the second term. First, assuming Hypothesis H3, we apply [6] (Prop. 3) to $u=a$, $v=-y_1$, $c=x$, $l_N\left(a\right)=S_2^N(t;a)$, $l\left(a\right)=S_2(t;a)$, $b_N(a,y_1)=x-y_1{S}_2^N(t;a)$ and $b(a,y_1)=x-y_1{S}_2(t;a)$, with $t\in {\widehat{T}}_0\cap [t_1,+\infty [$ fixed, and such that $(t,a)=(t,a\left(\omega \right))\in {\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$, $\omega \in \mathsf{\Omega }$. By Remark 3, $l_N\left(a\right)$ converges uniformly to $l\left(a\right)$ on ${\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$. Then, given $(x,t)\in \mathbb{R}\times {\widehat{T}}_0\cap [t_1,+\infty [$ fixed, $x-y_1{S}_2^N(t;a)$ converges uniformly to $x-y_1{S}_2(t;a)$ in ${\mathcal{D}}_{A,Y_1}$. Secondly, we apply [6] (Prop. 2) taking $z_1=a$, $z_2=y_1$, $g_N(a,y_1)=x-y_1{S}_2^N(t;a)$, $g(a,y_1)=x-y_1{S}_2(t;a)$, $h_N(a,y_1)=S_1^N(t;a)$, and $h(a,y_1)=S_1(t;a)$. Note that, previously we have proven, by applying [6] (Prop. 3), that $g_N(a,y_1)$ converges uniformly to $g(a,y_1)$. In addition, by Remark 3, $h_N(a,y_1)$ converges uniformly to $h(a,y_1)$. Let $M_g=|x|+\widehat{M}{M}_{S_2}$ and $M_h=M_{S_1}$, these bounds being the ones established in Hypothesis H3 and Remark 4, then:

$$\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)}\underset{N\to +\infty }{\overset{\mathrm{uniformly}{\mathcal{D}}_{A,Y_1}}{\to }}\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},(x,t)\in \mathbb{R}\times {\widehat{T}}_0\cap [t_1,+\infty [\mathrm{fixed}.$$

The next step is to apply [6] (Prop. 4) with the following identification: $z_1=a$, $z_2=y_1$, $(a,y_1)\in {\mathcal{D}}_{A,Y_1}$ such that $(t,a)\in {\mathcal{N}}_X^{*}(t_0;a_0\left(\omega \right))$ for every $a_0\left(\omega \right)\in {\mathcal{D}}_A$, $\omega \in \mathsf{\Omega }$, ${\gamma }_N(a,y_1)=\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)}$, $\gamma (a,y_1)=\frac{x-y_1{S}_2(t;a)}{S_1(t;a)}$, $\varphi (a,y_0,y_1)=f_{A,Y_0,Y_1}(a,y_0,y_1)$, ${\psi }_N(a,y_1)=f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)$ and $\psi (a,y_1)=f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},y_1\right)$. Finally, as it was previously shown, the sequence ${\gamma }_N(a,y_1)$ is uniformly convergent to $\gamma (a,y_1)$, and according to Hypothesis H1, the mapping $\varphi$ is continuous in $y_0$. Therefore, taking into account the lower bound of $|S_1(t;a)|$ given in Remark 4 and applying [6] (Prop. 4), for every $ϵ>0$ there exists $N_0$, which depends on $ϵ$, such that

$$\frac{1}{|S_1(t;a)|}\left|f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)-f_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},y_1\right)\right|<\frac{1}{m_{S_1}}\frac{m_{S_1}ϵ}{2}=\frac{ϵ}{2},\forall N\ge N_0,$$

independently of the values of $(a,y_1)$. Hence, we proved that $k^N(a,y_1)\underset{N\to +\infty }{\overset{\mathrm{uniformly}{\mathcal{D}}_{A,Y_1}}{\to }}k(a,y_1)$, and then, by [6] (Prop. 1) and Hypothesis H3, we can commute the limit and the integral. Therefore, applying the continuity of the joint PDF $f_{A,Y_1,Y_2}(a,y_1,y_2)$ with respect to the second variable (see Hypothesis H1) and taking into account Expression (23), one gets:

$$\begin{array}{ccc}{\displaystyle \underset{N\to \infty }{lim}{f}_1^N(x,t)}\hfill & =& {\displaystyle {\int }_{{\mathcal{D}}_{A,Y_1}}\underset{N\to \infty }{lim}{f}_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{S_1^N(t;a)}\right|\mathrm{d}a\mathrm{d}{y}_1,}\hfill \\ & =& {\displaystyle {\int }_{{\mathcal{D}}_{A,Y_1}}{f}_{A,Y_0,Y_1}\left(a,\underset{N\to \infty }{lim}\frac{x-y_1{S}_2^N(t;a)}{S_1^N(t;a)},y_1\right)\left|\frac{1}{{\displaystyle \underset{N\to \infty }{lim}{S}_1^N(t;a)}}\right|\mathrm{d}a\mathrm{d}{y}_1}\hfill \\ & =& {\displaystyle {\int }_{{\mathcal{D}}_{A,Y_1}}{f}_{A,Y_0,Y_1}\left(a,\frac{x-y_1{S}_2(t;a)}{S_1(t;a)},y_1\right)\left|\frac{1}{{\displaystyle S_1(t;a)}}\right|\mathrm{d}a\mathrm{d}{y}_1=f_1(x,t).}\hfill \end{array}$$

Summarizing, we establish the following result.

Theorem 2.

Let us consider the random IVP (10) and (11), and assume that:

(i)

The RV A satisfies Hypothesis H0.

(ii)

The joint PDF, $f_{A,Y_0,Y_1}(a,y_0,y_1)$, of the random input data $(a,y_0,y_1)$ satisfies Hypothesis H1.

(iii)

The coefficients $p(t;A)$ and $q(t;A)$ satisfy Hypothesis H2.

(iv)

The ${\mathcal{D}}_{A,Y_1}$ domain of the random vector $(A,Y_1)$ satisfies Hypothesis H3.

Let $X_1^N(t,a)$ and $X_2^N(t;a)$ be the finite sum defined in (20). Then, $f_1^N(x,t)$, defined by (19)–(21), is the 1-PDF of the truncated solution SP $X^N\left(t\right)$ given by (13)–(15) to the IVP (10) and (11). Furthermore, for each $(x,t)$, 1-PDF $f_1^N(x,t)$ converges to $f_1(x,t)$ given in (17) and (23).

4. Examples

In this section, we present two examples to illustrate the theoretical results previously established. In both examples, we calculate the 1-PDF of the approximate solution SP, and from it, we also approximate the mean and the variance of the solution by taking advantage of the expressions given in (1).

Example 1.

Let us consider the random IVP:

$$\left.\begin{array}{ccc}A{t}^2\ddot{X}\left(t\right)+t(t+1)\dot{X}\left(t\right)-X\left(t\right)& =& 0,\hfill \\ X\left(1\right)& =& Y_0,\hfill \\ \dot{X}\left(1\right)& =& Y_1.\hfill \end{array}\right\}$$

(27)

Notice that $t_0=0$ is a regular-singular point since according to Equation (2), $a_0\left(t\right)=A{t}^2$, $a_1\left(t\right)=t(t+1)$ and $a_2\left(t\right)=-1$, being $a_0\left(t_0\right)=0$ and, according to (4) and (5), $p(t;A)=t(t+1)/\left(A{t}^2\right)$ and $q(t;A)=-1/\left(A{t}^2\right)$, which are not analytic about $t_0=0$, while $\widehat{p}(t;A)=(t+1)/A$ and $\widehat{q}(t;A)=-1/A$ are both analytic about $t_0=0$. We will assume that A, $Y_0$, and $Y_1$ are independent RVs with the following distributions:

• A has a uniform distribution on the interval $[1,2]$, i.e., $A\sim \mathrm{U}\left(\right[1,2\left]\right)$.
• $Y_0$ has a Gaussian RV with mean zero and standard deviation $0.1$, i.e., $Y_0\sim \mathrm{N}(0;0.1)$.
• $Y_1$ has a beta distribution with parameters $\alpha =2$ and $\beta =3$, i.e., $Y_1\sim \mathrm{Be}(2;3)$.

First, we check that Hypotheses H0–H3 hold:

• H0: Since $A\sim \mathrm{U}\left(\right[1,2\left]\right)$, it is bounded. With the notation of Hypothesis H0, we can take $M_A=2$.
• H1: Since A, $Y_0$, and $Y_1$ are assumed to be independent RVs, their joint PDF is given by the product of the PDF of each RV; see Expression (12). In this case, proving the continuity of the joint PDF with respect to the second variable, $y_0$, is equivalent to proving the continuity of the PDF of the RV $Y_0$. Since $Y_0$ has a Gaussian distribution, then it is continuous and bounded on the whole real line. As the PDFs of A and $Y_1$ are also bounded, Hypothesis H1 holds,

  $$f_A\left(a\right)f_{Y_0}\left(y_0\right)f_{Y_1}\left(y_1\right)=\frac{1}{2-1}\frac{10}{\sqrt{2\pi }}{exp}^{-50y_0^2}\frac{4!}{1!2!}{y}_1{(1-y_1)}^2\le \frac{120}{\sqrt{2\pi }},\forall (a,y_0,y_1)\in [1,2]\times \mathbb{R}\times [0,1].$$
• H2: As it was previously indicated, $p(t;A)$ and $q(t;A)$ are not analytic functions at $t_0=0$, but $\widehat{p}(t;A)$ and $\widehat{q}(t;A)$ are, since both are polynomials as a function of t for every $a=a\left(\omega \right)$, $\omega \in \mathsf{\Omega }$.
• H3: The domain of random vector $(A,Y_1)$ is:

  $${\mathcal{D}}_{A,Y_1}={\mathcal{D}}_A\times {\mathcal{D}}_{Y_1}=[1,2]\times [0,1].$$
  Hence, ${\mathcal{D}}_{A,Y_1}$ is a Lebesgue measurable set of ${\mathbb{R}}^2$ and $\widehat{M}=1<\infty$. Therefore, Hypothesis H3 is fulfilled.

In this case, the solution SP of IVP (27) is given by:

$$X\left(t\right)=Y_0{S}_1(t;A)+Y_1{S}_2(t;A),$$

(28)

where $S_i(t;A)$, $i=1,2$, are defined in (17) with:

$$X_1(t;A)=t\left(1+\sum _{n=1}^{\infty }{(-1)}^n\prod _{j=1}^n\frac{1}{(j+1)A+1}{t}^n\right)and{X}_2(t;A)=t^{-1/A}\left(\sum _{n=0}^{\infty }\frac{{(-1)}^n}{n!A^n}{t}^n\right);$$

(29)

see Appendix A for further details. In this case, we can compute, via Mathematica${}^{\circledR }$, the exact value of both power series,

$${\displaystyle X_1(t;A)=-(1+A)e^{-t/A}{\left(-\frac{t}{A}\right)}^{-1/A}\left(\Gamma \left(1+\frac{1}{A}\right)-\Gamma \left(1+\frac{1}{A},-\frac{t}{A}\right)\right)and{X}_2\left(t\right)=e^{-t/A}{t}^{-1/A},}$$

(30)

where $\Gamma \left(z\right)$ and $\Gamma (b,z)$ are the gamma and the incomplete gamma function defined, respectively, by:

$$\Gamma \left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}\mathrm{d}tand\Gamma (b,z)={\int }_z^{\infty }{t}^{b-1}{e}^{-t}\mathrm{d}t.$$

(31)

Notice that in the above expression of $X_1(t;A)$, these gamma functions must be interpreted as random functions [22].

In Figure 1, Figure 2 and Figure 3, the 1-PDF of the exact solution SP, $f_1(x,t)$, and the 1-PDF of the truncation, $f_1^N(x,t)$, are plotted, respectively, for different values of the truncation order $N\in \{1,2,3,4,5\}$ at the time instants $t\in \{1.1,1.5,2\}$. In order to graphically check the convergence, in all of these figures, we plotted the approximation corresponding to the highest order ($N=5$) from Expression (21) together with the exact 1-PDF. In this latter case, the 1-PDF can been obtained by applying the RVT technique from the exact representation of the solution SP given by (28), (30), and (31) and using the software Mathematica${}^{\circledR }$. From this latter plot, we can observe in the three cases ($t\in \{1.1,1.5,2\}$) that $f_1(x,t)$ and $f_1^N(x,t)$, with $N=5$, practically coincide, so illustrating quick convergence. For the sake of completeness, in

Table 1. Error measure $e_N\left(t\right)$ defined by (32) for different time instants, $t\in \{1.1,1.5,2\}$, and orders of truncation of the approximate solution SP, $N\in \{1,2,3,4,5\}$. Example 1.

${\mathit{e}}_{\mathit{N}}\left(\mathit{t}\right)$	$\mathit{N}=1$	$\mathit{N}=2$	$\mathit{N}=3$	$\mathit{N}=4$	$\mathit{N}=5$
$t=1.1$	0.519832	0.166420	0.039233	0.007228	0.001129
$t=1.5$	0.437495	0.219520	0.071847	0.018712	0.004031
$t=2.0$	0.319464	0.249286	0.117028	0.0428193	0.0125606

, we show the values of the error measure $e_N\left(t\right)$, defined in (32), which has been calculated for each $t\in \{1.1,1.5,2\}$ using different orders of truncation $N\in \{1,2,3,4,5\}$:

$$e_N\left(t\right)={\int }_{-\infty }^{+\infty }|f_1^N(x,t)-f_1(x,t)|\mathrm{d}x.$$

(32)

From values collected in Table 1, we observe that for N fixed, the error increases as t goes far from the initial value $t=1$ (except for $N=1$ since the approximations are still rough), while for t fixed, the error decreases as the order of truncation N increases, as expected.

In Figure 4, we compare the expectation and the variance of the exact solution SP in the interval $t\in [1,2]$ with the corresponding approximations for different orders of truncation, $N\in \{1,2,3,4,5\}$. We measure the accuracy of the approximations for these two moments via the errors defined in (33),

$$e_N^{\mathbb{E}}={\int }_1^2\left|\mathbb{E}\left[X\left(t\right)\right]-\mathbb{E}\left[X^N\left(t\right)\right]\right|\mathrm{d}t,e_N^{\mathbb{V}}={\int }_1^2\left|\mathbb{V}\left[X\left(t\right)\right]-\mathbb{V}\left[X^N\left(t\right)\right]\right|\mathrm{d}t.$$

(33)

In

Table 2. Values of errors for the mean, $e_N^{\mathbb{E}}$, and the variance $e_N^{\mathbb{V}}$, given by (33) for $N\in \{1,2,3,4,5\}$. Example 1.

Error	$\mathit{N}=1$	$\mathit{N}=2$	$\mathit{N}=3$	$\mathit{N}=4$	$\mathit{N}=5$
$e_N^{\mathbb{E}}$	0.077488	0.037504	0.013075	0.003708	0.000903
$e_N^{\mathbb{V}}$	0.001776	0.001317	0.000343	0.000015	0.000012

, we collect the values for theses errors. For both moments, we observe that the corresponding error decreases as the order of truncation increases, as expected.

Example 2.

In this example, we consider the randomized Bessel equation:

$$t^2\ddot{X}\left(t\right)+t\dot{X}\left(t\right)+(t^2-A^2)X\left(t\right)=0,X\left(0.5\right)=Y_0,\dot{X}\left(0.5\right)=Y_1,$$

(34)

where $Y_0$ and $Y_1$ denote RVs fixed at the time instant $t_1=0.5$. In [12], the authors provided sufficient conditions on the input random data (A, $Y_0$, and $Y_1$) to extend the classical analysis for the Bessel equation to the random setting in the so-called mean squared sense [5]. In the aforementioned contribution, approximations for the mean and the variance of the solution SP to the Bessel RDE were given. In this example, we go further since, based on the results established in previous sections, we will construct approximations for the 1-PDF of the solution SP of this important RDE. To this end, let us first observe that, according to [12], the solution is given by $X\left(t\right)=Y_0{S}_1(t;A)+Y_1{S}_2(t;A)$, where $S_i(t;A)$, $i=1,2$ are defined by (17), being:

$${\displaystyle X_1(t;A)=t^A\left(1+\sum _{n=1}^{\infty }\frac{{(-1)}^n}{4^{n}n!{\prod }_{i=1}^n(A+i)}{t}^{2n}\right),X_2(t;A)=t^{-A}\left(1+\sum _{n=1}^{\infty }\frac{{(-1)}^n}{4^{n}n!{\prod }_{i=1}^n(-A+i)}{t}^{2n}\right).}$$

(35)

Hereinafter, we assume that A, $Y_0$, and $Y_1$ are independent RVs with the following distributions:

• A has a gamma distribution with shape and scale parameters $\alpha =2$ and $\beta =0.5$, respectively, truncated on the interval $I=[2.25,2.75]$, i.e., $A_I\sim \mathrm{Ga}(2;0.5)$.
• $Y_0$ has a beta distribution with parameters $\alpha =4$ and $\beta =3$, i.e., $Y_0\sim \mathrm{Be}(4;3).$
• $Y_1$ has a uniform distribution on the interval $[1,2]$, i.e., $Y_1\sim \mathrm{U}\left([1,2]\right)$.

Similarly to Example 1, we first check that Hypotheses H0–H3 are fulfilled:

• H0: Since A is a truncated RV on the interval $I=[2.25,2.75]$, A is bounded. With the notation of Hypothesis H0, notice that we can take, for instance, $M_A=2.75$.
• H1: Since A, $Y_0$, and $Y_1$ are assumed to be independent RVs, their joint PDF is given by the product of the PDF of each RV,

  $$f_{A,Y_0,Y_1}(a,y_0,y_1)=115.823a{e}^{-2a}\frac{6!}{3!2!}{y}_0^3{(1-y_0)}^2,\forall (a,y_0,y_1)\in [2.25,2.75]\times [0,1]\times [1,2].$$
  From this expression, it is clear that this function is continuous with respect to its second argument, $y_0$, and it is bounded. As we deal with a rectangular domain where the first factor is a decreasing function and the second factor, $y_0^3{(1-y_0)}^2$, has a maximum value at $y_0=0.6$ below $0.035$, the maximum of $f_{A,Y_0,Y_1}(a,y_0,y_1)$ is obtained taking this value and evaluating the first factor at $a=2.25$. Therefore, with the notation of H1, we can take, for instance, $M_f=6.1$. Therefore, Hypothesis H1 is fulfilled.
• H2: Notice that $t_0=0$ is a regular-singular point of the Bessel RDE (34). On the one hand, observe that according to Equation (2), $a_0\left(t\right)=t^2$, $a_1\left(t\right)=t$, and $a_2\left(t\right)=t^2-A^2$, being $a_0\left(t_0\right)=0$. On the other hand, according to (4) and (5), $p(t;A)=t/t^2$ and $q(t;A)=(t^2-A^2)/t^2$, which are not analytic functions at $t_0=0$, while $\widehat{p}(t;A)=1$ and $\widehat{q}(t;A)=t^2-A^2$ are both analytic at $t_0=0$ inasmuch as they are polynomials as functions of t, for every $a=a\left(\omega \right)$, $\omega \in \mathsf{\Omega }$.
• H3: Since A and $Y_1$ are assumed to be independent RVs, then the domain of the random vector $(A,Y_1)$ is:

  $${\mathcal{D}}_{A,Y_1}={\mathcal{D}}_A\times {\mathcal{D}}_{Y_1}=[2.25,2.75]\times [1,2].$$
  Hence, ${\mathcal{D}}_{A,Y_1}$ is a Lebesgue measurable set of ${\mathbb{R}}^2$, and with the notation of Hypothesis H3, we can take $\widehat{M}=2<\infty$. Therefore, Hypothesis H3 holds.

The 1-PDFs of the truncated solution SP, $f_1^N(x,t)$, for different orders of truncation, $N\in \{1,2,3,4,5\}$, at the time instants $t\in \{1.5,2.5,3.5\}$ are plotted in Figure 5, Figure 6 and Figure 7, respectively. As in Example 1, using the RVT method and the exact representation of series $X_1(t,A)$ and $X_2(t,A)$, defined in (35), in terms of the Bessel function $J_A\left(t\right)$ integrated in Mathematica${}^{\circledR }$ and the gamma function $\Gamma \left(z\right)$ defined in (31), namely $X_1(t,A)=2^A{J}_A\left(t\right)\Gamma (1+A)$ and $X_2(t,A)=2^{-A}{J}_{-A}\left(t\right)\Gamma (1-A)$, we can obtain the exact 1-PDF, $f_1(x,t)$. On the one hand, taking this benchmark, in the left panel of Figure 5, Figure 6 and Figure 7, we show the convergence of the approximations, $f_1^N(x,t)$, to the exact 1-PDF, $f_1(x,t)$, as N increases ($N\in \{1,2,3,4,5\}$) for different values of $t\in \{1.5,2.5,3.5\}$. On the other hand, to illustrate better the convergence, in the right panel of the aforementioned figures, we show the approximations corresponding to the highest order of truncation, $f_1^N(x,t)$, $N=5$, and the exact 1-PDF, $f_1(x,t)$. We observe that both plots overlap, thus showing quick convergence.

As in Example 1, we use the error $e_N\left(t\right)$ defined in (32) to measure the quality of approximations $f_1^N(x,t)$. In

Table 3. Values of the error measure, $e_N$, defined by (32) for different time instants, $t\in \{1.5,2.5,3.5\}$, and orders of truncation to the approximate solution SP of (34), $N\in \{1,2,3,4,5\}$. Example 2.

${\mathit{e}}_{\mathit{N}}$	$\mathit{N}=1$	$\mathit{N}=2$	$\mathit{N}=3$	$\mathit{N}=4$	$\mathit{N}=5$
$t=1.5$	0.03446	0.001593	0.000069	2.3272 $·{10}^{-6}$	5.5373 $·{10}^{-8}$
$t=2.5$	0.153784	0.030282	0.002078	0.000108	0.000020
$t=3.5$	1.518306	0.285614	0.038847	0.003435	0.000233

, we calculated, for each $t\in \{1.5,2.5,3.5\}$, the values of this error for the following orders of truncation $N\in \{1,2,3,4,5\}$. Notice that for each t, the error diminishes as the order of truncation increases, while the error increases as we move far from the initial data $t_1=0.5$, as expected.

Taking $f_1^N(x,t)$, given by (21), as an approximation to $f_1(x,t)$, in the expressions (1), we can calculate approximations for the mean ($\mathbb{E}\left[X^N\left(t\right)\right]$) and the variance ($\mathbb{V}\left[X^N\left(t\right)\right]$) of the solution of the random IVP (34). In Figure 8, we plotted the above-mentioned approximations on the interval $t\in [0.5,3.5]$ for different orders of truncation, $N\in \{1,2,3,4,5\}$, as well as their respective exact values, $\mathbb{E}\left[X\left(t\right)\right]$ and $\mathbb{V}\left[X\left(t\right)\right]$, using as the expression to $f_1(x,t)$ the one obtained via the RVT method. In these plots, we can clearly see that the approximations quickly improve on the whole interval as N increases. To finish, we calculated the errors of these approximations by means of the following expressions:

$$e_N^{\mathbb{E}}={\int }_{0.5}^{3.5}\left|\mathbb{E}\left[X^N\left(t\right)\right]-\mathbb{E}\left[X\left(t\right)\right]\right|\mathrm{d}t,e_N^{\mathbb{V}}={\int }_{0.5}^{3.5}\left|\mathbb{V}\left[X^N\left(t\right)\right]-\mathbb{V}\left[X\left(t\right)\right]\right|\mathrm{d}t,$$

(36)

Values of errors for the mean, $e_N^{\mathbb{E}}$, and the variance, $e_N^{\mathbb{V}}$, for each truncation are shown in

Table 4. Values of errors for the approximations of mean, $e_N^{\mathbb{E}}$, and the variance, $e_N^{\mathbb{V}}$, given by (36) for $N\in \{1,2,3,4,5\}$. Example 2.

Error	$\mathit{N}=1$	$\mathit{N}=2$	$\mathit{N}=3$	$\mathit{N}=4$	$\mathit{N}=5$
$e_N^{\mathbb{E}}$	6.766564	1.052105	0.108092	0.008096	0.000475
$e_N^{\mathbb{V}}$	20.769241	3.761294	0.352357	0.025039	0.001421

. From these values, we observe that errors decrease as N increases.

5. Conclusions

In this paper, we presented a methodology to construct reliable approximations to the first probability density function of the solution of second-order linear differential equations about a regular-singular point with a finite degree of randomness in the coefficients and also assuming that both initial conditions were random variables. Therefore, we studied this class of differential equations assuming randomness in all its input data, which provided great generality to our analysis. The results completed the ones already presented in a previous paper where the same class of random initial value problems was solved about an ordinary point. In this manner, the results obtained in both contributions permitted solving, from a probabilistic standpoint, a number of important randomized differential equations that appear in the field of mathematical physics like Airy, Hermite, Legendre, Laguerre, Chebyshev, Bessel, etc. Under our approach, we could compute the first probability density function of the solution. This is really advantageous since apart from allowing us to determine the mean and the variance, the knowledge of the first probability density function also permits determining higher one-dimensional moments, as well as the probability that the solution lies in any interval of interest, which might be key information in many practical problems.