Computing the two first probability density functions of the random Cauchy-Euler differential equation: Study about regular-singular points

Deterministic differential equations play a key role in many disciplines to model numerous phenomena. In practice, the application of models based on differential equations requires setting their inputs such as coefficients and initial conditions. These parameters are usually obtained by experiments where measurement errors are involved. In addition, there are external sources which can affect the physical system to be modelled. These facts motivate the treatment of inputs parameters as random variables (RV’s) or stochastic processes (SP’s) rather than deterministic constants or functions, respectively. In this contribution, we will focus on random differential equations (RDE’s), that is, differential equations whose coefficients are RV’s or SP’s.

As a main difference with respect to the deterministic framework, solving a RDE means not only to obtain its solution but also probabilistic information associated to the solution SP. Indeed, if W (u) denotes the solution SP of a RDE, then it is also important to compute its mean, μW(u)=E[W(u)]$\begin{array}{} \displaystyle \mu_W(u)=\mathbb{E}\left[W(u)\right] \end{array}$, and its variance, σW2(u)=V[W(u)]$\begin{array}{} \displaystyle \sigma_W^2(u)=\mathbb{V}\left[W(u)\right] \end{array}$. Additionally, the computation of the first probability density function (1-PDF), f^1(w,u)$\begin{array}{} \displaystyle \hat{f}_1(w,u) \end{array}$, is also important because from it one gets a full probabilistic description of the solution SP in each time instant u. Furthermore, from the 1-PDF one can compute all the one-dimensional statistical moments of W (u),

E[(W(u))k]=∫−∞∞wkf^1(w,u)dw,k=0,1,2,…$$\begin{array}{} \displaystyle \mathbb{E}\left[\left(W(u)\right)^k\right]=\int_{-\infty}^{\infty} w^k \hat{f}_1(w,u)\, \mathrm{d}w,\quad k=0,1,2,\dots \end{array}$$

Therefore, from it the mean and the variance can be straightforwardly obtained

E[W(u)]=∫−∞∞wf^1(w,u)dw,V[W(u)]=∫−∞∞w2f^1(w,u)dw−(E[W(u)])2.$$\begin{array}{} \displaystyle \mathbb{E}\left[W(u)\right]=\int_{-\infty}^{\infty} w\hat{f}_1(w,u) \,\mathrm{d}w,\quad \mathbb{V}\left[W(u)\right]=\int_{-\infty}^{\infty} w^2 \hat{f}_1(w,u)\, \mathrm{d} w-\left(\mathbb{E}\left[W(u)\right]\right)^2. \end{array}$$(1)

In general, to get more probabilistic information, the n-dimensional PDF’s of the solution SP might be computed, but it usually involves complex computations. For example, the 2-PDF, f^2(w1,u1;w2,u2)$\begin{array}{} \displaystyle \hat{f}_2(w_1,u_1; w_2, u_2) \end{array}$, provides a full probabilistic description of W (u) at every arbitrary pair of times, u₁ and u₂. In particular, from it the correlation function, Γ_W(u₁,u₂), can be computed. This function gives a measure of linear statistical interdependence between W (u₁) and W (u₂) and it is given by

ΓW(u1,u2)=E[W(u1)W(u2)]=∫−∞∞∫−∞∞w1w2f^2(w1,u1;w2,u2)dw1dw2.$$\begin{array}{} \displaystyle \Gamma_W(u_1,u_2)= \mathbb{E}\left[W(u_1)W(u_2)\right]= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} w_1w_2\hat{f}_2(w_1,u_1;w_2,u_2) \, \mathrm{d}w_1 \mathrm{d} w_2. \end{array}$$

Furthermore, Γ_W(u₁; u₂) allows us the computation of the covariance function

CW(u1,u2)=ΓW(u1,u2)−EW(u1)EW(u2).$$\begin{array}{}C_W(u_1,u_2)=\Gamma_W(u_1,u_2)-\mathbb{E}\left[W(u_1)\right]\mathbb{E}\left[W(u_2)\right]. \end{array}$$(2)

