## 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 [1–9], 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 $\u223c1/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. [22–26]. It is also generalized and used in constructing stochastic Galerkin methods [27–30] for problems with higher-dimensional random spaces. Other nonsampling numerical methods, including but not limited to the perturbation method [31–34] 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. [37–39]; 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. [41–48] 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 [51–58] 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. [61–72].

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

*σ*are defined as (see, e.g., Refs. [74] and [75])

*σ*, $n\u22121<\sigma \u2264n,\u2009n\u2208\mathbb{N}$, defined as

*σ*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 $\sigma \u2208(0,1)$ is given by

## 3 Forward Uncertainty Framework

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

*d*and stochastic function $u(t,x;\omega ):D\xd7\Omega \u2192\mathbb{R}$, where $\omega \u2208\Omega $ denotes the random input of the system in a properly defined complete probability space $(\Omega ,F,\mathbb{P})$. We consider the following SFPDE, subject to certain homogeneous Dirichlet initial/boundary conditions and stochastic process as additional force function, given as

where the fractional indices $\alpha (\omega )\u2208(0,1)$ and $\beta j(\omega )\u2208(1,2),\u2009\u2009j=1,2,\u2026d$ are mutually independent random variables, *k _{j}* 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 $\mathbb{P}$-almost everywhere $\omega \u2208\Omega $, 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 $\alpha (\omega )\u2208(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 [76–78], material damage [79], and also the effect of damping/stiffness in structural vibrations [80–83].

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

We approximate the additional random forcing term by representing $f(t;\omega )$ 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].

*σ*-algebra of sets in Ω, and $\mathbb{P}$ is the probability measure. The random field $f(t;\omega )$ has the ensemble mean $E{f(t;\omega )}=f\xaf(t)$, finite variance $E{[f(t;\omega )\u2212f\xaf(t)]2},$ and covariance $Cf(t1,t2)=E{[f(t1;\omega )\u2212f\xaf(t1)][f(t2;\omega )\u2212f\xaf(t2)]}$. The KL expansion of $f(t;\omega )$ takes the form

*λ*are the eigenfunction and eigenvalues of the covariance kernel $Cf(t1,t2)$. We obtain the covariance kernel

_{k}*C*and its eigenvalues and eigenfunctions, following Ref. [85] and by solving a stochastic Helmholtz equation

_{f}*A*is the correlation length of $f(t;\omega )$. To ensure that the random variables $Qk(\omega )$ have zero mean and unit variance, we define them on $Qk(\omega )\u2208[\u22123,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 ($\u224890%$) 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;\omega )=1\mu \u2009\u2211k=1Mak\u2009\u2009sin(2k\pi \u2009tT)\u2009Qk(\omega )$ denote the normalized truncated expansion, assuming $f\xafM(t;\omega )=0$, where $\mu =maxt{\sigma fM}$ and $\sigma fM$ is the standard deviation of $fM(t;\omega )$. Thus, we represent the random process to be employed in Eq. (8) as*

where *ϵ* is the amplitude of process.

holds $\mathbb{P}$-almost surely for $\omega \u2208\Omega $, subject to the homogeneous initial and boundary conditions.

### 3.3 Input Parametrization.

with the support $\Gamma =\u220fi=1N\Gamma i\u2282\mathbb{R}N$ constitutes a mapping of the sample space Ω onto the target space Γ. Therefore, a random vector $\xi =(\xi 1,\xi 2,\u2026,\xi N)\u2208\Gamma $ denotes an arbitrary point in the parametric space.

holds *ρ*-almost surely for $\xi (\omega )\u2208\Gamma $ and $\u2200t,x\u2208D$, 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.

Generating a set of random variables $\xi i,\u2009i=1,2,\u2026,K$ for a prescribed number of realizations

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

Computing the solution statistics, for example, $E[u]=1M\u2211i=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 $\Theta N={\u2009\xi i\u2009}i=1J$ be a set of prescribed sampling points. By employing the Lagrange interpolation polynomials

*L*, the polynomial approximation $I$ of the stochastic solution

_{i}*u*in the random space can be expressed as

*u*is

*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,88–90].

## 4 Forward Solver

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.

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

### 4.2 Petrov–Galerkin Spectral Method.

*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\u200a\u200a\tau (x)=(1+x)\tau Pn\u22121\u2212\tau ,\tau (x)$. Thus, we employ Legendre polynomials $\varphi mj(\xi ),\u2009j=1,2,\u2026,d$, and Jacobi polyfractonomial of first kind $\psi n\tau (\eta )$ [59,60], as the spatial and temporal bases, respectively, given in their corresponding standard domain as

*U*as

_{N}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.

*V*as

_{N}and $U$ is the matrix of unknown coefficients. The matrices *S _{T}* and

*M*denote the temporal stiffness and mass matrices, respectively; and the matrices

_{T}*S*and

_{j}*M*denote the spatial stiffness and mass matrices, respectively. We obtain the entries of spatial mass matrix

_{j}*M*analytically and employ proper quadrature rules to accurately compute the entries of other matrices

_{j}*S*,

_{T}*M*, and

_{T}*S*.

_{j}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).

*Let*${emj,\lambda mj\u2009}mj=1Mj$

*be the set of general Eigen solutions of the spatial stiffness matrix S*${en\u200a\u200a\tau ,\lambda n\tau \u2009}n=1N$

