*Article* **Quadrature Integration Techniques for Random Hyperbolic PDE Problems**

**Rafael Company 1,\*, Vera N. Egorova <sup>2</sup> and Lucas Jódar <sup>1</sup>**


**Abstract:** In this paper, we consider random hyperbolic partial differential equation (PDE) problems following the mean square approach and Laplace transform technique. Randomness requires not only the computation of the approximating stochastic processes, but also its statistical moments. Hence, appropriate numerical methods should allow for the efficient computation of the expectation and variance. Here, we analyse different numerical methods around the inverse Laplace transform and its evaluation by using several integration techniques, including midpoint quadrature rule, Gauss–Laguerre quadrature and its extensions, and the Talbot algorithm. Simulations, numerical convergence, and computational process time with experiments are shown.

**Keywords:** random hyperbolic model; random laplace transform; numerical integration; monte carlo method; numerical simulation; talbot algorithm

### **1. Introduction**

Random hyperbolic partial differential equations (PDEs) are mathematical models that describe wave phenomena with applications in various fields: fluid mechanics [1,2], electromagnetic radiation [3], geosciences [4], and many others. The theory of hyperbolic problems has been well developed based on the assumption that parameters of the model, such as coefficients or initial values are exactly known, which is not available in the real world, where error measurement and the unavailability of the measurement occur. It causes the increasing interest for the random models, which can estimate the impact of the uncertainty to the predicted solution.

The solution is found numerically due to the complexity of random models. Following the mean square approach [5], we can extend existing numerical methods for deterministic problems to the random case by applying the Monte Carlo method [6,7] in order to approximate the statistical moments of the solution. Nevertheless, iterative numerical methods require the storage of the preliminary results and huge number of repetitions, which leads to the the necessity of enormous computational resources and makes them not appropriated to deal with random models. Thus, it becomes urgent to search for an accurate and fast numerical algorithm. Integral transform is a good alternative, as it allows us to construct the solution at one fixed point, not necessarily in the whole domain as it occurs in the case of the finite difference methods, as it is shown in the literature [8].

Integral transform methods convert the original random PDE to an ordinary differential equation (ODE), which can be solved analytically, in some cases, or numerically. Once obtained the solution of the random ODE, the inverse transform is applied in order o restore the solution of the original problem. This inverse transform can be done by the definition, i.e., integrating over the infinite domain, or by using some numerical techniques [9]. There are several widely used methods: Fourier Series, Stehfest approach [10], and Talbot inverse algorithm [11]. Because the inverse Laplace transform is ill-posed problem,

**Citation:** Company, R.; Egorova, V.N.; Jódar, L. Quadrature Integration Techniques for Random Hyperbolic PDE Problems. *Mathematics* **2021**, *9*, 160. https://doi.org/10.3390/ math9020160

Received: 7 October 2020 Accepted: 12 January 2021 Published: 14 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the regularization property of the numerical algorithm is necessary. In this sense, the Talbot inverse becomes the best option, since it guarantees the regularization property, while other numerical inversion schemes fail in dealing with noisy data [12].

