## Abstract

Fractional calculus provides a rigorous mathematical framework to describe anomalous stochastic processes by generalizing the notion of classical differential equations to their fractional-order counterparts. By introducing the fractional orders as uncertain variables, we develop an operator-based uncertainty quantification framework in the context of stochastic fractional partial differential equations (SFPDEs), subject to additive random noise. We characterize different sources of uncertainty and then, propagate their associated randomness to the system response by employing a probabilistic collocation method (PCM). We develop a fast, stable, and convergent Petrov–Galerkin spectral method in the physical domain in order to formulate the forward solver in simulating each realization of random variables in the sampling procedure.

## 1 Introduction

Fractional models construct a tractable mathematical framework to describe and predict the behavior of multiscales multiphysics complex phenomena. Particularly, fractional differential equations, as a well-structured generalization of their integer-order counterparts, provide a rigorous mathematical tool to develop models that describe anomalous behavior in complex physical systems [19], where the anomaly manifests itself in heavy tail, sharp peaks, intermittency, and asymmetry in the distribution of corresponding underlying stochastic processes. Significant approximations as inherent part of assumptions upon which the model is built, lack of information about true values of parameters (incomplete data), and random nature of quantities being modeled pervade uncertainty in the corresponding mathematical formulations [10,11]. In this work, we develop an uncertainty quantification (UQ) framework in the context of stochastic fractional partial differential equations (SFPDEs), in which we characterize different sources of uncertainties and further propagate the associated randomness to the system response quantity of interest. The intention of this work is not to introduce new mathematical theories or methods for UQ, but rather to bring forward practical solutions using existing theories in an attempt to overcome the computational challenges of UQ in fractional models.

### 1.1 Types and Sources of Uncertainty.

The model uncertainties are in general being classified as aleatory and epistemic according to their fundamental essence. It is important to retain the separation between these two sources in order to assess the predictive efficiency of the model [12,13]. Aleatory uncertainty impacts the output of interest due to the natural variation of inputs and parameters; it is irreducible and commonly treated with probability theory. Epistemic uncertainty, however, results from a lack of knowledge about the system of interest and can be reduced by obtaining additional information. The epistemic uncertainties are broadly characterized as (i) model uncertainties, occurring in model inputs, numerical approximation errors, and model form uncertainty; and (ii) data uncertainties due to measurement inaccuracy and sparse or imprecise data. The model uncertainty encompasses all model parameters coming from geometry, constitutive laws, and fields equation, while also pertaining to surrounding interactions, such as boundary conditions and random forcing sources (noise). Numerical approximations, which are an essence of differential equations since they generally do not lend themselves to analytical solutions, introduce uncertainty by imposing different sources of discretization error, iterative convergence error, and round off error. In this work, we only consider the epistemic uncertainty in our fractional model and thus, introduce the fractional derivative orders as a new set of model parameters in addition to model coefficients. We note that the values of these new parameters are strongly tied to the distribution of underlying stochastic process and their statistics are estimated from experimental observations in practice, see, for example, Refs. [14] and [15].

### 1.2 Uncertainty Framework.

Conventional approaches in parametric UQ of differential equations are based around Monte Carlo sampling (MCS) [16], which performs ensemble of forward calculations to map the uncertain input space to the uncertain output space. This method enjoys from being embarrassingly parallelizable and can be implemented quite readily on high-dimensional random spaces. However, the key issue is the slow rate of convergence $∼1/N$ with N number of realization, which consequently imposes exhaustively so many operations of forward solver and makes the method not practical for expensive simulations. Other methods such as sequential MCS [17] and multilevel sequential MCS [18] are also developed and recently used in Ref. [19] to improve the parametric uncertainty assessment in elliptic nonlocal equations. An alternative to expensive MCS is to build surrogate models. An extensive comparison of two widely used ones, namely, polynomial chaos and Gaussian process, are provided in a recent work [20]. Polynomial chaos, in which the output of the stochastic model is represented as a series expansion of input parameters was initially applied in Ref. [21] and later extended and used in Refs. [2226]. It is also generalized and used in constructing stochastic Galerkin methods [2730] for problems with higher-dimensional random spaces. Other nonsampling numerical methods, including but not limited to the perturbation method [3134] and moment equation method [35,36] are also developed; however, their applications are restricted to stochastic systems with relatively low-dimensional random space. These so-called “intrusive” approaches typically do not treat the forward solver as a black box, rather require some knowledge and reformulation of the governing equations, and thus, may not be practical in many problems with complex codes.

A wide range of nonintrusive techniques mostly stretch oversampling, quadrature, and regression, see Ref. [20] and references therein. More recently, high-order probabilistic collocation methods (PCM), employing the idea of interpolation/collocation in the random spaces, are developed in Refs. [3739]; they are also known as Stochastic Collocation. These methods limit the sample points to an efficient subset of random space, while adequately sampling the necessary range. The excellence in the use of PCM is twofold; it has the benefit of easily sampling at nodal points that naturally leads to independent realizations of the problem as in MCS, and the advantage of fast convergence rate. The challenging postprocessing of solution statistics, which can essentially be described as a high-dimensional integration problem, can also be resolved by adopting sparse grid generators, such as the Smolyak algorithm [39,40]. The use of sparse grids, as opposed to full tensor product construction from one-dimensional quadrature rules, will effectively reduce the number of sampling, while preserving a fast convergence rate to a high level of accuracy.

### 1.3 Forward Solver.

A core task in computational forward UQ is to form an efficient numerical method, which for each realization of random variables can accurately solve and simulate the deterministic counterpart of the stochastic model in the physical domain. Such a numerical method is usually called forward solver or simulator. In the case of fractional partial differential equations (FPDEs), the excessive cost of numerical approximations becomes more challenging as FPDEs usually do not lend themselves to analytical solutions and more importantly, most of uncertainty propagation techniques instruct operations of forward solver many times. This requires the implementation of more efficient numerical schemes. In general, there are two main issues of nonlocality and end-points singularities in numerically solving FPDEs. The nonlocal feature of fractional derivatives causes the local methods such as finite difference method and finite element method to lose their advantage, see Refs. [4148] for some examples on finite difference method for solving FPDEs. In contrast, global and high-order schemes like spectral methods are proper techniques since their main disadvantage in treating standard differential equations becomes an advantage in the case of FPDEs. We refer to Refs. [49] and [50] for introduction to several spectral methods for integer-order differential equations and to Refs. [43] and [5158] for FPDEs. The solution of FPDEs also exhibits end-points singularities as another difficulty in developing high-accuracy numerical methods; an example is the inadequacy of a spectral method that uses only polynomial basis functions. More recently, Zayernouri et al. [59,60] developed two new spectral theories on fractional and tempered fractional Sturm–Liouville problems and introduced explicit corresponding eigenfunctions, namely, Jacobi polyfractonomials of first and second kind. These eigenfunctions are comprised of smooth and fractional parts, where the latter can be tuned to capture singularities of true solution. They are successfully employed in constructing discrete solution/test function spaces and developing a series of high-order and efficient Petrov–Galerkin spectral methods, see Refs. [6172].

The main focus of this work is to develop an operator-based computational forward UQ framework in the context of a stochastic fractional partial differential equation. Assuming that the mathematical model under consideration is well-posed and accounts in principle for all features of underlying phenomena, we identify three main sources of uncertainty, (i) parametric uncertainty, including fractional indices as new set of random parameters appeared in the operator, (ii) additive noises, which incorporates all intrinsic/extrinsic unknown forcing sources as lumped random inputs, and (iii) numerical approximations. Computational challenges arise when the admissible space of random inputs is infinite-dimensional, for example, problems subject to additive noise [73], and thus, the framework involves uncertainty parametrization via a finite number of random space basis. Unlike the classical approach in modeling random inputs, which considers idealized uncorrelated processes (white noises), we model the random inputs as more/fully correlated random processes (colored noises) and parametrize them via Karhunen–Loève (KL) expansion by assuming finite-dimensional noise assumption. This yields the problem in finite dimensional random space. We then, propagate the parametric uncertainties into the system response by applying PCM. We obtain the corresponding deterministic FPDE for each realization of random variables, using the Smolyak sparse grid generators for low to moderately high dimensions. In order to formulate the forward solver, we follow Ref. [65] and develop a high-order Petrov–Galerkin (PG) spectral method to solve for each realization of SFPDE in the physical domain by employing Jacobi polyfractonomials in addition to Legendre polynomials as temporal and spatial basis/test functions, respectively. The smart choice of coefficients in the construction of spatial basis/test functions yields symmetric properties in the resulting mass/stiffness matrix, which is then exploited to formulate an efficient fast solver. We also show that for each realization of random variables, the deterministic problem is mathematically well posed and the proposed forward solver is stable. By adopting a sufficient number of basis in the physical domain, we eliminate the epistemic uncertainty associated with numerical approximation and isolate the impact of parametric uncertainty on system response quantity of interest.

The organization of the article is as follows. We recall some preliminaries on fractional calculus in Sec. 2. Then, we formulate the stochastic system in Sec. 3 and parametrize the random inputs. We also develop the stochastic sampling, namely, PCM and MCS for our stochastic problem. We further construct the deterministic solver in Sec. 4 and provide the numerical results in Sec. 5. We end the article with a conclusion and summary.

## 2 Preliminaries on Fractional Calculus