The Cauchy-Euler differential equation is adequate to model a number of phenomena in Engineering, particularly in Mechanics and in Theory of Potential Fields. For example, in this latter context it is applied to model the electric potential field between two concentric spheres [8, Sec. 2.5]. The formulation of this problem depends on the potential fields of both spheres, which usually are not known in a deterministic way but randomly due the heterogeneity of the surrounding medium. This motivates the study of the random Cauchy-Euler differential equation. Besides, the deterministic Cauchy-Euler differential equation has two regular-singular points, zero and infinity. Therefore, the goal of this paper is to obtain the 1-PDF and the 2-PDF of the solution SP of each one of two random IVP’s, on the one hand we shall consider a random IVP based on a Cauchy-Euler differential equation, and on the other hand, using an adequate transformation, infinity is moved to the origin. With this aim the results presented in [1] will be applied. In [1], the 1-PDF and 2-PDF of the solution SP of a random homogeneous linear second-order differential equation has been computed. To conduct our analysis, we will take advantage of the Random Variable Transformation (RVT) method that is stated below in Theorem 1. This technique allows us to obtain the PDF of a random vector obtained from the mapping of another random vector whose PDF is known.

Let us consider the following random IVP (4) based on a Cauchy-Euler differential equation

u2W″(u)+BuW′(u)+A2W(u)=0,u>u0>0,W(u0)=X0,W'(u0)=C,}$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{u^2}W''(u) + BuW'(u) + {A_2}W(u) = 0,\quad u > {u_0} > 0,}\\ {W\left( {{u_0}} \right) = {X_0},}\\ {W'\left( {{u_0}} \right) = C,} \end{array}\} \end{array}$$(4)

where the input parameters X₀, C, B and A₂ are assumed to be absolutely continuous RV’s, defined on a common probability space (Ω,ℱ,ℙ)$\begin{array}{} \displaystyle (\Omega,\mathfrak{F},\mathbb{P}) \end{array}$, with a joint PDF fX0,C,B,A2(x0,c,b,a2)$\begin{array}{} \displaystyle {f_{{X_0},C,B,{A_2}}}({x_0},c,b,{a_2}) \end{array}$. Thus, for the sake of generality we are implicitly assuming that involved RV’s are probabilistically dependent. As usual, hereinafter the notation X ≡ X(ω) will be used indistinctly. Our subsequent analysis is strongly related to the problem studied in [1], this motivates the notation used in IVP (4) for coefficients, initial conditions and unknown. In order to make clearer our notation, we will assume that the domain of every random inputs is an interval, although it is not necessary since we are assuming that X₀, C, B and A₂ are probabilistically dependent RV’s

DX0={x0=X0(ω),ω∈Ω:x0,1<x0<x0,2},DC={c=C(ω),ω∈Ω:c1<c<c2},DB={b=B(ω),ω∈Ω:b1<b<b2},DA2=​​​{a2=A2(ω),ω∈Ω:a2,1<a2<a2,2}.$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{{\cal D}_{{X_0}}}}& = &{\{ {x_0} = {X_0}(\omega ),{\mkern 1mu} \omega \in \Omega :{\mkern 1mu} {x_{0,1}} < {x_0} < {x_{0,2}}\} ,}\\ {{{\cal D}_C}}& = &{\{ c = C(\omega ),{\mkern 1mu} \omega \in \Omega :{\mkern 1mu} {c_1} < c < {c_2}\} ,}\\ {{{\cal D}_B}}& = &{\{ b = B(\omega ),{\mkern 1mu} \omega \in \Omega :{\mkern 1mu} {b_1} < b < {b_2}\} ,}\\ {{{\cal D}_{{A_2}}}}& = &{\{ {a_2} = {A_2}(\omega ),{\mkern 1mu} \omega \in \Omega :{\mkern 1mu} {a_{2,1}} < {a_2} < {a_{2,2}}\} .} \end{array} \end{array}$$