In this work, we construct a numerical solution for random hyperbolic PDE models, not only by constructing the approximating stochastic process solution, but also while computing its expectation and variance. Thinking of practical applications, we deal with random models where the uncertainty is described by stochastic processes (s.p.'s) having a finite degree of randomness ([5], p. 37); this means that the involved s.p.'s take the form

$$\mathcal{G}(\mathbf{x}) = G(\mathbf{x}, V\_1, V\_2, \dots, V\_m), \tag{1}$$

where *Vi*, 1 ≤ *i* ≤ *m*, are mutually independent random variables (r.v.'s).

We propose an analytic-numerical approach that is based on random integral transform technique combined with various numerical integration methods, such as midpoint rule, Gauss-quadratures, and Talbot inverse [11]. The Monte Carlo method is used for the evaluations of the integrands involving the solution of random ordinary differential problems and also for the computation of the expectation and variance of the approximating stochastic process solution. The oscillatory nature of the appearing integrands deserves careful attention, because not all of the quadrature rules are advisable [13–15].

The proposed analytical–numerical approach for solving random hyperbolic PDE problems considered in this paper includes known state of the art of numerical integration methods, which are compared between themselves in terms of accuracy and computational time: the midpoint quadrature rule, the Talbot algorithm for Laplace inverse, the Gauss–Laguerre quadrature, the Exponential-Fitting Gauss-Laguerre quadrature, and the adaptative quadrature. This comparison is provided to highlight the advantages and drawbacks of each method. Moreover, this complex approach is compared with standard finite-difference methods for solving the random hyperbolic PDE problem. In all cases, the Monte Carlo simulations are used in order to calculate the statistical moments of the random solution process.

The rest of the paper is organized, as follows. In Section 2, the random hyperbolic PDE problem is formulated and the random Laplace transform method is briefly described. Section 3 proposes numerical integration methods for Laplace inverse, while Section 4 gives an algorithm for Monte Carlo simulations. All of the proposed methods are compared by the series of numerical tests in Section 5. Section 6 discusses the results.

All of the numerical tests have been carried out by MatLAB, version R2020a, for Windows 10 Home (64-bit), Intel(R) Core(TM) i5-8265U CPU, 1.60 GHz.

#### **2. Preliminaries and Integral Transform for Random Hyperbolic PDE**

This section begins by recalling previous results and definitions [8,16]. Let us consider a complete probability space (Ω, <sup>F</sup>, <sup>R</sup>) and the set *Lp* with the *<sup>p</sup>*-norm of a real-values random variable *<sup>Y</sup>* <sup>∈</sup> *Lp*(Ω), as defined by

$$||Y||\_p = \left(\mathbb{E}[|Y|^p]\right)^{1/p}, \quad p \ge 1,\tag{2}$$

where the expectation <sup>E</sup>[|*Y*<sup>|</sup> *p*] <sup>&</sup>lt; <sup>∞</sup>, and *Lp*(Ω) is a Banach space [17]. By using definition (2), the integrability, continuity, and differentiability of a function *<sup>Y</sup>*(*t*) <sup>∈</sup> *Lp*(Ω) can be defined straightforwardly.

Note that, if *<sup>p</sup>* = 2, then it is a mean square (m.s.) case. Let <sup>C</sup> be the class of all m.s. locally integrated two-stochastic processes (s.p.'s) *<sup>h</sup>*(*t*) defined in <sup>R</sup> such that *<sup>h</sup>*(*t*) = 0 , for all negative arguments and the two-norm satisfies

$$\exists \mathcal{c} \ge 0, \ M > 0: \quad \|h(t)\|\_2 \le M \exp(ct), \quad \forall t \ge 0. \tag{3}$$

Subsequently, for *<sup>h</sup>*(*t*) ∈ C, the m.s. integral

$$H(s) = \mathcal{L}[h(t)](s) = \int\_0^\infty h(t) \exp(-st) dt,\tag{4}$$

where *<sup>s</sup>* is a complex number with real part *Re*(*s*) <sup>&</sup>gt; *<sup>c</sup>*<sup>0</sup> <sup>≥</sup> 0, and it is called the random Laplace transform of 2-s.p. *<sup>h</sup>*(*t*). The constant *<sup>c</sup>*<sup>0</sup> is chosen, such that *Re*(*s*) <sup>&</sup>gt; *<sup>c</sup>*<sup>0</sup> specifies the region where *<sup>H</sup>*(*s*) is analytic and it has some form of singularity on the line *Re*(*s*) = *<sup>c</sup>*<sup>0</sup> [9]. If *<sup>H</sup>*(*s*) is known, then the random inverse transform for *<sup>t</sup>* <sup>&</sup>gt; 0 is defined, as follows

$$h(t) = \frac{1}{2\pi i} \int\_{a - i\infty}^{a + i\infty} H(s) \exp(st) ds,\tag{5}$$

where *i* stands for the imaginary unit and *α* > *c*<sup>0</sup> [16].

For the purposes of present study, we recall some of the important properties of the random Laplace transform (4): if s.p. *<sup>h</sup>*(*t*) is twice m.s. differentiable and *<sup>h</sup>* (*t*), *<sup>h</sup>* (*t*) belong to C, then

$$
\mathcal{L}[h'(t)](s) = sH(s) - h(0+), \qquad \mathcal{L}[h''(t)](s) = s^2H(s) - sh(0+) - h'(0+). \tag{6}
$$

In this paper, we consider a one-dimensional random hyperbolic PDE modelling the s.p. of the vibrating string motion *<sup>u</sup>*(*x*, *<sup>t</sup>*), depending on the spatial variable *<sup>x</sup>* and time *<sup>t</sup>*,

$$u\_{tt}(\mathbf{x},t)(\underline{\xi}) = a(\mathbf{x})(\underline{\xi})u\_{xx}(\mathbf{x},t)(\underline{\xi}) + b(\mathbf{x})(\underline{\xi})u\_x(\mathbf{x},t)(\underline{\xi}) + c(\underline{\xi})u(\mathbf{x},t)(\underline{\xi}), \quad \mathbf{x} \in [0,L], \ t > 0, \ \underline{\xi} \in \Omega,\tag{7}$$

$$u(\mathbf{x},0)(\xi) = f\_0(\mathbf{x})(\xi), \quad u\_t(\mathbf{x},0)(\xi) = f\_1(\mathbf{x})(\xi), \tag{8}$$

$$f\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha\_{\alpha}}\right)}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}\dots}}}}}}}$$
}}\}}}}}}}}}}}}}}}}}}

$$
\mu(0,t)(\xi) = \mathfrak{g}\_0(t)(\xi), \qquad \mathfrak{u}(L,t)(\xi) = \mathfrak{g}\_1(t)(\xi), \tag{9}
$$

where *<sup>a</sup>*(*x*)(*ξ*) <sup>&</sup>gt; 0, *<sup>b</sup>*(*x*)(*ξ*) are m.s.-continuous stochastic processes with a finite degree of randomness and absolutely integrable with respect to the spatial variable in <sup>R</sup>; *<sup>c</sup>*(*ξ*) is a random variable (r.v.). The s.p.'s *<sup>f</sup>*0(*x*)(*ξ*), *<sup>f</sup>*1(*x*)(*ξ*), *<sup>g</sup>*0(*t*)(*ξ*), and *<sup>g</sup>*1(*t*)(*ξ*) are functions depending on a finite number of r.v. that represent random initial and boundary conditions with a finite degree of randomness.

The random hyperbolic partial differential equation (PDE) (7) is solved using an analytic-numerical method that is based on Laplace transform combined with an appropriate numerical integration technique. In this paper, we consider various quadratures for inverse Laplace transform.

Following the ideas of [8,18], let us define the random Laplace transform with respect to the temporal variable, as

$$
\mathcal{U}(\mathbf{x}, \mathbf{s})(\xi) = \mathcal{L}[u(\mathbf{x}, t)(\xi)].\tag{10}
$$

Because *<sup>u</sup>*(*x*, *<sup>t</sup>*)(*ξ*) is a twice m.s. differentiable s.p., one gets

$$\mathcal{L}[u\_{\rm tr}(\mathbf{x},t)(\xi)] = \mathbf{s}^2 \mathcal{U}(\mathbf{x},\mathbf{s})(\xi) - \mathbf{s}u(\mathbf{x},0)(\xi) - u\_{\rm tr}(\mathbf{x},0)(\xi) = \mathbf{s}^2 \mathcal{U}(\mathbf{x},\mathbf{s})(\xi) - sf\_0(\mathbf{x})(\xi) - f\_1(\mathbf{x})(\xi). \tag{11}$$

Subsequently, (7) is transformed to the following random non-homogeneous ordinary differential equation (ODE) with respect to the spatial variable

$$a(\mathbf{x})(\xi)\mathcal{U}\_{\mathbf{x}\mathbf{x}}(\mathbf{x},\mathbf{s})(\xi) + b(\mathbf{x})(\xi)\mathcal{U}\_{\mathbf{x}}(\mathbf{x},\mathbf{s})(\xi) + (\mathbf{c} - \mathbf{s}^2)\mathcal{U}(\mathbf{x},\mathbf{s})(\xi) = -[sf\_0(\mathbf{x})(\xi) + f\_1(\mathbf{x})(\xi)],\tag{12}$$

for *<sup>x</sup>* <sup>∈</sup> [0, *<sup>L</sup>*], *<sup>ξ</sup>* <sup>∈</sup> <sup>Ω</sup>.