Let $ξ∈[−1,1]$. The left- and right-sided fractional derivatives of order σ are defined as (see, e.g., Refs. [74] and [75])
$(−1RLDξσ)u(ξ)=1Γ(n−σ)dndξn∫−1ξu(s)ds(ξ−s)σ+1−n, ξ>−1$
(1)

$(ξRLD1σ)u(ξ)=1Γ(n−σ)(−d)ndξn∫ξ1u(s)ds(s−ξ)σ+1−n, ξ<1$
(2)
respectively. An alternative approach in defining the fractional derivatives is the left- and right-sided Caputo derivatives of order σ, $n−1<σ≤n, n∈ℕ$, defined as
$(−1CDξσu)(ξ)=1Γ(n−σ)∫−1ξu(n)(s)ds(ξ−s)σ+1−n, ξ>−1$
(3)

$(ξCD1σu)(ξ)=1Γ(n−σ)∫ξ1u(n)(s)ds(s−ξ)σ+1−n, ξ<1$
(4)
By performing an affine mapping from the standard domain $[−1,1]$ to the interval $t∈[a,b]$, we obtain
$aRLDtσu=(2b−a)σ(−1RLDξσ u)(ξ)$
(5)

$aCDtσu=(2b−a)σ(−1CDξσ u)(ξ)$
(6)
Hence, we can perform the operations in the standard domain only once for any given σ and efficiently utilize them on any arbitrary interval without resorting to repeating the calculations. Moreover, the corresponding relationship between the Riemann–Liouville and Caputo fractional derivatives in $[a,b]$ for any $σ∈(0,1)$ is given by
$(aRLDtσ u)(t)=u(a)Γ(1−σ)(t−a)σ+(aCDtσ u)(t)$
(7)

## 3 Forward Uncertainty Framework

### 3.1 Formulation of Stochastic Fractional Partial Differential Equation.

Let $D=[0,T]×[a1,b1]×[a2,b2]×⋯×[ad,bd]$ be the physical computational domain for some positive integer d and stochastic function $u(t,x;ω):D×Ω→ℝ$, where $ω∈Ω$ denotes the random input of the system in a properly defined complete probability space $(Ω,F,ℙ)$. We consider the following SFPDE, subject to certain homogeneous Dirichlet initial/boundary conditions and stochastic process as additional force function, given as
$Lq(ω) u(t,x;ω)=F(t,x;ω)$
(8)

$u |t=0=0$
(9)

$u |x=aj=u |x=bj=0$
(10)
such that for $ℙ$-almost everywhere $ω∈Ω$ the equation holds. The stochastic fractional operator and force term are given, respectively, as
$Lq(ω) =0Dtα(ω)−∑j=1d kj [ajDxjβj(ω)+xjDbjβj(ω)]$
(11)

$F(t,x;ω)=h(t,x)+f(t;ω)$
(12)

where the fractional indices $α(ω)∈(0,1)$ and $βj(ω)∈(1,2), j=1,2,…d$ are mutually independent random variables, kj are real positive constant coefficients, and the fractional derivatives are taken in the Riemann–Liouville sense. We assume that the driving terms h and f are properly posed, such that Eqs. (8)(10) is well-posed $ℙ$-almost everywhere $ω∈Ω$, and also the solution in the physical domain $D$ is smooth enough such that we can construct a numerical scheme to solve each realization of SFPDE. As an extension to future works, the stochastic operator Eq. (11) can be extended to $α(ω)∈(1,2)$ for the case of wave equations, and thus applied in formulating fractional models to study complex time-varying nonlinear fluid–solid interaction phenomena [7678], material damage [79], and also the effect of damping/stiffness in structural vibrations [8083].

### 3.2 Representation of the Noise: Dimension Reduction.

We approximate the additional random forcing term by representing $f(t;ω)$ into its finite-dimensional version and thus, reduce the infinite-dimensional probability space to a finite-dimensional space. This is achieved via truncating KL expansion with the desired accuracy [84].

Let $(Ω,F,ℙ)$ be a complete probability space, where Ω is the space of events, $F⊂2Ω$ denotes the σ-algebra of sets in Ω, and $ℙ$ is the probability measure. The random field $f(t;ω)$ has the ensemble mean $E{f(t;ω)}=f¯(t)$, finite variance $E{[f(t;ω)−f¯(t)]2},$ and covariance $Cf(t1,t2)=E{[f(t1;ω)−f¯(t1)][f(t2;ω)−f¯(t2)]}$. The KL expansion of $f(t;ω)$ takes the form
$f(t;ω)=f¯(t)+∑k=1∞λk ψk(t) Qk(ω)$
(13)
where $Q(ω)={Qk(ω)}|k=1k=∞$ is a set of mutually uncorrelated random variables with zero mean and unit variance, while $ψk(t)$ and λk are the eigenfunction and eigenvalues of the covariance kernel $Cf(t1,t2)$. We obtain the covariance kernel Cf and its eigenvalues and eigenfunctions, following Ref. [85] and by solving a stochastic Helmholtz equation
$△f(t;ω)−m2f(t;ω)=g(t;ω)$
(14)
where the random forcing $g(t;ω)$ is a white-noise process with zero mean and unit variance. The eigenvalues and eigenvectors of Eq. (14) form a Fourier series, so that the KL expansion Eq. (13) is replaced with its sine Fourier series version
$f(t;ω)=f¯(t;ω)+∑k=1∞ ak sin(2kπ tT) Qk(ω)$
(15)
in which the random variables $Qk(ω)$ are chosen to be uniformly distributed with probability density function $ρk(qk)$. T is the length of the process along the t-axis, and the coefficients
$ak=2Tℓ2[1+(2πkTℓ)2]−1$
(16)
where $ℓ=T/A$ and A is the correlation length of $f(t;ω)$. To ensure that the random variables $Qk(ω)$ have zero mean and unit variance, we define them on $Qk(ω)∈[−3,3]$. We note that this process is consistent with the zero-Dirichlet initial condition given in Eq. (9). Next, in order to render Eq. (15) computable, we truncate the infinite series with a prescribed ($≈90%$) fraction of the energy of the process, following the finite-dimensional noise assumption in stochastic computations. To this end, we set T =1, the correlation length $A=T/2$, and consider only the first four terms in the KL expansion. Let $fM(t;ω)=1μ ∑k=1Mak sin(2kπ tT) Qk(ω)$ denote the normalized truncated expansion, assuming $f¯M(t;ω)=0$, where $μ=maxt{σfM}$ and $σfM$ is the standard deviation of $fM(t;ω)$. Thus, we represent the random process to be employed in Eq. (8) as
$f(t;ω)=ϵfM(t;ω)$
(17)

where ϵ is the amplitude of process.

Therefore, the formulation of SFPDE Eq. (8) can be posed as follows: Find $u(t,x;ω):D×Ω→ℝ$ such that $∀t,x∈D$
$0Dtα(ω)u(t,x;ω)−∑j=1d kj [ajDxjβj(ω)+xjDbjβj(ω)]u(t,x;ω)=h(t,x)+f(t;Q1(ω),Q2(ω),…,QM(ω))$
(18)

holds $ℙ$-almost surely for $ω∈Ω$, subject to the homogeneous initial and boundary conditions.

### 3.3 Input Parametrization.

Let $Z:Ω→ℝN$ be the set of $N=1+d+M$ independent random parameters, given as
$Z={ Zi }i=1N={ α(ω),β1(ω),β2(ω),…,βd(ω),Q1(ω),Q2(ω),…,QM(ω) }$
with probability density functions $ρi:Γi→ℝ, i=1,2,…,N$, where their images $Γi≡Zi(Ω)$ are bounded intervals in $ℝ$. The joint probability density function
$ρ(ξ)=∏i=1Nρi(Zi), ∀ξ∈Γ$
(19)

with the support $Γ=∏i=1NΓi⊂ℝN$ constitutes a mapping of the sample space Ω onto the target space Γ. Therefore, a random vector $ξ=(ξ1,ξ2,…,ξN)∈Γ$ denotes an arbitrary point in the parametric space.

According to the Doob–Dynkin lemma [86], the solution $u(t,x;ω)$ can be expressed as $u(t,x;ξ)$, which provides a very useful tool to work in the target space rather than the abstract sample space. Thus, the formulation of SFPDE Eq. (8) can be posed as: Find $u(t,x;ξ):D×Γ→ℝ$ such that $∀t,x∈D$
$0Dtα(ξ)u(t,x;ξ)−∑j=1d kj [ajDxjβj(ξ)+xjDbjβj(ξ)]u(t,x;ξ)=h(t,x)+f(t;ξ)$
(20)

holds ρ-almost surely for $ξ(ω)∈Γ$ and $∀t,x∈D$, subject to proper initial and boundary conditions.

### 3.4 Stochastic Sampling.

We expound the two sampling methods, MCS and PCM to sample from random space and, then propagate the associated uncertainties by computing the statistics of stochastic solutions via postprocessing.

Monte Carlo Sampling: MCS. The general procedure in statistical Monte Carlo sampling is the three following steps.

1. Generating a set of random variables $ξi, i=1,2,…,K$ for a prescribed number of realizations K.

2. Solving the deterministic problem Eq. (20) and obtaining the solution $ui=u(t,x;ξi)$ for each $i=1,2,…,K$.

3. Computing the solution statistics, for example, $E[u]=1M∑i=1Mui.$