With the transformation u = e^tu₀, ∀u > u₀, IVP (4) is equivalent to IVP (5)

Z″(t)+A1Z′(t)+A2Z(t)=0,t>0,Z(0)=X0,Z'(0)=X1,}$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {Z''(t) + {A_1}Z'(t) + {A_2}Z(t) = 0,\quad t > 0,}\\ {Z\left( 0 \right) = {X_0},}\\ {Z'\left( 0 \right) = {X_1},} \end{array}\} \end{array}$$(5)

where A₁ = B − 1, X₁ = Cu₀ and Z(t) = W (e^tu₀). Then, the domains of RV’s A₁ and X₁ are defined as follows

DX1​​​​=​​{x1=X1(ω),ω∈Ω:c1u0<x1<c2u0},DA1​​​​=​​{a1=A1(ω),ω∈Ω:b1−1<a1<b2−1}.$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{{\cal D}_{{X_1}}} = \{ {x_1} = {X_1}(\omega ),{\mkern 1mu} \omega \in \Omega :{\mkern 1mu} {c_1}{u_0} < {x_1} < {c_2}{u_0}\} ,}\\ {{{\cal D}_{{A_1}}} = \{ {a_1} = {A_1}(\omega ),{\mkern 1mu} \omega \in \Omega {\mkern 1mu} :{\mkern 1mu} {b_1} - 1 < {a_1} < {b_2} - 1\} .} \end{array} \end{array}$$

This section is addressed to compute the 1-PDF and the 2-PDF of the solution SP to IVP (4), W (u). With this aim, first we will apply the theoretical results obtained in [1] to compute the 1-PDF and 2-PDF of the solution SP to IVP (5), Z(t). In [1], the 1-PDF and the 2-PDF of the random homogeneous linear second-order differential equation was determined and the mean, the variance and the covariance functions as well. Secondly, we will take advantage of Corollary 2, to get probabilistic information in terms of the joint PDF fX0,C,B,A2(x0,c,b,a2)$\begin{array}{} \displaystyle {f_{{X_0},C,B,{A_2}}}({x_0},c,b,{a_2}) \end{array}$, which is known. Finally, we will undo the change of variable u = e^tu₀ to obtain the 1-PDF and 2-PDF of the solution SP to the random IVP (4).

This section is addressed in the study of the infinity point of the randomized Cauchy-Euler differential equation. In the deterministic theory, in order to analyse u = +∞, the first step is to introduce the change of variable u = 1/s, and then to study a neighbourhood about the point s = 0 in the resulting expression. In our case, the transformed random IVP is the following

s2V″(s)+DsV′(s)+A2V(s)=0,0<s<s0,V(s0)=X0,V'(s0)=E,}$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{s^2}V''(s) + DsV'(s) + {A_2}V(s) = 0,\quad 0 < s < {s_0},}\\ {V\left( {{s_0}} \right) = {X_0},}\\ {V'\left( {{s_0}} \right) = E,} \end{array}\} \end{array}$$(17)

where D = 2 − B, E=−C/s02$\begin{array}{} \displaystyle E=-C/s_0^2 \end{array}$ and V (s) = W (1/s), ∀s : 0 < s < s₀.

The 1-PDF and 2-PDF of the solution SP, V (s), of the IVP (17) can be obtained following the same strategy exhibited in the previous section, because we are dealing with a random Cauchy-Euler differential equation too. We use the change of variable s = e⁻_ts₀, ∀s : 0 < s < s₀, and we obtain the following random IVP

Z″(t)+A1Z′(t)+A2Z(t)=0,t>0,Z(0)=X0,Z'(0)=X1,}$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {Z''(t) + {A_1}Z'(t) + {A_2}Z(t) = 0,\quad t > 0,}\\ {Z\left( 0 \right) = {X_0},}\\ {Z'\left( 0 \right) = {X_1},} \end{array}\} \end{array}$$(18)

where A₁ = (1 − D) = B − 1, X₁ = C/s₀, ∀t = log(s₀/s), ∀t > 0 and Z(t) = V (s₀ e^−t). Then, the domains of RV’s A₁ and X₁ are defined as follows

DA1​​​=​​​{a1=A1(ω),ω∈Ω:b1−1<a1<b2−1},DX1​​​=​​​{x1=X1(ω),ω∈Ω:c1s0<x1<c2s0}.$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{{\cal D}_{{A_1}}} = \{ {a_1} = {A_1}(\omega ),{\mkern 1mu} \omega \in \Omega {\mkern 1mu} :{\mkern 1mu} {b_1} - 1 < {a_1} < {b_2} - 1\} ,}\\ {{{\cal D}_{{X_1}}} = \{ {x_1} = {X_1}(\omega ),{\mkern 1mu} \omega \in \Omega {\mkern 1mu} :{\mkern 1mu} \frac{{{c_1}}}{{{s_0}}} < {x_1} < \frac{{{c_2}}}{{{s_0}}}\} .} \end{array} \end{array}$$