Assuming *<sup>a</sup>*(*x*)(*ξ*) <sup>&</sup>gt; 0 for each event *<sup>ξ</sup>* <sup>∈</sup> <sup>Ω</sup>, one gets

$$\mathcal{U}\_{\mathbf{X}\mathbf{X}}(\mathbf{x},\mathbf{s})(\boldsymbol{\xi}) + \frac{b(\mathbf{x})(\boldsymbol{\xi})}{a(\mathbf{x})(\boldsymbol{\xi})} \mathcal{U}\_{\mathbf{x}}(\mathbf{x},\mathbf{s})(\boldsymbol{\xi}) + \frac{\mathbf{c} - \mathbf{s}^{2}}{a(\mathbf{x})(\boldsymbol{\xi})} \mathcal{U}(\mathbf{x},\mathbf{s})(\boldsymbol{\xi}) = -\frac{sf\_{0}(\mathbf{x})(\boldsymbol{\xi}) + f\_{1}(\mathbf{x})(\boldsymbol{\xi})}{a(\mathbf{x})(\boldsymbol{\xi})}.\tag{13}$$

Equation (13) is a linear second order ODE with respect to the spatial variable, which can be analytically solved in some cases, or numerically in other cases. Because the boundary conditions (9) for the PDE are functions on *t*, the boundary conditions for (13) are the corresponding Laplace transforms of (9):

$$\mathcal{U}(0,\mathbf{s})(\xi) = \mathcal{L}[\mathcal{g}\_0(t)(\xi)], \quad \mathcal{U}(L,\mathbf{s})(\xi) = \mathcal{L}[\mathcal{g}\_1(t)(\xi)].\tag{14}$$

Once obtaining the solution *<sup>U</sup>*(*x*,*s*)(*ξ*), a real-valued *<sup>u</sup>*(*x*, *<sup>t</sup>*)(*ξ*) is restored by while using random inverse Laplace transform that is given by (5). Taking advantage of the relationship between the inverse Laplace transform and Fourier cosine integrals, see [9], the following formula is used

$$u(\mathbf{x},t)(\xi) = \frac{2e^{\alpha t}}{\pi} \int\_0^\infty \text{Re}[\ell I(\mathbf{x}, \mathbf{a} + i\nu)(\xi)] \cos(wt) dw, \quad \xi \in \Omega,\tag{15}$$

where *Re*[·] stands for the real part of a complex number. Note that the integrand appearing in (15) has an oscillatory kernel that deserves special care for the numerical integration.

#### **3. Numerical Integration Methods**

This section describes briefly acknowledged integration methods for the integrals of the type (15).

THe numerical solution of Equation (7) is constructed in the domain <sup>Δ</sup> = [0; *<sup>L</sup>*] <sup>×</sup> [0; *<sup>T</sup>*] for each fixed event *ξ*. Let us introduce a uniform grid {*xj*, *t <sup>n</sup>*}, such that

$$\mathbf{x}\_{j} = j \mathbf{h}\_{\prime} \; h = \frac{L}{N\_{x}} \; j = 0, \dots, N\_{x}; \quad \mathbf{t}^{n} = n \mathbf{k}\_{\prime} \; k = \frac{T}{N\_{\mathbf{t}}}, \; n = 0, \dots, N\_{\mathbf{t}}.\tag{16}$$

At each node (*xj*, *<sup>t</sup> n*), the numerical solution is defined by *<sup>u</sup><sup>n</sup> <sup>j</sup>* (*ξ*) for each realization of *<sup>ξ</sup>* and it is obtained by approximating the integral (15). Hence, at every fixed (*xj*, *<sup>t</sup> n*), the following function is defined

$$f\_{\vec{j},n}(w) = f\_{\vec{j},n}(w,\vec{\xi}) = \text{Re}\left[lI(\mathbf{x}\_{\vec{j}}, a+iw)(\vec{\xi})\right] \cos(wt^n),\tag{17}$$

where *<sup>U</sup>*(*xj*, *<sup>α</sup>* <sup>+</sup> *iw*)(*ξ*) is the numerical solution of ODE (13) at the point *xj* for fixed value of *<sup>s</sup>* = *<sup>α</sup>* + *iw*. Now, we briefly describe all of the considered methods for numerical integration.

#### *3.1. Midpoint Quadrature Rule*

The midpoint quadrature rule is a method of approximation of integral (15) based on the Riemann sums, the simplest case of Newton–Cotes open formulas, for truncated domain [0, *<sup>R</sup>*]. In the general case, the midpoint quadrature rule is written, as follows

$$\int\_0^\infty f(w)dw \approx \int\_0^R f(w)dw = \sum\_{k=0}^N f(w\_{k+1/2})h\_{MP} + O(h^2),\tag{18}$$

where *wk*+1/2 <sup>=</sup> *<sup>k</sup>* <sup>+</sup> <sup>1</sup> 2 *hMP*, *hMP* <sup>=</sup> *<sup>R</sup> <sup>N</sup>* , *<sup>k</sup>* <sup>=</sup> 0, . . . , *<sup>N</sup>* <sup>−</sup> 1.

It is well known that the main advantage of this method is its simplicity of implementation and the consideration of all the information regarding the integrand, which makes it applicable for a wider class of integrand functions [14]. However, the high accuracy of the quadrature requires large enough value of *N*, leading to the increasing computational cost. In the case of improper integral (in the infinite domain), the method can also be sensitive to the choice of *R*.

#### *3.2. Gauss-Laguerre Quadrature*

The novelty of Gauss quadratures is to choose nodes where the integrand is evaluated in order to minimize the error of approximation. It is a good alternative to Newton–Cotes formulas, especially when the evaluation of function itself requires a lot of computational resources, because good accuracy can be reached with a small number, four or five, of nodes if the integrand is well conditioned. This is not the case when the integrand is of oscillatory type [19].

The improper integral is approximated by Gauss–Laguerre (GL) quadrature of *NGL* nodes by the following sum, see [20],

$$\int\_0^\infty f(w)dw \approx \sum\_{k=1}^{N\_{GL}} \gamma\_k f(w\_k) \epsilon^{w\_k} \tag{19}$$