We note that steps 1 and 3 are preprocessing and postprocessing steps, respectively. Step 2 requires repetitive simulation of the deterministic counterpart of the problem, which we obtain by developing a Petrov–Galerkin spectral method in the physical domain. Although MCS is relatively easy to implement once a deterministic forward solver is developed, it requires too many samplings for the solution statistics to converge, and yet the extra numerical cost due to nonlocality and memory effect in fractional operators are still remained. In addition, the number of required sampling also grows rapidly as the dimension of the problem increases, resulting in an exhaustively long run time for the statistics to converge.

Probability Collocation Method: PCM. We employ a high-order stochastic discretization in the random space following Refs. [37] and [87] in order to construct a PCM, which yields a high convergence rate with much fewer number of sampling. The idea of PCM is based on polynomial interpolation, however, in the random space. Let $ΘN={ ξi }i=1J$ be a set of prescribed sampling points. By employing the Lagrange interpolation polynomials Li, the polynomial approximation $I$ of the stochastic solution u in the random space can be expressed as
$û(t,x;ξ)=Iu(t,x;ξ)=∑i=1Ju(t,x;ξi)Li(ξ)$
(21)
Therefore, the collocation procedure of solving Eq. (20) to obtain the stochastic solution u is
$R(û(t,x;ξ))|ξi=(Lq(ξ) û(t,x;ξ))−F(t,x;ξ))|ξi=0$
(22)
for $i=1,2,…,J$, where $Lq$ is given in Eq. (11). By using the property of Lagrange interpolants that satisfy the Kronecker delta at the interpolation points, we obtain
$Lq(ξi) u(t,x;ξi))=F(t,x;ξi), i=1,2,…,J$
(23)
subject to proper initial/boundary conditions. Thus, the probabilistic collocation procedure is equivalent to solving $J$ deterministic problems Eq. (23) with conditions (9) and (10). Once the deterministic solutions are obtained at each sampling point, the numerical stochastic solution is interpolated, using Eq. (21) to construct a global approximate $û(t,x;ξ)$. We then obtain the solution statistics as
$E[û]=∫Γû(t,x;ξ) ρ(ξ) dξ, σ[u]=E[û2]−E[û]2$
(24)
The above integrals can be computed efficiently by letting the interpolation/collocation points to be the same as a set of cubature rules $ΘN={ ξi }i=1J$ on the parametric space with integration weights ${wi}i=1J$, which are employed in computing the integral. By property of Kronecker delta of Lagrange interpolant and use of any quadrature rule over the above integral yields
$E[û(t,x:ξ)]≈∑i=1J wi u(t,x;ξi)$
(25)

Choice of Collocation/Interpolation Points. A natural choice of the sampling points is the tensor-product of one-dimensional sets, which is efficient for low-dimensional random spaces. However, in the high-dimensional multivariate case, where $N>6$, the tensor-product interpolation operators are computationally expensive due to the increasing nested summation loops. In addition, the total number of sampling points grows rapidly by the increase of dimension by $JN$, where $J$ is the number of points in each direction.

Another choice that provides an alternative to the more costly full tensor product rule is the isotropic Smolyak sparse grid operator $A(w,N)$ [39,40] with two input parameters dimension size $N$ and the level of grid w. The Smolyak algorithm significantly reduces the total number of sampling points; see Fig. 1 for comparison of A(2, 2), A(4, 2), and A(6, 2) with full tensor product rule for a two-dimensional random spaces. The total number of sampling points for each case is also listed in Table 1. More research has also been devoted to the analysis and construction of Smolyak sparse grids [37,8890].

## 4 Forward Solver

For each realization of random variables in the employed sampling methods, the stochastic model yields a deterministic FPDE, left to be solved in the physical domain. We recall that for every $ξi, i=1,2,…$ in SFPDE Eq. (20), the deterministic problem is recast as
$0Dtαu(t,x)−∑j=1d kj [ajDxjβj+xjDbjβj]u(t,x)=h(t,x)+f(t)$
(26)

subject to the same initial/boundary conditions as Eqs. (9) and (10). In the sequel, we develop a Petrov–Galerkin spectral method to numerically solve the deterministic problem in the physical domain. We also show the well posedness of the deterministic problem in a weak sense and further investigate the stability of the proposed numerical scheme.

### 4.1 Mathematical Framework.

The comprehensive development of mathematical framework and stability analysis of the numerical method are given in Appendix. Here, we only discuss the essential requirement for the rest of the article.

We denote the fractional Sobolev space on $ℝ$ by $Hσ(ℝ), σ≥0$, which is endowed with some proper norms. Let $I=[0,T], Λ1=(a1,b1), Λj=(aj,bj)×Λj−1$ for $j=2,…,d$. We define $X1=H0β12(Λ1)$, and accordingly, $Xj, j=2,…,d$ as
$X2=H0β22((a2,b2);L2(Λ1))∩L2((a2,b2);X1)⋮$
(27)

$Xd=H0βd2((ad,bd);L2(Λd−1))∩L2((ad,bd);Xd−1)$
(28)
and thus, define the “solution space” U and “test space” V as
$U=0lHα2(I;L2(Λd))∩L2(I;Xd),V=0rHα2(I;L2(Λd))∩L2(I;Xd)$
(29)
respectively, where
$0lHα2(I;L2(Λd))= {u | ||u(t,·)||L2(Λd)∈Hα2(I),u|t=0=u|x=aj=u|x=bj=0},0rHα2(I;L2(Λd))= {v | ||v(t,·)||L2(Λd)∈Hα2(I),v|t=T=v|x=aj=v|x=bj=0 }$
For any realization of Eq. (20), we obtain the weak system, that is, the variational form of the deterministic counterpart of the problem, subject to the given initial/boundary conditions, by multiplying the equation with proper test functions and integrating over the whole computational domain $D$. By using Lemmas A.3–A.5, the bilinear form can be written as
$a(u,v)=(0Dtα2 u,tDTα2 v)D−∑j=1dkj[(ajDxjβj2 u, xjDbjβj2 v)D+(xjDbjβj2 u, ajDxjβj2v)D]$
(30)
and thus, by letting U and V be the proper solution/test spaces, the problem reads as: find $u∈U$ such that
$a(u,v)=(f,v)D, ∀v∈V$
(31)

where $f=h(t,x)+f(t)$.

### 4.2 Petrov–Galerkin Spectral Method.

We define the following finite-dimensional solution and test spaces. We first denote the Jacobi polynomials of order n with parameters A, B by $PnA,B(x)$, where $A=B=0$ gives the Legendre polynomial $Pn(x)$ of order n. We also denote the Jacobi polyfractonomial of first kind with parameter τ by $(1)Pn τ(x)=(1+x)τPn−1−τ,τ(x)$. Thus, we employ Legendre polynomials $ϕmj(ξ), j=1,2,…,d$, and Jacobi polyfractonomial of first kind $ψnτ(η)$ [59,60], as the spatial and temporal bases, respectively, given in their corresponding standard domain as
$ϕmj(ξ)=σmj (Pmj+1(ξ)−Pmj−1(ξ))$
(32)

$ψnτ(η)=σn(1)Pn τ(η)=σn(1+η)τPn−1−τ,τ(η)$
(33)
in which $ξ∈[−1,1], mj=1,2,…, σmj=2+(−1)mj, η∈[−1,1], n=1,2,…,$ and $σn=2+(−1)n$. Therefore, by performing affine mappings $η=2t/T−1$ and $ξ=2x−aj/bj−aj−1$ from the computational domain to the standard domain, we construct the solution space UN as
$UN=span {(ψn τ°η)(t) ∏j=1d(ϕmj°ξ)(xj) :n=1,2,…,N, mj=1,2,…,Mj}$
(34)

We note that the choice of temporal and spatial basis functions naturally satisfies the initial and boundary conditions, respectively. The parameter τ in the temporal basis functions plays the role of fine-tuning parameter, which can be chosen properly to capture the singularity of an exact solution.

Moreover, we employ Legendre polynomials $Φrj(ξ),j=1,2,…,d$, and Jacobi polyfractonomial of second kind $Ψkτ(η)$, as the spatial and temporal test functions, respectively, given in their corresponding standard domain as
$Φrj(ξ)=σ̃rj (Prj+1(ξ)−Prj−1(ξ))$
(35)

