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.
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.  and .
1.2 Uncertainty Framework.
Conventional approaches in parametric UQ of differential equations are based around Monte Carlo sampling (MCS) , 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 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  and multilevel sequential MCS  are also developed and recently used in Ref.  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 . Polynomial chaos, in which the output of the stochastic model is represented as a series expansion of input parameters was initially applied in Ref.  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.  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.  and  for introduction to several spectral methods for integer-order differential equations and to Refs.  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 , 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.  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
3 Forward Uncertainty Framework
3.1 Formulation of Stochastic Fractional Partial Differential Equation.
where the fractional indices and 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 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 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 , 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 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 .
where ϵ is the amplitude of process.
holds -almost surely for , subject to the homogeneous initial and boundary conditions.
3.3 Input Parametrization.
with the support constitutes a mapping of the sample space Ω onto the target space Γ. Therefore, a random vector denotes an arbitrary point in the parametric space.
holds ρ-almost surely for and , 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 for a prescribed number of realizations K.
Solving the deterministic problem Eq. (20) and obtaining the solution for each .
Computing the solution statistics, for example,
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.
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 , 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 , where 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 [39,40] with two input parameters dimension size 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.
4.2 Petrov–Galerkin Spectral Method.
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.
and 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).
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.
subject to zero initial condition, where . We let for each realization of α. In this case, by choosing the tuning parameter τ in the temporal basis function to be , 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 . 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.
with two cases and . Figure 3 shows the mean value and variance of solution for 104 sampling of MCS compared to 625 realizations in PCM.
subject to zero initial/boundary conditions, where , and the only random variables are the fractional indices α and β. We let , and choose and . For each realization of α and β, we obtain the force function h(t, x) by substituting the corresponding uext to Eq. (47). Defining , 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.
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 and , 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 and , the normalized error in computing the expectation and standard deviation of the solution are of orders and , 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.
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.
Lemma A.1. Letand. Then, the normsandare equivalent to.
We also define as the space of smooth functions with compact support in (a, b). We denote by , and as the closure of with respect to the norms , and . It is shown in Refs.  and  that these Sobolev spaces are equal and their seminorms are also equivalent to . Therefore, we can prove that and , in which β is a positive constant.
Moreover, by letting and be the space of smooth functions with compact support in and , respectively, we define and as the closure of and with respect to the norms and . Other equivalent useful seminorms associated with are also introduced in Refs.  and , as , where .
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 . For all, ifsuch that, and, then, whererepresents the standard inner product in.
Lemma A.4 . Let, a and b be arbitrary finite or infinite real numbers. Assumesuch that u(a) = 0, alsois integrable in (a, b) such that v(b) = 0. Then,.
A.3 Stability Analysis.
We show the well posedness of the deterministic problem and prove the stability of the proposed PG scheme.
Proof. The proof directly concludes from Eq. (A10) and Lemma A.6.▪
Theorem A.9 (well-posedness). For alland, and, 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 and independent of N, where .