_{j}with respect to the mass matrix M_{j}. Moreover, let*be the set of general eigen-solutions of the temporal mass matrix M*$U$

_{T}with respect to the stiffness matrix S_{T}. Then, the matrix of unknown coefficients*is explicitly obtained as*

## 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.

*h*(

*t*) as the external forcing term. Therefore, we obtain

subject to zero initial condition, where $u(t,\xi ):(0,T]\xd7\Lambda \u2192\mathbb{R}$. We let $uext(t)=\alpha /2\u2009\u2009t3+\alpha 2,\u2009h(t)=0Dt\alpha (\xi )uext(t)$ for each realization of *α*. In this case, by choosing the tuning parameter *τ* in the temporal basis function to be $\alpha /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 *L*^{2}-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.

*M*=

*4, as*

with two cases $h(t)=t2$ and $h(t)=sin(\pi t)$. Figure 3 shows the mean value and variance of solution for 10^{4} sampling of MCS compared to 625 realizations in PCM.

*d*=

*1 and the diffusion coefficient is*

*k*. We ignore the additional random input and consider

_{l}*h*(

*t*,

*x*) as the only external forcing term. Therefore, we obtain

subject to zero initial/boundary conditions, where $u(t,x;\xi ):(0,T]\xd7(\u22121,1)\xd7\Lambda \u2192\mathbb{R}$, and the only random variables are the fractional indices *α* and *β*. We let $uext(t,x)=t3+\tau \u2009((1+x)3+\mu \u22121/2(1+x)4+\mu )$, and choose $\tau =\alpha /2$ and $\mu =\beta /2$. For each realization of *α* and *β*, we obtain the force function *h*(*t*, *x*) by substituting the corresponding *u*^{ext} to Eq. (47). Defining $Eext[u]=E[uext]$, Fig. 4 shows the *L*^{2}-norm convergence of solution expectation as compared to the exact expectation. We observe that PCM converges accurately with only a few number of realizations.

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\u22127)$ and $O(10\u22123)$, 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.

*u*. Subsequently, we denote by $H\sigma (\Lambda )={u\u2208L2(\Lambda )\u2009|\u2009\u2203\u2009u\u0303\u2208H\sigma (\mathbb{R})\u2009\u2009s.t.\u2009\u2009u\u0303|\Lambda =u\u2009},\u2009\sigma \u22650$, the fractional Sobolev space on any finite closed interval, for example, $\Lambda =(a,b)$, with norm $||u||H\sigma (\Lambda )=infu\u0303\u2208H\mathbb{R}\sigma ,\u200au\u0303|\Lambda =u\u2009||u\u0303||H\sigma (\mathbb{R})$. We define the following useful norms as:

where the equivalence of $||\xb7||lH\sigma (\Lambda )$ and $||\xb7||rH\sigma (\Lambda )$ are shown in Refs. [54], [55], and [91].

Lemma A.1. *Let*$\sigma \u22650$*and*$\sigma \u2260n\u221212$*. Then, the norms*$||\xb7||lH\sigma (\Lambda )$*and*$||\xb7||rH\sigma (\Lambda )$*are equivalent to*$||\xb7||cH\sigma (\Lambda )$.

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

Moreover, by letting $0C\u221e(I)$ and $C0\u221e(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\u221e(I)$ and $C0\u221e(I)$ with respect to the norms $||\xb7||lHs(I)$ and $||\xb7||rHs(I)$. Other equivalent useful seminorms associated with $Hs(I)$ are also introduced in Refs. [54] and [91], as $|\xb7|lHs(I)=||0Dts(\xb7)||L2(I),\u2009|\xb7|rHs(I)=||tDTs(\xb7)||L2(I),\u2009|\xb7|Hs(I)*=|(0Dts(\xb7),tDTs(\xb7))I|1/2$, where $|\xb7|Hs(I)*\u2261|\xb7|lHs(I)1/2\u2009|\xb7|rHs(I)1/2$.

for $j=2,3,\u2026,d$.

###### Solution and Test Spaces.

##### 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*$\alpha \u2208(0,1)$*, if*$u\u2208H1([0,T])$*such that*$u(0)=0$*, and*$v\u2208H\alpha /2([0,T])$*, then*$(0Dt\u200a\u200a\alpha u,v)\Omega =(\u20090Dt\u200a\u200a\alpha /2u\u2009,\u2009tDT\u200a\u200a\alpha /2v\u2009)\Omega $*, where*$(\xb7,\xb7)\Omega $*represents the standard inner product in*$\Omega =[0,T]$.

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

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.

*The bilinear form Eq. (A8) is continuous*, that is,

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

*The following inf-sup condition holds for the bilinear form Eq. (A8)*, that is,

*where*$\Omega =I\xd7\Lambda d$*and*$supu\u2208U\u2009\u2009|a(u,v)|>0$.

Theorem A.9 (well-posedness). *For all*$0<\alpha <1$*and*$1<\beta j<2$*, and*$j=1,\u2026,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.▪

with $\gamma >0$ and independent of *N*, where $supuN\u2208UN|a(uN,vN)|>0$.

## References

^{2}N) Finite Difference Method for Fractional Diffusion Equations

*Spectral and High Order Methods for Partial Differential Equations (ICOSAHOM 2016)*(Lecture Notes in Computational Science and Engineering, Vol. 119), Springer, Cham, Switzerland, pp.

^{d}