$Ψkτ(η)=σ̃k(2)Pk τ(η)=σ̃k(1−η)τ Pk−1τ,−τ(η)$
(36)
where $ξ∈[−1,1], rj=1,2,…, σ̃rj=2 (−1)rj+1, η∈[−1,1], k=1,2,…$, and $σ̃k=2 (−1)k+1$. Therefore, by similar affine mapping, we construct the test space VN as
$VN=span {(Ψkτ°η)(t) ∏j=1d(Φrj°ξj)(xj) :k=1,2,…,N, rj=1,2,…,Mj}$
(37)
Thus, since $UN⊂U$ and $VN⊂V$, the problems Eq. (31) read as: find $uN∈UN$ such that
$ah(uN,vN)=l(vN), ∀vN∈VN$
(38)
where $l(vN)=(f,vN)$. The discrete bilinear form $ah(uN,vN)$ can be written as
$ah(uN,vN)=(0Dtα2 uN,tDTα2 vN)D−∑j=1dkj[(ajDxjβj2 uN, xjDbjβj2 vN)D+(xjDbjβj2 uN, ajDxjβj2vN)D]$
(39)
We expand the approximate solution $uN∈UN$, satisfying the discrete bilinear form Eq. (39), in the following form:
$uN(t,x)=∑n=1N∑m1=1M1⋯∑md=1Md ûn,m1,…,md [(ψn τ°η)(t) ∏j=1d(ϕmj°ξ)(xj)]$
(40)
and obtain the corresponding Lyapunov system by substituting Eq. (40) into Eq. (39) by choosing
$vN(t,x)=(Ψkτ°η)(t) ∏j=1d(Φrj°ξj)(xj), k=1,2,…,N, rj=1,2,…,Mj$
Therefore,
$[ST⊗M1⊗M2⋯⊗Md+∑j=1dMT⊗M1⊗⋯⊗Mj−1⊗Sj⊗Mj+1⋯⊗Md] U=F$
(41)
in which ⊗ represents the Kronecker product, F denotes the multidimensional load matrix whose entries are given as
$Fk,r1,…,rd=∫Df(t,x) (Ψk τ°η)(t)∏j=1d(Φrj°ξj)(xj) dD$
(42)

and $U$ is the matrix of unknown coefficients. The matrices ST and MT denote the temporal stiffness and mass matrices, respectively; and the matrices Sj and Mj denote the spatial stiffness and mass matrices, respectively. We obtain the entries of spatial mass matrix Mj analytically and employ proper quadrature rules to accurately compute the entries of other matrices ST, MT, and Sj.

We note that the choices of basis/test functions, employed in developing the PG scheme leads to symmetric mass and stiffness matrices, providing useful properties to further develop a fast solver. The following Theorem 4.1 provides a unified fast solver, developed in terms of the generalized eigensolutions in order to obtain a closed-form solution to the Lyapunov system Eq. (41).

Theorem 4.1 (Unified Fast FPDE Solver [65]). Let${emj,λmj }mj=1Mj$be the set of general Eigen solutions of the spatial stiffness matrix Sj with respect to the mass matrix Mj. Moreover, let${en τ,λnτ }n=1N$be the set of general eigen-solutions of the temporal mass matrix MT with respect to the stiffness matrix ST. Then, the matrix of unknown coefficients$U$is explicitly obtained as
$U=∑n=1N ∑m1=1M1⋯∑md=1Mdκn,m1,…, md en τ ⊗ em1 ⊗⋯⊗ emd$
(43)
where$κn,m1,…, md$is given by
$κn,m1,…, md=( en τ em1⋯ emd)F[(en τT ST en τ) ∏j=1d(emjT Mj emj) ]Λn,m1,…,md$
(44)
in which the numerator represents the standard multidimensional inner product, and$Λn,m1,…,md$is obtained in terms of the eigenvalues of all mass matrices as
$Λn,m1,…,md=[1+λnτ∑j=1d(λmj)]$

## 5 Numerical Results

We investigate the performance of developed numerical methods by considering a couple of numerical simulations. We compare MCS and PCM in random space discretization while using the PG method in the physical domain. We note that by several numerical examples, we make sure that the developed PG method is stable and accurate in solving each deterministic problem; the results are not provided here.

### 5.1 Low-Dimensional Random Inputs.

As the first case, we consider a stochastic fractional initial value problem with random fractional index by letting the diffusion coefficient to be zero, and also ignoring the additional random input and only taking h(t) as the external forcing term. Therefore, we obtain
$0Dtα(ξ)u(t;ξ)=h(t)$
(45)

subject to zero initial condition, where $u(t,ξ):(0,T]×Λ→ℝ$. We let $uext(t)=α/2 t3+α2, h(t)=0Dtα(ξ)uext(t)$ for each realization of α. In this case, by choosing the tuning parameter τ in the temporal basis function to be $α/2$, we can efficiently employ PG numerical scheme and also obtain the exact expectation by rendering the exact solution to be random with similar distribution as the random fractional index. Figure 2 shows the L2-norm convergence rate of MCS and PCM in comparison to solution expectation with $Eext[u]=E[uext]$. The results confirm converges rate of 0.5 for MCS, while in PCM, the statistics of solution converges accurately very fast, using only a few numbers of realizations. In this example, by ignoring the additional random input to the system, we take the advantage of having the exact random solution to be available.

As another example, we also consider Eq. (45) with additional random input, expanded by KL expansion with M =4, as
$0Dtα(ξ)u(t;ξ)=h(t)+∑k=1M ak sin(2kπ tT) ξk$
(46)

with two cases $h(t)=t2$ and $h(t)=sin(πt)$. Figure 3 shows the mean value and variance of solution for 104 sampling of MCS compared to 625 realizations in PCM.

Moreover, we consider (1 + 1)-D one-sided SFPDE given in Eq. (20), where d =1 and the diffusion coefficient is kl. We ignore the additional random input and consider h(t, x) as the only external forcing term. Therefore, we obtain
$0Dtα(ξ1)u(t,x;ξ)−kl−1Dx β(ξ2)u(t,x;ξ)=h(t,x)$
(47)

subject to zero initial/boundary conditions, where $u(t,x;ξ):(0,T]×(−1,1)×Λ→ℝ$, and the only random variables are the fractional indices α and β. We let $uext(t,x)=t3+τ ((1+x)3+μ−1/2(1+x)4+μ)$, and choose $τ=α/2$ and $μ=β/2$. For each realization of α and β, we obtain the force function h(t, x) by substituting the corresponding uext to Eq. (47). Defining $Eext[u]=E[uext]$, Fig. 4 shows the L2-norm convergence of solution expectation as compared to the exact expectation. We observe that PCM converges accurately with only a few number of realizations.

Considering additional random input, expanded by KL expansion with M =4, the problem can be recast as
$0Dtα(ξ)u(t,x;ξ)−kl−1Dx β(ξ)u(t,x;ξ)=h(t,x)+∑k=1M ak sin(2kπ tT) ξk$
(48)

subject to zero initial/boundary conditions. Figure 5 shows the mean value of the solution for MCS and PCM at different times.

Remark 5.1. We note that generally use of the sparse grid operators in obtaining solution statistics is more effective when the dimension of the random space is higher than 6. Thus, in the numerical examples for low-dimensional random inputs, we employ the easy-to-implement tensor product nodal sets.▪

### 5.2 Moderate- to High-Dimensional Random Inputs.

We render the problem with a higher number of terms in KL expansion of random inputs in Eq. (48) by choosing M =10 and M =20. This yields the dimension of random space $N=12$ and $N=22$, respectively. As mentioned in Remark 5.1, in the case of high-dimensional random space constructing a grid based on tensor product rule results in very expensive computation of solution statistics due to exhaustive increase of forward solver instruction. Table 1 shows the comparison between different levels of the Smolyak algorithm and tensor product rule. Therefore, to obtain the solution statistics, we employ the Smolyak sparse grid generator in the developed PCM. For each cases of KL expansion, we generate the sparse grid on two levels w =1 and w =2, that is, A(1, 12), A(2, 12), A(1, 22), and A(2, 22), where we let the higher resolution case be a benchmark value to the solution statistics, based on which we compute and normalize the error. We observe that for both cases $N=12$ and $N=22$, the normalized error in computing the expectation and standard deviation of the solution are of orders $O(10−7)$ and $O(10−3)$, respectively.

## 6 Summary and Discussion

We developed a mathematical framework to numerically quantify the solution uncertainty of a stochastic FPDE, associated with the randomness of model parameters. The stochastic FPDE is reformulated by rendering the problem with random fractional indices, subject to additional random noise. We used the truncated Karhunen–Loéve expansion to parametrize the additive noise. Then, by employing a nonintrusive PCM, we propagated the associated randomness to the system response, by using Smolyak sparse grid generator to construct the set of sample points in the random space. We also formulated a forward solver to simulate the deterministic counterpart of the stochastic problem for each realization of random variables. We showed that the deterministic problem is mathematically well posed in a weak sense. Furthermore, by employing Jacobi polyfractonomials and Legendre polynomials as the temporal and spatial basis/test functions, respectively, we developed a Petrov–Galerkin spectral method to solve the deterministic problem in the physical domain. We also proved that the inf-sup condition holds for the proposed numerical scheme, and thus, it is stable. By considering several numerical examples with low- to high-dimensional random spaces, we examined the performance of our stochastic discretization. We showed that in each case, PCM converges very fast to a very high level of accuracy with very few number of sampling.

## Funding Data

• ARO Young Investigator Program Award (W911NF-19-1-0444; Funder ID: 10.13039/100000183).

• National Science Foundation Award (DMS-1923201; Funder ID: 10.13039/100000001).

• MURI/ARO (W911NF-15-1-0562; Funder ID: 10.13039/100014036).

• AFOSR Young Investigator Program Award (FA9550-17-1-0150; Funder ID: 10.13039/100000181).

### Appendix: Mathematical Framework of Deterministic Solver

Here, we extensively discuss the mathematical framework of the developed method in detail.