In this section we will show an example where the results obtained in Sections 2 and 3 are illustrated. As we have pointed out in (8), the 1-PDF and 2-PDF in both regular-singular problems depend on real and complex nature of the roots of the involved characteristic equation, having as associated probabilities p₁ and p₂, respectively. These probabilities are defined in (7). To account for the most interesting cases regarding our previous analysis, we shall consider the following three possibles scenarios

In the three cases we will consider a joint Gaussian distribution for the input parameters J_i = (X₀, C, B, A₂)^T ~ N(μ_i;Σ), i = 1,2,3 (i corresponding to Cases I, II and III, respectively), where

μi={​​(2,1,4,1)⊤ if i=1,​​(1,1,3,1)⊤ if i=2,​​(1,1,3,2)⊤ if i=3,Σ=110(4111141111211113).$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{\mu _i} = \{ \begin{array}{*{20}{c}} {{{(2,1,4,1)}^ \top }{\rm{ if }}i = 1,}\\ {{{(1,1,3,1)}^ \top }{\rm{ if }}i = 2,}\\ {{{(1,1,3,2)}^ \top }{\rm{ if }}i = 3,} \end{array}}&{\Sigma = \frac{1}{{10}}\left( {\begin{array}{*{20}{c}} 4&1&1&1\\ 1&4&1&1\\ 1&1&2&1\\ 1&1&1&3 \end{array}} \right)} \end{array}. \end{array}$$(24)

Probabilities p₁ and p₂ = 1 − p₁ are collected in Table 1. Notice that

Table 1

Columns p₁ and p₂ = 1 − p₁ collect the values of the probabilities given by (25) corresponding to Cases I–III, when J_i ~ N(μ_i,Σ), being μ_i and Σ specified in (24). Values of p_s represent the probabilities associated with asymptotic stability according to (26).

p1=ℙ[ω∈Ω:Δ(ω)>0]=ℙ[ω∈Ω:(B(ω)−1)2−4A2(ω)>0]=ℙ[ω∈Ω:(B(ω)−1)24>A2(ω)]=∫−∞∞[∫−∞(b−1)24fB,A2(b,c)da2]db,$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{p_1}}& = &{\left[ {\omega \in \Omega :\Delta (\omega ) > 0} \right] = \left[ {\omega \in \Omega :{{\left( {B(\omega ) - 1} \right)}^2} - 4{A_2}(\omega ) > 0} \right]}\\ {}& = &{\left[ {\omega \in \Omega :\frac{{{{\left( {B(\omega ) - 1} \right)}^2}}}{4} > {A_2}(\omega )} \right] = \int\limits_{ - \infty }^\infty {\left[ {\int\limits_{ - \infty }^{\frac{{{{(b - 1)}^2}}}{4}} {{f_{B,{A_2}}}} (b,c){\rm{d}}{a_2}} \right]} {\rm{d}}b,} \end{array} \end{array}$$(25)