where *wk* is the *<sup>k</sup>*-th root of Laguerre polynomial *LNGL* (*w*), *<sup>γ</sup><sup>k</sup>* is the weight of the quadrature given by

$$\gamma\_k = \frac{w\_k}{(N\_{\rm GL} + 1)^2 \left[L\_{\rm NGL} + 1 \, (w\_k)\right]^2}, \quad k = 1, \ldots, N\_{\rm GL}. \tag{20}$$

#### *3.3. Exponentially-Fitted Gauss-Laguerre Quadrature*

Exponential fitting is an approach that is used in numerical differentiation, interpolation, and integration for improving the accuracy of the methods. Because integrand in (15) is oscillating, Exponentially-fitted Gauss–Laguerre quadrature (EF-GL), as proposed in [21], could be a good option. For EF-GL, nodes and weights depend on integrand and cannot be defined a priori. The computation of these *NGL* pairs of nodes and weights is based on the solution of a nonlinear system of *NGL* equations, which leads to additional computational cost. In [21], the numerical algorithm is described in details. Further, in Section 5, we compare the accuracy and computational time of GL and EF-GL quadrature rules.

#### *3.4. Talbot Inverse*

The method of Talbot for the Laplace inversion problem [11] is based on numerical contour integration. Instead of formula (15), the Bromwich integral is used