##### A.1 Mathematical Framework.
We define the useful functional spaces and their associated norms [54,69]. By $Hσ(ℝ)= {u(t)|u∈L2(ℝ); (1+|ω|2)σ/2F(u)(ω)∈L2(ℝ)}, σ≥0$, we denote the fractional Sobolev space on $ℝ$, endowed with the norm $||u||Hℝσ=||(1+|ω|2)σ/2F(u)(ω)||L2(ℝ)$, where $F(u)$ represents the Fourier transform of u. Subsequently, we denote by $Hσ(Λ)={u∈L2(Λ) | ∃ ũ∈Hσ(ℝ) s.t. ũ|Λ=u }, σ≥0$, the fractional Sobolev space on any finite closed interval, for example, $Λ=(a,b)$, with norm $||u||Hσ(Λ)=infũ∈Hℝσ, ũ|Λ=u ||ũ||Hσ(ℝ)$. We define the following useful norms as:
$||·||lHσ(Λ)=(||aDxσ (·)||L2(Λ)2+||·||L2(Λ)2)12,||·||rHσ(Λ)=(||xDbσ (·)||L2(Λ)2+||·||L2(Λ)2)12,||·||cHσ(Λ)=(||xDbσ (·)||L2(Λ)2+||aDxσ (·)||L2(Λ)2+||·||L2(Λ)2)12,$

where the equivalence of $||·||lHσ(Λ)$ and $||·||rHσ(Λ)$ are shown in Refs. [54], [55], and [91].

Lemma A.1. Let$σ≥0$and$σ≠n−12$. Then, the norms$||·||lHσ(Λ)$and$||·||rHσ(Λ)$are equivalent to$||·||cHσ(Λ)$.

We also define $C0∞(Λ)$ as the space of smooth functions with compact support in (a, b). We denote by $lH0σ(Λ), rH0σ(Λ)$, and $cH0σ(Λ)$ as the closure of $C0∞(Λ)$ with respect to the norms $||·||lHσ(Λ), ||·||rHσ(Λ)$, and $||·||cHσ(Λ)$. It is shown in Refs. [55] and [91] that these Sobolev spaces are equal and their seminorms are also equivalent to $|·|Hσ(Λ)*=|(aDxσ (·),xDbσ (·))|Λ1/2$. Therefore, we can prove that $|(aDxσ u,xDbσ v)Λ|≥β |u|lHσ(Λ) |v|rHσ(Λ)$ and $|(xDbσ u,aDxσ v)Λ |≥β |u|rHσ(Λ) |v|lHσ(Λ)$, in which β is a positive constant.

Moreover, by letting $0C∞(I)$ and $C0∞(I)$ be the space of smooth functions with compact support in $(0,T]$ and $[0,T)$, respectively, we define $lHs(I)$ and $rHs(I)$ as the closure of $0C∞(I)$ and $C0∞(I)$ with respect to the norms $||·||lHs(I)$ and $||·||rHs(I)$. Other equivalent useful seminorms associated with $Hs(I)$ are also introduced in Refs. [54] and [91], as $|·|lHs(I)=||0Dts(·)||L2(I), |·|rHs(I)=||tDTs(·)||L2(I), |·|Hs(I)*=|(0Dts(·),tDTs(·))I|1/2$, where $|·|Hs(I)*≡|·|lHs(I)1/2 |·|rHs(I)1/2$.

Borrowing definitions from Ref. [65], we define the following spaces, which we use later in the construction of the corresponding solution and test spaces of our problem. Thus, by letting $Λ1=(a1,b1), Λj=(aj,bj)×Λj−1$ for $j=2,…,d$, we define $X1=H0β12(Λ1)$, which is associated with the norm $||·||cHβ12(Λ1)$, and accordingly, $Xj, j=2,…,d$ as
$X2=H0β22((a2,b2);L2(Λ1))∩L2((a2,b2);X1)⋮$
(A1)

$Xd=H0βd2((ad,bd);L2(Λd−1))∩L2((ad,bd);Xd−1)$
(A2)
associated with norms
$||·||Xj= {||·||H0βj2((aj,bj);L2(Λj−1))2+||·||L2((aj,bj);Xj−1)2}12$

for $j=2,3,…,d$.

Lemma A.2. Let$βj≥0$and$βj≠n−1/2$. Then, for$j=1,2,…,d$
$||·||Xj≡{∑i=1j(||xiDbiβi/2 (·)||L2(Λj)2+||aiDxiβi/2 (·)||L2(Λj)2)+||·||L2(Λj)2}12$
###### Solution and Test Spaces.
We define the “solution space” U and “test space” V, respectively, as
$U=0lHα2(I;L2(Λd))∩L2(I;Xd),V=0rHα2(I;L2(Λd))∩L2(I;Xd),$
endowed with norms
$||u||U={||u||lHα2(I;L2(Λd))2+||u||L2(I;Xd)2}12,||v||V={||v||rHα2(I;L2(Λd))2+||v||L2(I;Xd)2}12$
(A3)
where $I=[0,T]$, and
$0lHα2(I;L2(Λd))= {u | ||u(t,·)||L2(Λd)∈Hα2(I),u|t=0=u|x=aj=u|x=bj=0},0rHα2(I;L2(Λd))= {v | ||v(t,·)||L2(Λd)∈Hα2(I),v|t=T=v|x=aj=v|x=bj=0}$
equipped with norms $||u||lHα/2(I;L2(Λd))$ and $||u||rHα/2(I;L2(Λd))$, respectively. We can show that these norms take the following forms:
$||u||lHα2(I;L2(Λd))=|| ||u(t,·)||L2(Λd) ||lHα2(I)=(||0Dtα2 (u)||L2(Ω)2+||u||L2(Ω)2)12,||u||rHα2(I;L2(Λd))=|| ||u(t,·)||L2(Λd) ||rHα2(I)=(||tDTα2 (u)||L2(Ω)2+||u||L2(Ω)2)12$
(A4)
Also, using Lemma A.2, we can show that
$||u||L2(I;Xd)=|| ||u(t,.)||Xd ||L2(I)={||u||L2(Ω)2+∑j=1d(||xjDbjβj2 (u)||L2(Ω)2+||ajDxjβj2 (u)||L2(Ω)2)}12$
(A5)
Therefore, Eq. (A3) can be written as
$||u||U={||u||L2(Ω)2+||0Dtα2 (u)||L2(Ω)2+∑j=1d(||xjDbjβj2 (u)||L2(Ω)2+||ajDxjβj2 (u)||L2(Ω)2)}12$
(A6)

$||v||V={||v||L2(Ω)2+||tDTα2 (v)||L2(Ω)2 +∑j=1d(||xjDbjβj2 (v)||L2(Ω)2+||ajDxjβj2 (v)||L2(Ω)2)}12$
(A7)
##### A.2 Weak Formulation.

The following lemmas help us obtain the weak formulation of the deterministic problem in the physical domain and construct the numerical scheme.

Lemma A.3 [54]. For all$α∈(0,1)$, if$u∈H1([0,T])$such that$u(0)=0$, and$v∈Hα/2([0,T])$, then$(0Dt αu,v)Ω=( 0Dt α/2u , tDT α/2v )Ω$, where$(·,·)Ω$represents the standard inner product in$Ω=[0,T]$.

Lemma A.4 [69]. Let$1<β<2$, a and b be arbitrary finite or infinite real numbers. Assume$u∈Hβ(a,b)$such that u(a) = 0, also$xDbβ/2v$is integrable in (a, b) such that v(b) = 0. Then,$(aDxβu , v)=(aDxβ/2u , xDbβ/2v)$.

Lemma A.5. Let$1<βj<2$for$j=1,2,…,d$, and$u,v∈Xd$. Then,
$(ajDxjβju,v)Λd=(ajDxjβj2u,xjDbjβj2v)Λd,(xjDbjβju,v)Λd=(xjDbjβj2u,ajDxjβj2v)Λd$
For any realization of Eq. (20), we obtain the weak system, that is, the variational form of the deterministic counterpart of the problem, subject to the given initial/boundary conditions, by multiplying the equation with proper test functions and integrate over the whole computational domain $D$. Using Lemmas A.3–A.5, the bilinear form can be written as
$a(u,v)=(0Dtα2 u,tDTα2 v)D−∑j=1dkj[(ajDxjβj2 u, xjDbjβj2 v)D+(xjDbjβj2 u, ajDxjβj2v)D]$
(A8)
and thus, by letting U and V be the proper solution/test spaces, the problem reads as: find $u∈U$ such that
$a(u,v)=(f,v)D, ∀v∈V$
(A9)

where $f=h(t,x)+f(t).$

##### A.3 Stability Analysis.

We show the well posedness of the deterministic problem and prove the stability of the proposed PG scheme.

Lemma A.6. Let$α∈(0,1), Ω=I×Λd$, and$u∈0lHα/2(I;L2(Λd))$. Then,
$|(0Dtα/2u,tDTα/2v)Ω|≡||u||lHα/2(I;L2(Λd)) ||v||rHα/2(I;L2(Λd)), ∀v∈0rHα/2(I;L2(Λd))$
Moreover,
$|(adDxdβd/2u,xdDbdβd/2v)Λd|≡|u|cHβd/2((ad,bd);L2(Λd−1)) |v|cHβd/2((ad,bd);L2(Λd−1))$
(A10)
and
$|(xdDbdβd/2u,adDxdβd/2v)Λd|≡|u|cHβd/2((ad,bd);L2(Λd−1)) |v|cHβd/2((ad,bd);L2(Λd−1))$
(A11)
Lemma A.7 (Continuity). The bilinear form Eq. (A8) is continuous, that is,
$∀u∈U, ∃ γ>0, s.t. |a(u,v)|≤γ ||u||U ||v||V, ∀v∈V$
(A12)