where fB,A2(b,a2)=∫−∞∞∫−∞∞fX0,C,B,A2(x0,c,b,a2)dx0dc$\begin{array}{} \displaystyle {f_{B,{A_2}}}(b,{a_2}) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {{f_{{X_0},C,B,{A_2}}}} } ({x_0},c,b,{a_2}){\mkern 1mu} {\rm{d}}{x_0}{\mkern 1mu} {\rm{d}}c \end{array}$.

Based on the well-known condition that characterizes the asymptotic stability of the zero-steady state solution, Z(t) ≡ 0, to the deterministic counterpart of random IVP’s (5) and (18), that is A₁ > 0 and A₂ > 0, we will study the asymptotic stability of the solution W(u) ≡ 0 to the random IVP (4), hence also of the solution V (s) ≡ 0. This analysis relies on the computation of the following probability, p_s,

ps=ℙ[ω∈Ω:A1(ω)>0,A2(ω)>0]=ℙ[ω∈Ω:(B−1)(ω)>0,A2(ω)>0]=∫−∞∞∫−∞∞∫1∞∫0∞fX0,C,B,A2(x0,c,b,a2)da2dbdx0dc.$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {{p_s}}& = &{\left[ {\omega \in \Omega {\mkern 1mu} :{\mkern 1mu} {A_1}(\omega ) > 0,{\mkern 1mu} {A_2}(\omega ) > 0} \right] = \left[ {\omega \in \Omega {\mkern 1mu} :{\mkern 1mu} (B - 1)(\omega ) > 0,{\mkern 1mu} {A_2}(\omega ) > 0} \right]}\\ {}& = &{\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {\int\limits_1^\infty {\int\limits_0^\infty {{f_{{X_0},C,B,{A_2}}}} } } } ({x_0},c,b,{a_2}){\mkern 1mu} {\rm{d}}{a_2}{\mkern 1mu} {\rm{d}}b{\mkern 1mu} {\rm{d}}{x_0}{\mkern 1mu} {\rm{d}}c.} \end{array} \end{array}$$(26)

Values of p_s in the Cases I–III are shown in Table 1.

In Figure 1, the graphical representations for the 1-PDF, f^1(w,u)$\begin{array}{} \displaystyle \hat{f}_1(w,u) \end{array}$, at u ∊ {2,3,...,10} and f¯(v,s)$\begin{array}{} \displaystyle \bar{f}(v,s) \end{array}$, at v ε {0.05,0.1,...,0.5} in Cases I–III are shown. In these graphical representations, we can see how the PDFs evolves over the times u and s in each Case I–III. The behaviour of the PDFs is in full agreement with the plots shown in Figure 2 where the means μ_W(u) and μ_V(s) plus/minus the standard deviations σ_W(u) and σ_V(s) have been plotted for both problems. Regarding the IVP (4), in the three cases we observe as the mean function decreases as u increases and the standard deviation increases slowly. To the IVP (17) we can observe that the mean increase at first and then decrease. In addition, in Figure 3 the covariance function given by the expression (2) has been plotted to Case I for both problems.

Fig. 1

Fig. 2

Fig. 3

In this work we have given a full probabilistic solution to the randomized Cauchy-Euler differential equation under very general conditions. Specifically, we have obtained closed explicit formulas for the first and the second probability density functions of the solution stochastic process of the random Cauchy-Euler differential equation assuming that input parameters (coefficients and initial conditions) are absolutely continuous random variables with a joint probability density function. In this manner, all one-dimensional statistical moments of the solution can be computed including the mean, the variance and the covariance functions. The method allows us to construct punctual and probabilistic predictions in practical applications where the random Cauchy-Euler model appears. Furthermore, the technique applied in our study could be extended to perform an analogous analysis of other important randomized differential equations.