$$
\mu\_j^n(\xi) = \frac{1}{2\pi i} \int\_{a-i\infty}^{a+i\infty} e^{s\xi^n} \mathcal{U}(x\_{j'}s) ds. \tag{21}
$$

The contour deformation is used in order to obtain the Hankel contour and exploit the exponential factor, which makes the integral suitable for further application of a Newton– Cotes formula [22]. The Talbot inversion quadrature for *NT I* nodes is written, as follows

$$\mu\_{\dot{\jmath}}^{n}(\xi) = \frac{2}{5t^{n}} \sum\_{k=0}^{N\_{\text{II}}-1} \text{Re}\left[\gamma\_{k} \text{II}(\mathbf{x}\_{\dot{\jmath}\prime} \frac{\mathbf{w}\_{k}}{t^{n}})\right],\tag{22}$$

where *wk* are the nodes and *γ<sup>k</sup>* are the weights defined by

$$w\_0 = \frac{2N\_{TI}}{5}, \quad w\_k = \frac{2\pi k}{5} \left( \cot \left( \frac{k\pi}{N\_{TI}} + i \right) \right), \tag{23}$$

$$\gamma\_0 = 0.5 \exp(w\_0), \quad \gamma\_k = \exp(w\_k) \cdot \left[ 1 + \frac{k\pi}{N\_{TI}} \left( 1 + \cot^2\left(\frac{k\pi}{N\_{TI}}\right) - i \cot\left(\frac{k\pi}{N\_{TI}}\right) \right) \right]. \tag{24}$$

Here, the number of nodes *NT I* should be chosen in accordance with desired accuracy: for *<sup>n</sup>* significant digits *NT I* <sup>=</sup> 1.7*n*. It shows the flexibility of the method and the high degree of accuracy with fast convergence. Moreover, as stated in [12], the main advantage of the Talbot algorithm is the regularization property, which means the ability to handle noisy data. It is important for the inverse Laplace transform problem due to its ill-posedness and

it becomes even more urgent in the random case dealing with perturbed initial conditions or parameters of the problem.

Summarizing, a numerical solution is constructed following the steps of Algorithm 1 for all of the described methods.

#### **Algorithm 1:** Numerical solution for deterministic string vibrating problem

Initialization: set the mesh {*xj*, *t <sup>n</sup>*} by (16); Set initial conditions *<sup>u</sup>*(*xj*, 0) = *<sup>f</sup>*0(*xj*), *<sup>j</sup>* = 0, . . . , *Nx*; Set *α* > *c*0; Set number of nodes of the quadrature *N*; Set *<sup>n</sup>* = 0; **while** *t <sup>n</sup>* < *T* **do** Increment *n*; **for** *<sup>j</sup>* <sup>=</sup> 0, . . . , *Nx* **do** Compute nodes and weigths {*wk*, *γk*} of the chosen quadrature - Midpoint rule: uniform grid with *N* nodes ; - GL quadrature: nodes *wk* are the roots of the Laguerre polynomial of *<sup>N</sup>*-th order, *<sup>k</sup>* <sup>=</sup> 1, . . . , *<sup>N</sup>*; - EF-GL quadrature [21]: nodes *wk* and weights *γ<sup>k</sup>* are found by solving nonlinear system of 2*N* equations; - Talbot inverse: nodes *wk* and weights *<sup>γ</sup>k*, *<sup>k</sup>* <sup>=</sup> 0, . . . , *<sup>N</sup>* <sup>−</sup> 1 are defined by (23)–(24); Get the approximated value *u<sup>n</sup> j* : - Midpoint rule: integral in (15) is approximated by (18); - GL and EF-GL quadratures: integral in (15) is approximated by (19)–(20); - Talbot inverse: formula (22) ; **end end**

#### **4. Monte Carlo Method for Random Hyperbolic PDE**

The coefficients of the random m.s. Equation (7) and corresponding initial and boundary conditions (9) are stochastic processes (s.p's) that are defined in a complete probability space (Ω, <sup>F</sup>, <sup>P</sup>), i.e., s.p.'s *<sup>a</sup>*(*x*), *<sup>b</sup>*(*x*), *<sup>f</sup>*0(*x*), *<sup>f</sup>*1(*x*), *<sup>g</sup>*0(*x*) and *<sup>g</sup>*1(*x*) are described as continuous s.p.'s with with one-degree of randomness.

The solution of the random m.s. problem is approximated by using the the Monte Carlo approach [6,7], when the expectation <sup>E</sup>[*u*(*x*, *<sup>t</sup>*)] is approximated by the average of a sufficiently large number of realizations *ξ* ∈ Ω of the corresponding deterministic realized transformed random ordinary differential problem. The Algorithm 2 describes the steps of the numerical solution.

#### **Algorithm 2:** Numerical solution for random hyperbolic PDE problem

Initialization: set the mesh {*xj*, *t <sup>n</sup>*} by (16); Set number of the MC realizations *NMC*; Generate *NMC* random variables for s.p.'s *<sup>a</sup>*(*x*), *<sup>b</sup>*(*x*), *<sup>f</sup>*0(*x*), *<sup>f</sup>*1(*x*), *<sup>g</sup>*0(*x*); Choose the method of numerical integration **for** *<sup>m</sup>* <sup>=</sup> 1, . . . *NMC* **do** Define s.p.'s *<sup>a</sup>*(*x*), *<sup>b</sup>*(*x*), *<sup>f</sup>*0(*x*), *<sup>f</sup>*1(*x*), *<sup>g</sup>*0(*x*) for fixed realization; Run Algorithm 1 to obtain the numerical solution *um* of the deterministic problem; Increment *m*; **end** Compute <sup>E</sup>[*u*] = <sup>∑</sup>*NMC <sup>m</sup>*=<sup>1</sup> *um NMC* ; Compute <sup>E</sup>[*u*2] = <sup>∑</sup>*NMC u*2 *m NMC* ;

*<sup>m</sup>*=<sup>1</sup> Compute <sup>0</sup>Var[*u*] = <sup>0</sup>E[*u*2] <sup>−</sup> (E[*u*])<sup>2</sup>

#### **5. Numerical Results**

This section deals with the comparison of the above-described methods of numerical integration and Laplace inversion for several test problems.

#### *5.1. Deterministic PDE Problem with Constant Coefficients*

We start with simple one dimensional deterministic problem with a known analytical solution in order to check the viability of the proposed numerical integration techniques. The deterministic example corresponds to one fixed event *ξ* ∈ Ω. Instead of the bounded spatial domain [0; *<sup>L</sup>*], the whole real axis <sup>R</sup> is considered. Thus, no boundary conditions are needed. We also assume that *a* > 0, *b*, and *c* are constants, i.e., the following wave equation is considered

$$u\_{tt}(\mathbf{x},t) = a^2 u\_{xx}(\mathbf{x},t) + b u\_x(\mathbf{x},t) + c u(\mathbf{x},t), \quad \mathbf{x} \in \mathbb{R}, t > 0,\tag{25}$$

subject to initial conditions *<sup>u</sup>*(*x*, 0) = *<sup>f</sup>*0(*x*), *ut*(*x*, 0) = *<sup>f</sup>*1(*x*).

This problem admits an analytical solution that can be written in terms of Bessel function of the first kind, see [23], p. 574, Equation 6.1.5, as follows

• for *<sup>c</sup>* <sup>−</sup> <sup>1</sup> <sup>4</sup> *<sup>a</sup>*−2*b*<sup>2</sup> <sup>=</sup> *<sup>σ</sup>*<sup>2</sup> <sup>&</sup>gt; 0:

$$\begin{split} u(\mathbf{x},t) &= \frac{1}{2}f(\mathbf{x}+at)\exp\left(\frac{bt}{2a}\right) + \frac{1}{2}f(\mathbf{x}-at)\exp\left(-\frac{bt}{2a}\right) \\ &+ \frac{\sigma t}{2a}\exp\left(-\frac{bx}{2a^2}\right)\int\_{\mathbf{x}-at}^{\mathbf{x}+at} \exp\left(\frac{b^x\_{\xi}}{2a^2}\right) \frac{I\_1\left(\sigma\sqrt{t^2-(\mathbf{x}-\xi)^2/a^2}\right)}{\sqrt{t^2-(\mathbf{x}-\xi)^2/a^2}}f(\xi)d\xi \\ &+ \frac{1}{2a}\exp\left(-\frac{bx}{2a^2}\right)\int\_{\mathbf{x}-at}^{\mathbf{x}+at} \exp\left(\frac{b^x\_{\xi}}{2a^2}\right)I\_0\left(\sigma\sqrt{t^2-(\mathbf{x}-\xi)^2/a^2}\right)g(\xi)d\xi,\end{split} \tag{26}$$

where *<sup>I</sup>*0(*z*) and *<sup>I</sup>*1(*z*) are the modified Bessel function of the first kind; • for *<sup>c</sup>* <sup>−</sup> <sup>1</sup> <sup>4</sup> *<sup>a</sup>*−2*b*<sup>2</sup> <sup>=</sup> <sup>−</sup>*σ*<sup>2</sup> <sup>&</sup>lt; 0:

$$\begin{split} u(\mathbf{x},t) &= \frac{1}{2}f(\mathbf{x}+at)\exp\left(\frac{bt}{2a}\right) + \frac{1}{2}f(\mathbf{x}-at)\exp\left(-\frac{bt}{2a}\right) \\ &- \frac{\sigma t}{2a}\exp\left(-\frac{bx}{2a^2}\right)\int\_{\mathbf{x}-at}^{\mathbf{x}+at} \exp\left(\frac{b\mathbf{y}^{\mathbf{x}}}{2a^2}\right) \frac{l\_1\left(\sigma\sqrt{t^2-(\mathbf{x}-\boldsymbol{\xi})^2/a^2}\right)}{\sqrt{t^2-(\mathbf{x}-\boldsymbol{\xi})^2/a^2}}f(\boldsymbol{\xi})d\boldsymbol{\xi} \\ &+ \frac{1}{2a}\exp\left(-\frac{bx}{2a^2}\right)\int\_{\mathbf{x}-at}^{\mathbf{x}+at} \exp\left(\frac{b\mathbf{y}^{\mathbf{x}}}{2a^2}\right)l\_0\left(\sigma\sqrt{t^2-(\mathbf{x}-\boldsymbol{\xi})^2/a^2}\right)g(\boldsymbol{\xi})d\boldsymbol{\xi},\end{split} \tag{27}$$

where *<sup>J</sup>*0(*z*) and *<sup>J</sup>*1(*z*) are Bessel function of the first kind.

In order to test the proposed numerical integration methods we apply Laplace transform, as described in Section 2, and obtain a deterministic version of Equation (13):

$$d\mathcal{U}\_{xx}(\mathbf{x},\mathbf{s}) + \frac{b}{a^2} l \mathcal{U}\_{\mathbf{x}}(\mathbf{x},\mathbf{s}) + \frac{c - s^2}{a^2} l \mathcal{U}(\mathbf{x},\mathbf{s})(\boldsymbol{\zeta}) = -\frac{s f\_0(\mathbf{x}) + f\_1(\mathbf{x})}{a^2}.\tag{28}$$

Applying the non-unitary Fourier transform with angular frequency

$$
\hat{\mathcal{U}}(w,s) = \mathcal{F}[\mathcal{U}(\mathbf{x},s)] = \int\_{-\infty}^{\infty} \mathcal{U}(\mathbf{x},s) \exp(-i\mathbf{x}w)d\mathbf{x},\tag{29}
$$

Equation (13) takes the following form

$$-w^2\hat{\mathcal{U}}(w,s) + iw\frac{b}{a^2}\hat{\mathcal{U}}(w,s) + \frac{c-s^2}{a^2}\hat{\mathcal{U}}(w,s) = -\mathcal{F}\left[\frac{sf\_0(\mathbf{x}) + f\_1(\mathbf{x})}{a^2}\right].\tag{30}$$

Algebraic Equation (30) is solved directly

$$\hat{\mathcal{U}}(w,s) = \frac{-\mathcal{F}\left[\frac{sf\_0(x) + f\_1(x)}{a^2}\right]}{-w^2 + iw\frac{b}{a^2} + \frac{c - s^2}{a^2}}.\tag{31}$$

Hence, the solution *<sup>U</sup>*(*x*,*s*) of (28) can be obtained by applying inverse Fourier transform to (31).

In the next Example 1, we consider a particular case of Equation (25) with constant coefficients and trigonometric initial conditions.

**Example 1.** *Let us consider deterministic problem* (25) *with coefficients <sup>a</sup>* = <sup>2</sup>*, <sup>b</sup>* = <sup>1</sup>*, <sup>c</sup>* = <sup>3</sup>*, and initial conditions f*0(*x*) = cos(*x*) *and f*1(*x*) = sin(*x*)*.*

Numerical solution is constructed in the truncated domain [0, *<sup>L</sup>*] <sup>×</sup> [0, *<sup>T</sup>*], *<sup>L</sup>* = 5, *<sup>T</sup>* = 1, for discrete uniformly distributed nodes (16), *Nx* <sup>=</sup> 5, *Nt* <sup>=</sup> 10, by applying the described in previous section methods, see Algorithm 1. We set *<sup>α</sup>* = 1.

Applying the inverse Fourier transform to (31), one obtains

$$\mathcal{U}(\mathbf{x},s) = \frac{1}{2\pi} \left[ \frac{\pi e^{-\mathbf{x}i}(\mathbf{s}+i)}{a^2 + bi + (\mathbf{s}^2 - \mathbf{c})} + \frac{\pi e^{\mathbf{x}i}(\mathbf{s}-i)}{a^2 - bi + (\mathbf{s}^2 - \mathbf{c})} \right], \quad \mathbf{s} = \mathbf{a} + i\mathbf{w}, \tag{32}$$

where *i* is the imaginary unit. Once the solution of ODE (28) is obtained, formula (15) is used to restore the solution of the PDE while using various numerical integration techniques.

Note that Equation (25) admits the analytical solution, as described above. Because the function *<sup>u</sup>*(*x*, *<sup>t</sup>*) is close to zero, we compute the relative error of the discrete numerical solution at the mesh nodes in order to estimate the accuracy of the methods

$$\text{RelError}(j, n) = \frac{|u\_{\text{ref}}(\mathbf{x}\_{j'}, t^n) - \mathcal{U}\_{\text{num}}(\mathbf{x}\_{j'}, t^n)|}{|u\_{\text{ref}}(\mathbf{x}\_{j'}, t^n)|},\tag{33}$$

where *<sup>U</sup>*num is the matrix of numerical solution *<sup>U</sup>*num <sup>=</sup> {*u<sup>n</sup> <sup>j</sup>* }, *<sup>j</sup>* <sup>=</sup> 0, ... , *Nx*, *<sup>n</sup>* <sup>=</sup> 0, ... , *Nt*, as computed by Algorithm 1; *<sup>u</sup>*ref(*xj*, *<sup>t</sup> n*) is the reference value at the point (*xj*, *<sup>t</sup> n*). In this example, as the exact solution is known, the reference value is equal to this exact solution. For other cases where the exact solution is not available, a reference value is obtained using accurate finite difference method (FDM) for solving the original PDE (7). The total computational time for the proposed methods are presented in Table 1, together with the maximum of RelErr(*j*, *<sup>n</sup>*).

The adaptative quadrature (MatLAB function integral [24]) has the same order of accuracy as the midpoint rule, but it requires greater computational resources. Thus, it will not be considered in further more complicated examples.

For the Talbot algorithm *<sup>M</sup>* = 17 is chosen to guarantee the accuracy up to 10 significant digits [22]. Even in that case, the method performs much faster than standard numerical integration methods for (15). Thus, the Talbot inverse method is found to be the most effective method for the deterministic case with constant coefficients.

**Table 1.** Comparison of various numerical methods for problem (25) with *<sup>a</sup>* = 2, *<sup>b</sup>* = 1 and *<sup>c</sup>* = <sup>3</sup> (Example 1).


The relative errors for Midpoint rule and Talbot inverse methods are plotted in Figure 1. Because no boundary conditions are posed for the problem, the largest values of the relative errors are situated at the boundary *<sup>x</sup>* = *<sup>L</sup>*.

**Figure 1.** Distribution of the relative error among space and time for midpoint rule (**left**) and Talbot inverse (**right**) methods in Example 1.

Table 2 presents a comparison of GL and EF-GL quadratures in terms of the maximum relative error and the CPU-time varying the number of nodes *NGL*. It is important to notice that the CPU-time may vary from simulation to simulation, thus only the order should be taken into account. In the case of GL quadrature, we find out that the computational time is similar with increasing number of nodes, while the CPU-time for EF-GL method is increasing exponentially. The convergence of the GL quadrature is shown, while taking the results shown in Table 1 into account: the error reduces significantly with an increasing number of nodes. The potential improvement of the GL method by the exponential fitting expectedly has higher computational cost, due to the solution of the nonlinear system at each point of the computational domain. However, the accuracy of the EF-Gl quadrature for this example with oscillating integrand has not been improved when comparing with the standard GL rule. Thus, it will not be considered in further more complicated examples.


**Table 2.** Gauss–Laguerre (GL) and Exponentially-fitted Gauss–Laguerre quadrature (EF-GL) methods results, depending on number of nodes of the quadrature for Example 1.

The accuracy of the midpoint rule depends on the truncation *R* and step size *hMP*. A bigger domain, as well as smaller step size, lead to an increased computational time. Figure <sup>2</sup> presents the plots of errors and the CPU-time for fixed step size *hMP* <sup>=</sup> <sup>10</sup>−<sup>1</sup> with respect to increasing domain. The accuracy in dependence on the step size *hMP* is also studied. In Table 3, the maximum relative error is reported for various *hMP* and fixed *<sup>R</sup>* = 104. The maximum relative error is decreasing with step size until 4.4699 <sup>×</sup> <sup>10</sup>−<sup>7</sup> (*hMP* <sup>=</sup> 1/16); further fragmentation of the step size does not reduce the error for *<sup>R</sup>* <sup>=</sup> <sup>10</sup>4.


**Table 3.** The maximum relative error and computational time of the midpoint quadrature rule with respect to step size *hMP* for Example 1.

**Figure 2.** Error and total computational time of the midpoint rule with respect to the domain size *<sup>R</sup>* with fixed *<sup>h</sup>* = <sup>10</sup>−<sup>1</sup> for Example 1.

#### *5.2. Deterministic PDE with Non-Constant Coefficients*

In the case of non-constant coefficients in (13), the analytical solution is not always available; thus, FDM is applied to construct a reference numerical solution. Note that the function *<sup>U</sup>*(*x*,*s*) that is used in expression (17) means the value of the numerical solution of the ODE (13) at the fixed point *x* for fixed parameter *s*.

Equation (13) is discretized by the central differences on the same mesh {*xj*}, *<sup>j</sup>* = 0, . . . , *Nx*, as follows

$$\frac{l\mathcal{U}\_{j+1} - 2\mathcal{U}\_j + \mathcal{U}\_{j-1}}{h^2} + \frac{b(\mathbf{x}\_j)}{a(\mathbf{x}\_j)} \frac{\mathcal{U}\_{j+1} - \mathcal{U}\_{j-1}}{2h} + \frac{c - s^2}{a(\mathbf{x}\_j)} \mathcal{U}\_j = -\frac{sf\_0(\mathbf{x}\_j) + f\_1(\mathbf{x}\_j)}{a(\mathbf{x}\_j)}, \quad j = 1, \dots, N\_\mathbf{x} - 1,\tag{34}$$

where *Uj* stands for the approximated value of *<sup>U</sup>*(*x*,*s*) at the node *xj*. The values at the boundaries are found from the boundary conditions by applying the Laplace transform

$$\mathcal{U}l\_0 = \mathcal{L}[\mathcal{g}\_0(t)], \quad \mathcal{U}\_{\mathcal{N}\_x} = \mathcal{L}[\mathcal{g}\_1(t)]. \tag{35}$$

Hence, the integrand (17) has to be evaluated at each fixed node of the computational grid in order to approximate integral (15), which provokes a significant augment of the CPU-time. In the next example, we increase the complexity by regarding a variable coefficients deterministic problem.

Because the analytical solution for the deterministic PDE problem in general form (7) is not available, a numerical method has to be employed to obtain the reference numerical

solution. We consider an explicitly centred in time and space finite difference scheme for the mesh function *u<sup>n</sup> <sup>j</sup>* <sup>≈</sup> *<sup>u</sup>*(*xj*, *<sup>t</sup> n*):

$$\frac{u\_j^{n+1} - 2u\_j^n + u\_j^{n-1}}{(\Delta t)^2} = a(x\_j) \frac{u\_{j+1}^n - 2u\_j^n + u\_{j-1}^n}{(\Delta x)^2} + b(x\_j) \frac{u\_{j+1}^n - u\_{j-1}^n}{2\Delta x} + c u\_j^n,\tag{36}$$

where *<sup>j</sup>* = 1, ... , *Nx*, *<sup>n</sup>* = 2, ... , *Nt*. The initial conditions (8) are used in order to obtain the solution at the first time levels *t* <sup>0</sup> and *t* 1. The derivative in (8) is approximated by the forward difference. Because the considered scheme is conditionally stable, the step sizes Δ*t* and Δ*x* are chosen to guarantee the stability. In order to obtain a good approximation, which could be considered as the reference solution, the mesh should be chosen appropriately fine.

**Example 2.** *Let us consider a deterministic vibrating string problem* (7) *on rectangle* [0, *<sup>L</sup>*] <sup>×</sup> [0, *<sup>T</sup>*]*, <sup>L</sup>* = 0.5*, <sup>T</sup>* = 0.2*. We set non-constant coefficients <sup>a</sup>*(*x*) = <sup>9</sup>*<sup>x</sup>* + <sup>1</sup>*, <sup>b</sup>* = <sup>−</sup>*ex, <sup>c</sup>* = <sup>−</sup>5*, initial conditions <sup>f</sup>*0(*x*) = *<sup>x</sup>*(*<sup>x</sup>* <sup>−</sup> *<sup>L</sup>*) *and <sup>f</sup>*1(*x*) = <sup>0</sup>*, and boundary conditions <sup>g</sup>*0(*t*) = *<sup>g</sup>*1(*t*) = <sup>0</sup>*.*

The numerical solution is constructed by the Algorithm 1, choosing *Nx* <sup>=</sup> 10, *Nt* <sup>=</sup> 5. For the midpoint rule, *<sup>N</sup>* = 100 and *<sup>R</sup>* = 100 are used. Table <sup>4</sup> presents the comparison of the methods in terms of maximum relative error and computational time. The reference solution is the numerical solution that is computed by the FDM (36) in refined mesh (*Nx* <sup>=</sup> 100, *Nt* = 16,000), which preserves the stability of the scheme. Because an explicit method is used and no iterative procedures are needed for solving nonlinear system at each time-level, the total computational time is comparably small: 0.15 s. Figure 3 plots the reference solution.

**Table 4.** A comparison of various methods of numerical integration for Example 2.


Figure <sup>4</sup> plots the solution at the moment *<sup>t</sup>* = *<sup>T</sup>*. The midpoint rule and Talbot inverse method perform more accurately than GL quadrature of nine nodes, but they require more computational time due to larger number of calls of integrand (17). However, taking 25 nodes in the GL quadrature, the accuracy has been improved significantly.

**Figure 3.** Reference solution for Example 2 computed by the finite difference method (FDM) (36).

**Figure 4.** Numerical solution for Example <sup>2</sup> at the moment *<sup>t</sup>* = *<sup>T</sup>* obtained by considered methods.

#### *5.3. Random PDE with Constant Coefficients*

In this subsection, we deal with random models with constant coefficient random variables. It is remarkable that, in this case, we need not only the computation of the approximation s.p. solution, but also the computation of its statistical moments.

**Example 3.** *We consider a random version of problem* (25)*, with <sup>a</sup>* ∼ N (2, 0.25)*, <sup>b</sup>*, *<sup>c</sup>* <sup>∼</sup> *Beta*(2, 5)*. In order to approximate the mean and variance of the solutions, the Monte Carlo method with NMC simulations is used.*

Expectation and variance of the exact solution for the random hyperbolic PDE (25) are plotted in Figure 5. As in previous examples, we compare the proposed methods of integration and Laplace inverse in terms of maximum relative error and computational time. Table 5 presents the results for various *NMC*. The CPU-time refers to the total computational

time for all *NMC* simulations. Note that, for 1000 simulations, the exact solution (26)–(27) requires 28.41 s to perform the simulations. Thus, Midpoint rule (*<sup>R</sup>* = 100, *<sup>h</sup>* = 0.1), Talbot inverse and GL quadrature require less computational time than calculation by the exact formula. As expected, the computational time is increasing with the number of simulations linearly, but errors preserve the order in most cases.

**Figure 5.** Expectation and variance of the exact solution for the random hyperbolic partial differential equation (PDE) (25) with *<sup>a</sup>* ∼ N (2, 0.25), *<sup>b</sup>*, *<sup>c</sup>* <sup>∼</sup> Beta(2, 5), performed using the Monte-Carlo method with *NMC* <sup>=</sup> <sup>103</sup> simulations.


**Table 5.** Comparison of various methods of numerical integration for the random hyperbolic PDE (25) with *<sup>a</sup>* ∼ N (2, 0.25), *<sup>b</sup>*, *<sup>c</sup>* <sup>∼</sup> Beta(2, 5).

*5.4. Random PDE with Non-Constant Coefficients*

To complete the study, a random variable coefficient problem is considered.

**Example 4.** *The vibration of the string in* [0, *<sup>L</sup>*] *is described by Equation* (7)*, subject to the initial conditions <sup>f</sup>*0(*x*) = *<sup>x</sup>*(*<sup>x</sup>* <sup>−</sup> *<sup>L</sup>*) *and <sup>f</sup>*1(*x*) = <sup>0</sup>*; and boundary conditions <sup>g</sup>*0(*t*) = *<sup>g</sup>*1(*t*) = <sup>0</sup>*. We set up the parameters:*

*<sup>L</sup>* <sup>=</sup> 0.5, *<sup>T</sup>* <sup>=</sup> 0.2, *<sup>a</sup>*(*x*) = *<sup>ϕ</sup><sup>x</sup>* <sup>+</sup> 1, *<sup>ϕ</sup>* ∼ N (9, 0.5), *<sup>b</sup>*(*x*) = <sup>−</sup>*ex*, *<sup>c</sup>* <sup>∼</sup> *Beta*(2, 5). (37)

Unlike the deterministic Example 2 with non-constant coefficients where FDM provides a reference analytical solution, reference values are not available here due to the computational complexity that arises in the evaluation of the statistical moments of the approximate stochastic process when time step advances [18]. A survival reference FDM solution is taking the Monte Carlo method for an appropriate set of realizations. In this case, the number of realizations is *NMC* <sup>=</sup> <sup>10</sup><sup>3</sup> and CPU-time is 16,212 s.

Figure 6 plots the numerical solution. The zero-variance at the boundaries is caused by the boundary conditions. Similar plots are obtained for the considered methods. Thus, we compare them in terms of the maximum relative error, see Table 6. As it is expected from the previous examples, the most accurate solution is obtained by the midpoint rule and Talbot inverse, although this advantage pays the price of additional computational cost.

**Figure 6.** Expectation and variance of the numerical solution for the random hyperbolic PDE (25) with *ϕ* ∼ <sup>N</sup> (9, 0.5), *<sup>b</sup>*(*x*) = <sup>−</sup>*ex*, *<sup>c</sup>* <sup>∼</sup> Beta(2, 5), performed by using the Monte-Carlo method with *NMC* <sup>=</sup> <sup>103</sup> simulations.

**Table 6.** Comparison of various methods of numerical integration for Example 4, *NMC* <sup>=</sup> 1000.


#### **6. Conclusions**

The solution of a random hyperbolic PDE problem is a challenging task that is demanded in many practical applications. Computing an expression of the approximating stochastic process makes the computation of its statistical moments available. In this paper, we propose a combination of the random Laplace transform with the numerical integration

techniques for its inverse, and the Monte Carlo method for the evaluation of numerical solution of the transformed problem at a particular required point.

The Monte Carlo simulations require a fast and efficient basis numerical algorithm for solving deterministic hyperbolic PDE problem, for every fixed realization. FDM could not be an option due to the high computational cost and memory requirements. In order to avoid the numerical differentiation of the PDE, Laplace transform is applied, which results in ODE equation. In some cases, as it has been shown in present paper, the analytical solution of ODE is known; thus, we use numerical integration methods for inverse Laplace transform. If the solution of ODE is not available, then numerical techniques for boundary value problem have to be employed.

Several numerical integration methods have been considered: midpoint rule and GL-quadrature for improper integrals. However, due to the oscillatory behaviour of the integrand function GL quadrature with a small number of nodes shows comparatively poor results, while the midpoint rule is comparable with Talbot's Laplace inverse for random hyperbolic PDEs. The proposed complex analytic-numerical approach is compared with the classical explicit FDM scheme for the original random PDE problem.

**Author Contributions:** R.C., V.N.E. and L.J. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been funded by the Spanish Ministerio de Economía, Industria y Competitividad (MINECO), the Agencia Estatal de Investigación (AEI) and Fondo Europeo de Desarrollo Regional (FEDER UE) grant MTM2017-89664-P.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** We state that data are available to the readers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