Proof. The proof directly concludes from Eq. (A10) and Lemma A.6.▪

Theorem A.8 (Stability). The following inf-sup condition holds for the bilinear form Eq. (A8), that is,
$infu≠0∈U supv≠0∈V|a(u,v)| ||v||V ||u||U≥γ>0$
(A13)

where$Ω=I×Λd$and$supu∈U |a(u,v)|>0$.

Theorem A.9 (well-posedness). For all$0<α<1$and$1<βj<2$, and$j=1,…,d$, there exists a unique solution to Eq. (A9), continuously dependent on f, where f belongs to the dual space of U.

Proof. Lemmas A.7 (continuity) and A.8 (stability) yield the well posedness of weak form Eq. (A9) in (1 + d)-dimension due to the generalized Babuška–Lax–Milgram theorem.▪

Since the defined basis and test spaces are Hilbert spaces, and $UN⊂U$ and $VN⊂V$, we can prove that the developed Petrov–Galerkin spectral method is stable and the following condition holds:
$infuN≠0∈UN supv≠0∈VN|a(uN,vN)|||vN||V ||uN||U≥γ>0$
(A14)

with $γ>0$ and independent of N, where $supuN∈UN|a(uN,vN)|>0$.

## References

References
1.
Zhang
,
Y.
,
Sun
,
H.
,
Stowell
,
H. H.
,
Zayernouri
,
M.
, and
Hansen
,
S. E.
,
2017
, “
A Review of Applications of Fractional Calculus in Earth System Dynamics
,”
Chaos, Solitons Fractals
,
102
, pp.
29
46
.10.1016/j.chaos.2017.03.051
2.
Jaishankar
,
A.
, and
McKinley
,
G. H.
,
2014
, “
A Fractional K-BKZ Constitutive Formulation for Describing the Nonlinear Rheology of Multiscale Complex Fluids
,”
J. Rheol.
,
58
(
6
), pp.
1751
1788
.10.1122/1.4892114
3.
Jha
,
R.
,
Kaw
,
P. K.
,
Kulkarni
,
D. R.
,
Parikh
,
J. C.
, and
,
2003
, “
Evidence of Lévy Stable Process in Tokamak Edge Turbulence
,”
Phys. Plasmas
,
10
(
3
), pp.
699
704
.10.1063/1.1541607
4.
del Castillo-Negrete
,
D.
,
Carreras
,
B. A.
, and
Lynch
,
V. E.
,
2004
, “
Fractional Diffusion in Plasma Turbulence
,”
Phys. Plasmas
,
11
(
8
), pp.
3854
3864
.10.1063/1.1767097
5.
Jaishankar
,
A.
, and
McKinley
,
G. H.
,
2013
, “
Power-Law Rheology in the Bulk and at the Interface: Quasi-Properties and Fractional Constitutive Equations
,”
Proc. R. Soc. A
,
469
(
2149
), p.
20120284
.10.1098/rspa.2012.0284
6.
Naghibolhosseini
,
M.
,
2015
, “
Estimation of Outer-Middle Ear Transmission Using DPOAEs and Fractional-Order Modeling of Human Middle Ear
,” Ph.D. thesis, City University of New York, New York.
7.
Naghibolhosseini
,
M.
, and
Long
,
G. R.
,
2017
, “
Fractional-Order Modelling and Simulation of Human Ear
,”
Int. J. Comput. Math.
,
95
(
6–7
), pp.
1257
1273
.10.1080/00207160.2017.1404038
8.
Samiee
,
M.
,
Akhavan-Safaei
,
A.
, and
Zayernouri
,
M.
,
2019
, “
A Fractional Subgrid-Scale Model for Turbulent Flows: Theoretical Formulation and a Priori Study
,” preprint arXiv: 1909.09943.
9.
Meral
,
F. C.
,
Royston
,
T. J.
, and
Magin
,
R.
,
2010
, “
Fractional Calculus in Viscoelasticity: An Experimental Study
,”
Commun. Nonlinear Sci. Numer. Simul.
,
15
(
4
), pp.
939
945
.10.1016/j.cnsns.2009.05.004
10.
Cullen
,
A. C.
, and
Frey
,
H. C.
,
1999
,
Probabilistic Techniques in Exposure Assessment: A Handbook for Dealing With Variability and Uncertainty in Models and Inputs
,
, New York, p.
336
.
11.
Roy
,
C. J.
, and
Oberkampf
,
W. L.
,
2011
, “
A Comprehensive Framework for Verification, Validation, and Uncertainty Quantification in Scientific Computing
,”
Comput. Methods Appl. Mech. Eng.
,
200
(
25–28
), pp.
2131
2144
.10.1016/j.cma.2011.03.016
12.
Mullins
,
J.
, and
,
S.
,
2016
, “
Bayesian Uncertainty Integration for Model Calibration, Validation, and Prediction
,”
ASME J. Verif., Validation Uncertainty Quantif.
,
1
(
1
), p.
011006
.10.1115/1.4032371
13.
Der Kiureghian
,
A.
, and
Ditlevsen
,
O.
,
2009
, “
Aleatory or Epistemic? Does It Matter?
,”
Struct. Saf.
,
31
(
2
), pp.
105
112
.10.1016/j.strusafe.2008.06.020
14.
Benson
,
D. A.
,
Wheatcraft
,
S. W.
, and
Meerschaert
,
M. M.
,
2000
, “
Application of a Fractional Advection-Dispersion Equation
,”
Water Resour. Res.
,
36
(
6
), pp.
1403
1412
.10.1029/2000WR900031
15.
Baeumer
,
B.
,
Benson
,
D. A.
,
Meerschaert
,
M. M.
, and
Wheatcraft
,
S. W.
,
2001
, “
Subordinated Advection-Dispersion Equation for Contaminant Transport
,”
Water Resour. Res.
,
37
(
6
), pp.
1543
1550
.
16.
Fishman
,
G. S.
,
1996
,
Monte Carlo: Concepts, Algorithms, and Applications
, Springer, New York.
17.
Del Moral
,
P.
,
Doucet
,
A.
, and
Jasra
,
A.
,
2006
, “
Sequential Monte Carlo Samplers
,”
J. R. Stat. Soc.: Ser. B
,
68
(
3
), pp.
411
436
.10.1111/j.1467-9868.2006.00553.x
18.
Beskos
,
A.
,
Jasra
,
A.
,
Law
,
K.
,
Tempone
,
R.
, and
Zhou
,
Y.
,
2017
, “
Multilevel Sequential Monte Carlo Samplers
,”
Stochastic Processes Their Appl.
,
127
(
5
), pp.
1417
1440
.10.1016/j.spa.2016.08.004
19.
Jasra
,
A.
,
Law
,
K. J.
, and
Zhou
,
Y.
,
2016
, “
Forward and Inverse Uncertainty Quantification Using Multilevel Monte Carlo Algorithms for an Elliptic Nonlocal Equation
,”
Int. J. Uncertainty Quantif.
,
6
(
6
), pp.
501
514
.10.1615/Int.J.UncertaintyQuantification.2016018661
20.
Owen
,
N. E.
,
Challenor
,
P.
,
Menon
,
P. P.
, and
Bennani
,
S.
,
2017
, “
Comparison of Surrogate-Based Uncertainty Quantification Methods for Computationally Expensive Simulators
,”
SIAM/ASA J. Uncertainty Quantif.
,
5
(
1
), pp.
403
435
.10.1137/15M1046812
21.
Ghanem
,
R. G.
, and
Spanos
,
P. D.
,
2003
,
Stochastic Finite Elements: A Spectral Approach
,
Courier Corporation
, Springer, New York.
22.
Xiu
,
D.
, and
,
G. E.
,
2002
, “
Modeling Uncertainty in Steady State Diffusion Problems Via Generalized Polynomial Chaos
,”
Comput. Methods Appl. Mech. Eng.
,
191
(
43
), pp.
4927
4948
.10.1016/S0045-7825(02)00421-8
23.
Xiu
,
D.
, and
,
G. E.
,
2002
, “
The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations
,”
SIAM J. Sci. Comput.
,
24
(
2
), pp.
619
644
.10.1137/S1064827501387826
24.
Knio
,
O. M.
, and
Le Maître
,
O. P.
,
2006
, “
Uncertainty Propagation in CFD Using Polynomial Chaos Decomposition
,”
Fluid Dyn. Res.
,
38
(
9
), pp.
616
640
.10.1016/j.fluiddyn.2005.12.003
25.
Najm
,
H. N.
,
2009
, “
Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics
,”
Annu. Rev. Fluid Mech.
,
41
(
1
), pp.
35
52
.10.1146/annurev.fluid.010908.165248
26.
Van Dam
,
N.
, and
Rutland
,
C.
,
2016
, “
Uncertainty Quantification of Large-Eddy Spray Simulations
,”
ASME J. Verif., Validation Uncertainty Quantif.
,
1
(
2
), p.
021006
.10.1115/1.4032196
27.
Babuska
,
I.
,
Tempone
,
R.
, and
Zouraris
,
G. E.
,
2004
, “
Galerkin Finite Element Approximations of Stochastic Elliptic Partial Differential Equations
,”
SIAM J. Numer. Anal.
,
42
(
2
), pp.
800
825
.10.1137/S0036142902418680
28.
Babuška
,
I.
,
Tempone
,
R.
, and
Zouraris
,
G. E.
,
2005
, “
Solving Elliptic Boundary Value Problems With Uncertain Coefficients by the Finite Element Method: The Stochastic Formulation
,”
Comput. Methods Appl. Mech. Eng.
,
194
(
12–16
), pp.
1251
1294
.10.1016/j.cma.2004.02.026
29.
Le Maître
,
O. P.
,
Knio
,
O. M.
,
Najm
,
H. N.
, and
Ghanem
,
R. G.
,
2004
, “
Uncertainty Propagation Using Wiener–Haar Expansions
,”
J. Comput. Phys.
,
197
(
1
), pp.
28
57
.10.1016/j.jcp.2003.11.033
30.
Le Maître
,
O. P.
,
Najm
,
H. N.
,
Ghanem
,
R. G.
, and
Knio
,
O. M.
,
2004
, “
Multi-Resolution Analysis of Wiener-Type Uncertainty Propagation Schemes
,”
J. Comput. Phys.
,
197
(
2
), pp.
502
531
.10.1016/j.jcp.2003.12.020
31.
Schuss
,
Z.
,
1980
, “
Singular Perturbation Methods in Stochastic Differential Equations of Mathematical Physics
,”
SIAM Rev.
,
22
(
2
), pp.
119
155
.10.1137/1022024
32.
Babuška
,
I.
, and
Chatzipantelidis
,
P.
,
2002
, “
On Solving Elliptic Stochastic Partial Differential Equations
,”
Comput. Methods Appl. Mech. Eng.
,
191
(
37–38
), pp.
4093
4122
.10.1016/S0045-7825(02)00354-7
33.
Todor
,
R. A.
,
2005
, “
Sparse Perturbation Algorithms for Elliptic PDE's With Stochastic Data
,” Ph.D. thesis, ETH Zurich, Zurich, Switzerland.
34.
Winter
,
C. L.
, and
Tartakovsky
,
D. M.
,
2002
, “
Groundwater Flow in Heterogeneous Composite Aquifers
,”
Water Resour. Res.
,
38
(
8
), pp.
23
1
23-11
.10.1029/2001WR000450
35.
Liu
,
W. K.
,
Belytschko
,
T.
, and
Mani
,
A.
,
1986
, “
Probabilistic Finite Elements for Nonlinear Structural Dynamics
,”
Comput. Methods Appl. Mech. Eng.
,
56
(
1
), pp.
61
81
.10.1016/0045-7825(86)90136-2
36.
Liu
,
W. K.
,
Belytschko
,
T.
, and
Mani
,
A.
,
1986
, “
Random Field Finite Elements
,”
Int. J. Numer. Methods Eng.
,
23
(
10
), pp.
1831
1845
.10.1002/nme.1620231004
37.
Xiu
,
D.
, and
Hesthaven
,
J. S.
,
2005
, “
High-Order Collocation Methods for Differential Equations With Random Inputs
,”
SIAM J. Sci. Comput.
,
27
(
3
), pp.
1118
1139
.10.1137/040615201
38.
Babuška
,
I.
,
Nobile
,
F.
, and
Tempone
,
R.
,
2007
, “
A Stochastic Collocation Method for Elliptic Partial Differential Equations With Random Input Data
,”
SIAM J. Numer. Anal.
,
45
(
3
), pp.
1005
1034
.10.1137/050645142
39.
Nobile
,
F.
,
Tempone
,
R.
, and
Webster
,
C. G.
,
2008
, “
A Sparse Grid Stochastic Collocation Method for Partial Differential Equations With Random Input Data
,”
SIAM J. Numer. Anal.
,
46
(
5
), pp.
2309
2345
.10.1137/060663660
40.
Smolyak
,
S.
,
1963
, “
Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions
,”
Sov. Math. Dokl.
,
4
(
5
), pp.
240
243
.http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=dan&paperid=27586&option_lang=eng
41.
Gorenflo
,
R.
,
Mainardi
,
F.
,
Moretti
,
D.
, and
,
P.
,
2002
, “
Time Fractional Diffusion: A Discrete Random Walk Approach
,”
Nonlinear Dyn.
,
29
(
1–4
), pp.
129
143
.10.1023/A:1016547232119
42.
Sun
,
Z.
, and
Wu
,
X.
,
2006
, “
A Fully Discrete Difference Scheme for a Diffusion-Wave System
,”
Appl. Numer. Math.
,
56
(
2
), pp.
193
209
.10.1016/j.apnum.2005.03.003
43.
Lin
,
Y.
, and
Xu
,
C.
,
2007
, “
Finite Difference/Spectral Approximations for the Time-Fractional Diffusion Equation
,”
J. Comput. Phys.
,
225
(
2
), pp.
1533
1552
.10.1016/j.jcp.2007.02.001
44.
Wang
,
H.
,
Wang
,
K.
, and
Sircar
,
T.
,
2010
, “
A Direct O(N log2N) Finite Difference Method for Fractional Diffusion Equations
,”
J. Comput. Phys.
,
229
(
21
), pp.
8095
8104
.10.1016/j.jcp.2010.07.011
45.
Wang
,
K.
, and
Wang
,
H.
,
2011
, “
A Fast Characteristic Finite Difference Method for Fractional Advection–Diffusion Equations
,”
,
34
(
7
), pp.
810
816
46.
Cao
,
J.
, and
Xu
,
C.
,
2013
, “
A High Order Schema for the Numerical Solution of the Fractional Ordinary Differential Equations
,”
J. Comput. Phys.
,
238
(
1
), pp.
154
168
.10.1016/j.jcp.2012.12.013
47.
Zeng
,
F.
,
Li
,
C.
,
Liu
,
F.
, and
Turner
,
I.
,
2015
, “
Numerical Algorithms for Time-Fractional Subdiffusion Equation With Second-Order Accuracy
,”
SIAM J. Sci. Comput.
,
37
(
1
), pp.
A55
A78
.10.1137/14096390X
48.
Zayernouri
,
M.
, and
Matzavinos
,
A.
,
2016
, “
Fractional Adams-Bashforth/Moulton Methods: An Application to the Fractional Keller–Segel Chemotaxis System
,”
J. Comput. Phys.
,
317
, pp.
1
14
.10.1016/j.jcp.2016.04.041
49.
,
G.
, and
Sherwin
,
S.
,
2013
,
Spectral/hp Element Methods for Computational Fluid Dynamics
,
Oxford University Press
, New York.
50.
Canuto
,
C.
,
Hussaini
,
M. Y.
,
Quarteroni
,
A.
, and
Zang
,
T. A.
,
2006
,
Spectral Methods
,
Springer
, New York.
51.
Rawashdeh
,
E.
,
2006
, “
Numerical Solution of Fractional Integro-Differential Equations by Collocation Method
,”
Appl. Math. Comput.
,
176
(
1
), pp.
1
6
.10.1016/j.amc.2005.09.059
52.
,
M.
,
2011
, “
On the Numerical Solutions for the Fractional Diffusion Equation
,”
Commun. Nonlinear Sci. Numer. Simul.
,
16
(
6
), pp.
2535
2542
.10.1016/j.cnsns.2010.09.007
53.
,
M.
, and
Hendy
,
A.
,
2012
, “
The Approximate and Exact Solutions of the Fractional-Order Delay Differential Equations Using Legendre Pseudospectral Method
,”
Int. J. Pure Appl. Math.
,
74
(
3
), pp.
287
297
.https://ijpam.eu/contents/2012-74-3/1/
54.
Li
,
X.
, and
Xu
,
C.
,
2009
, “
A Space-Time Spectral Method for the Time Fractional Diffusion Equation
,”
SIAM J. Numer. Anal.
,
47
(
3
), pp.
2108
2131
.10.1137/080718942
55.
Li
,
X.
, and
Xu
,
C.
,
2010
, “
Existence and Uniqueness of the Weak Solution of the Space-Time Fractional Diffusion Equation and a Spectral Method Approximation
,”
Commun. Comput. Phys.
,
8
(
5
), pp.
1016
1051
.10.4208/cicp.020709.221209a
56.
Chen
,
S.
,
Shen
,
J.
, and
Wang
,
L.
,
2014
, “
Generalized Jacobi Functions and Their Applications to Fractional Differential Equations
,”
Math. Comput.
85
, pp.
1603
1638
.https://www.ams.org/journals/mcom/2016-85-300/S0025-5718-2015-03035-X/home.html
57.
Wang
,
H.
, and
Zhang
,
X.
,
2015
, “
A High-Accuracy Preserving Spectral Galerkin Method for the Dirichlet Boundary-Value Problem of Variable-Coefficient Conservative Fractional Diffusion Equations
,”
J. Comput. Phys.
,
281
, pp.
67
81
.10.1016/j.jcp.2014.10.018
58.
Bhrawy
,
A. H.
,
Doha
,
E. H.
,
Baleanu
,
D.
, and
Ezz-Eldien
,
S. S.
,
2015
, “
A Spectral Tau Algorithm Based on Jacobi Operational Matrix for Numerical Solution of Time Fractional Diffusion-Wave Equations
,”
J. Comput. Phys.
,
293
, pp.
142
156
.10.1016/j.jcp.2014.03.039
59.
Zayernouri
,
M.
, and
,
G. E.
,
2013
, “
Fractional Sturm–Liouville Eigen-Problems: Theory and Numerical Approximations
,”
J. Comput. Phys
.,
252
, pp.
495
517
.https://www.sciencedirect.com/science/article/pii/S0021999113004610
60.
Zayernouri
,
M.
,
Ainsworth
,
M.
, and
,
G. E.
,
2015
, “
Tempered Fractional Sturm–Liouville Eigenproblems
,”
SIAM J. Sci. Comput.
,
37
(
4
), pp.
A1777
A1800
.10.1137/140985536
61.
Lischke
,
A.
,
Zayernouri
,
M.
, and
Zhang
,
Z.
, “
Spectral and Spectral Element Methods for Fractional Advection–Diffusion–Reaction Equations
,”
Handbook of Fractional Calculus With Applications, Numerical Methods
, Walter de Gruyter GmbH & Co KG, Berlin.
62.
Lischke
,
A.
,
Zayernouri
,
M.
, and
,
G. E.
,
2017
, “
A Petrov–Galerkin Spectral Method of Linear Complexity for Fractional Multiterm ODEs on the Half Line
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
A922
A946
.10.1137/17M1113060
63.
Suzuki
,
J. L.
,
Zayernouri
,
M.
,
Bittencourt
,
M. L.
, and
,
G. E.
,
2016
, “
Fractional-Order Uniaxial Visco-Elasto-Plastic Models for Structural Analysis
,”
Comput. Methods Appl. Mech. Eng.
,
308
, pp.
443
467
.10.1016/j.cma.2016.05.030
64.
Suzuki
,
J. L.
, and
Zayernouri
,
M.
,
2018
, “
An Automated Singularity-Capturing Scheme for Fractional Differential Equations
,” preprint arXiv:1810.12219.
65.
Samiee
,
M.
,
Zayernouri
,
M.
, and
Meerschaert
,
M. M.
,
2019
, “
A Unified Spectral Method for PDEs With Two-Sided Derivatives; Part I: A Fast Solver
,”
J. Comput. Phys.
,
385
, pp.
225
243
.10.1016/j.jcp.2018.02.014
66.
Samiee
,
M.
,
Kharazmi
,
E.
, and
Zayernouri
,
M.
,
2017
, “
Fast Spectral Methods for Temporally-Distributed Fractional PDEs
,” Spectral and High Order Methods for Partial Differential Equations (ICOSAHOM 2016) (Lecture Notes in Computational Science and Engineering, Vol. 119), Springer, Cham, Switzerland, pp.
651
667
.
67.
Samiee
,
M.
,
Zayernouri
,
M.
, and
Meerschaert
,
M. M.
,
2019
, “
A Unified Spectral Method for FPDEs With Two-Sided Derivatives—Part II: Stability, and Error Analysis
,”
J. Comput. Phys.
,
385
, pp.
244
261
.10.1016/j.jcp.2018.07.041
68.
Samiee
,
M.
,
Kharazmi
,
E.
,
Zayernouri
,
M.
, and
Meerschaert
,
M. M.
,
2018
, “
Petrov–Galerkin Method for Fully Distributed-Order Fractional Partial Differential Equations
,” preprint arXiv:1805.08242.
69.
Kharazmi
,
E.
,
Zayernouri
,
M.
, and
,
G. E.
,
2017
, “
Petrov–Galerkin and Spectral Collocation Methods for Distributed Order Differential Equations
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
A1003
A1037
.10.1137/16M1073121
70.
Kharazmi
,
E.
,
Zayernouri
,
M.
, and
,
G. E.
, “
A Petrov–Galerkin Spectral Element Method for Fractional Elliptic Problems
,”
Comput. Methods Appl. Mech. Eng.
,
324
, pp.
512
536
.10.1016/j.cma.2017.06.006
71.
Kharazmi
,
E.
, and
Zayernouri
,
M.
,
2018
, “
Fractional Pseudo-Spectral Methods for Distributed-Order Fractional PDEs
,”
Int. J. Comput. Math.
,
95
(
6–7
), pp.
1340
1361
.10.1080/00207160.2017.1421949
72.
Kharazmi
,
E.
, and
Zayernouri
,
M.
,
2019
, “
Fractional Sensitivity Equation Method: Application to Fractional Model Construction
,”
J. Sci. Comput.
,
80
(
1
), pp.
110
140
.10.1007/s10915-019-00935-0
73.
Rizzi
,
F.
,
Najm
,
H. N.
,
Debusschere
,
B. J.
,
Sargsyan
,
K.
,
Salloum
,
M.
,
,
H.
, and
Knio
,
O. M.
,
2012
, “
Uncertainty Quantification in MD Simulations. Part I: Forward Propagation
,”
Multiscale Model. Simul.
,
10
(
4
), pp.
1428
1459
.10.1137/110853169
74.
Miller
,
K. S.
, and
Ross
,
B.
,
1993
,
An Introduction to the Fractional Calculus and Fractional Differential Equations
,
Wiley
,
New York
.
75.
Podlubny
,
I.
,
1999
,
Fractional Differential Equations
,
,
San Diego, CA
.
76.
Atanackovic
,
T. M.
,
Pilipovic
,
S.
,
Stankovic
,
B.
, and
Zorica
,
D.
,
2014
,
Fractional Calculus With Applications in Mechanics: Vibrations and Diffusion Processes
,
Wiley
, Hoboken, NJ.
77.
Afzali
,
F.
,
Kapucu
,
O.
, and
Feeny
,
B. F.
,
2016
, “
Vibrational Analysis of Vertical-Axis Wind-Turbine Blades
,”
ASME
Paper No: DETC2016-60374.10.1115/DETC2016-60374
78.
Afzali
,
F.
,
Acar
,
G. D.
, and
Feeny
,
B. F.
,
2017
, “
Analysis of the Periodic Damping Coefficient Equation Based on Floquet Theory
,”
ASME
Paper No. DETC2017-68450
.10.1115/DETC2017-68450
79.
de Moraes
,
E. A. B.
,
Zayernouri
,
M.
, and
Meerschaert
,
M. M.
, “
An Integrated Sensitivity-Uncertainty Quantification Framework for Stochastic Phase-Field Modeling of Material Damage
,”
Int. J. Numer. Methods Eng.
(submitted).
80.
Zamani
,
V.
,
Kharazmi
,
E.
, and
Mukherjee
,
R.
,
2015
, “
Asymmetric Post-Flutter Oscillations of a Cantilever Due to a Dynamic Follower Force
,”
J. Sound Vib.
,
340
, pp.
253
266
.10.1016/j.jsv.2014.11.020
81.
Varghaei
,
P.
,
Kharazmi
,
E.
,
Suzuki
,
J. L.
, and
Zayernouri
,
M.
, “
Nonlinear Vibration of Fractional Viscoelastic Cantilever Beam: Application to Structural Health Monitoring
,”
ASME J. Vib. Acoust.
(submitted).
82.
Ari
,
M.
,
Faal
,
R. T.
, and
Zayernouri
,
M.
,
2019
, “
Vibrations Suppression of Fractionally Damped Plates Using Multiple Optimal Dynamic Vibration Absorbers
,”
Int. J. Comput. Math.
, pp.
1
24
.10.1080/00207160.2019.1594792
83.
Sapmaz
,
A.
,
Acar
,
G. D.
, and
Feeny
,
B. F.
,
2019
, “
Approximate General Responses of Tuned and Mistuned 4-Degree-of-Freedom Systems With Parametric Stiffness
,”
Topics in Modal Analysis & Testing
, Vol.
9
, Springer, Cham, Switzerland, pp.
315
324
84.
Loéve
,
M.
,
1977
,
Probability Theory
,
Springer-Verlag
,
New York
.
85.
Su
,
C. H.
, and
Lucor
,
D.
,
2006
, “
Covariance Kernel Representations of Multidimensional Second-Order Stochastic Processes
,”
J. Comput. Phys.
,
217
(
1
), pp.
82
99
.10.1016/j.jcp.2006.02.006
86.
Oksendal
,
B.
,
1998
,
Stochastic Differential Equations
,
Springer-Verlag
, New York.
87.
Foo
,
J.
,
Wan
,
X.
, and
,
G. E.
,
2008
, “
The Multi-Element Probabilistic Collocation Method (ME-PCM): Error Analysis and Applications
,”
J. Comput. Phys.
,
227
(
22
), pp.
9572
9595
.10.1016/j.jcp.2008.07.009
88.
Barthelmann
,
V.
,
Novak
,
E.
, and
Ritter
,
K.
,
2000
, “
High Dimensional Polynomial Interpolation on Sparse Grids
,”
,
12
(
4
), pp.
273
288
.10.1023/A:1018977404843
89.
Novak
,
E.
, and
Ritter
,
K.
,
1999
, “
Simple Cubature Formulas With High Polynomial Exactness
,”
Constr. Approximation
,
15
(
4
), pp.
499
522
.10.1007/s003659900119
90.
Novak
,
E.
, and
Ritter
,
K.
,
1996
, “
High Dimensional Integration of Smooth Functions Over Cubes
,”
Numer. Math.
,
75
(
1
), pp.
79
97
.10.1007/s002110050231
91.
Ervin
,
V. J.
, and
Roop
,
J. P.
,
2007
, “
Variational Solution of Fractional Advection Dispersion Equations on Bounded Domains in Rd
,”
Numer. Methods Partial Differ. Equations
,
23
(
2
), pp.
256
281
.10.1002/num.20169