Fractional Stochastic Partial Differential Equations: Numerical Advances and Practical Applications—A State of the Art Review

: This article aims to provide a comprehensive review of the latest advancements in numerical methods and practical implementations in the field of fractional stochastic partial differential equations (FSPDEs). This type of equation integrates fractional calculus, stochastic processes


Introduction
The study of stochastic partial differential equations (SPDEs) has emerged as a captivating area of interest, engaging a diverse community of academics and investigators from various disciplines.At its essence, SPDEs elegantly merge the unpredictability of stochastic processes with the foundational principles of partial differential equations (PDEs), providing a robust framework for the mathematical modeling of physical phenomena with randomness, alongside the spatial and temporal dynamics inherent in them.The versatility and depth of SPDEs have led to their application in a diverse range of fields, from the intricate patterns of physics to the innovative designs of engineering [1][2][3][4][5], the fluctuating dynamics of financial markets [6][7][8][9][10], and the complex processes within biological sciences [11][12][13][14][15][16].This cross-disciplinary enthusiasm highlights the crucial importance of SPDEs in deepening our understanding of the world, connecting theoretical exploration with pragmatic resolutions to complex problems.
The groundwork for the mathematical exploration of SPDEs was established by pioneers such as Pardoux [17], Bensoussan [18], and Temam [19] in the early 1970s.Their groundbreaking work opened new avenues for the rigorous analysis of SPDEs, catalyzing subsequent research into the intricate behaviors of systems influenced by randomness and uncertainty.This period marked the beginning of an in-depth investigation into the nature of random dynamical systems, setting the stage for a rich field of study that continues to evolve and expand its applications in today's complex and interconnected world.
Over the past decades, fractional calculus has proven adept at managing non-local interactions and effects, including modeling memory effects and capturing the hereditary characteristics of dynamical systems [20].These advancements have expanded the scope of mathematical modeling, allowing for more accurate representations of complex phenomena influenced by long-range dependencies and historical influences.In the realm of physics, they provide insights into the chaotic behavior of turbulent flows [21,22] and the thermal properties of materials with memory [23][24][25].The financial sector benefits from their application in the modeling of market dynamics and the valuation of financial derivatives, where long-term correlations and market volatilities present significant challenges [26,27].Similarly, in biology and epidemiology, fractional models are used to understand the spread of diseases and the movement of biological populations, taking into account both spatial dispersion and random environmental influences [28].
In recent years, fractional stochastic partial differential equations (FSPDEs) have risen to the forefront of mathematical innovation, owing to their sophisticated approach to depicting the dynamics of processes that are often oversimplified by conventional mathematical models [29].FSPDEs provide the mathematical community with a powerful toolkit to investigate, model, and address problems that involve not only randomness and extensive interactions but also memory effects.They include fractional calculus, which extends the concept of derivatives and integrals to non-integer orders to offer a more accurate and comprehensive representation of phenomena characterized by long-range dependencies, where interactions span vast temporal or spatial distances, and anomalous diffusion, which describes the spread of particles in a manner that deviates from classical Brownian motion [30].
Solving SPDEs is inherently computationally intensive due to their solutions existing in an infinite-dimensional or functional space over sets of random variables.This complexity increases with the incorporation of fractional operators that introduce memory effects into the model.In classical SPDE scenarios, leveraging the symmetry properties of the process can mitigate computational demands.Symmetry-based models are prevalent in various applications such as combustion models [31], SAR data analysis [32], autocatalytic networks [33], and stochastic models for the Navier-Stokes equations [34].Noteworthy studies include the investigation of symmetry breaking in an annular combustor model affected by deterministic electroacoustic feedback and stochastic forcing, which showed significant changes in mode behavior due to noise [31], and the evaluation of spherical symmetry in stochastic models for high-resolution PolSAR images, using statistical hypothesis testing to substantiate the model's assumptions [32].Moreover, the discussion on chiral symmetry breaking in autocatalytic networks emphasized fluctuations induced by autocatalysis [33], while explorations of symmetries and martingales in models for the Navier-Stokes equations revealed invariances under assumed symmetries [34].However, fractional operators typically challenge the preservation of periodic or symmetric solutions due to their memory characteristics, necessitating careful consideration and specific assumptions like infinite memory length or the use of left-right terminal fractional operators.To date, no published studies in the realm of FSPDEs have thoroughly investigated these phenomena, leaving the discussion on symmetric solutions in FSPDEs out of scope for this review paper.
As interest in FSPDEs grows, it is a necessity to advance the creation of advanced numerical methods and analysis techniques to tackle these equations.The proposed approaches for tackling FSPDEs can be classified into three main categories: spectral methods, finite difference methods, and meshless methods, as depicted in Figure 1.Spectral methods leverage the orthogonal properties of polynomial bases to approximate solutions of fractional operators.This approximation involves interpolating solutions across non-equidistant grid points spanning the entire domain, employing techniques such as Tau, Galerkin, or collocation approaches.Conversely, finite difference methods approximate solutions on equidistant grid points, focusing on local points where convergence is linearly related to the grid density.Distinct from both finite difference and spectral methods, meshless methods forgo the need for a predefined mesh or grid.Instead, they utilize a collection of domain-distributed points, applying interpolation techniques to these points to approximate solutions.This approach offers flexibility for complex geometries and dynamic domains.The reader can gain insight into the comparative strengths and weaknesses of the various numerical methods utilized in solving complex problems by referring to Table 1.
Nonetheless, for enhanced coherence and readability, this paper adopts a chronological organization based on publication year.While this paper does not delve into the comparative advantages and drawbacks of these methods, it offers a systematic review of their application across a spectrum of FSPDE scenarios.Through a detailed survey, our aim is to provide researchers and practitioners with a comprehensive understanding of the existing numerical methods for solving FSPDEs.This exploration encourages the pursuit of new and creative strategies to address the challenges posed by these equations.Our goal is to stimulate innovation in developing numerical methods for FSPDEs, contributing to the field's advancement.

Spatial Dimension Temporal Dimension Temporal Dimension Spatial Dimension Temporal Dimension Spatial Dimension
References: Identified as ST.

References:
Identified as SS.

References:
Identified as FT.

References:
Identified as FS.

References:
Identified as MS.

Preliminaries: Foundations of FSPDEs
This section serves as an introduction to the fundamental concepts underlying FSPDEs, which are typically formulated as: where T denotes the rule of the function u's progression over time, L defines its spatial behavior, F represents the deterministic external forces, and the product of σ and the stochastic operator Ξ quantifies the magnitude of the stochastic element involved.We shall review the basic definitions in fractional calculus, delve into stochastic processes with multiplicative noise, and explore Q-Wiener processes for modeling non-Gaussian phenomena, along with the concept of additive fractional Brownian motion.
Definition 1 (Fractional Derivative and Integral [53]).The fractional integral of the function f (t), denoted by I α 0,t f (t), is defined as where α is the fractional integral order and Γ(•) denotes the Gamma function.In addition, its fractional derivative in the sense of Caputo is defined as where n is the smallest integer greater than α.
Definition 2 (Stochastic Process [54]).A stochastic process X(t) is a collection of random variables representing a system's evolution over time.
Definition 3 (Adapted Process [54]).Let Ω, F , P, {F t } t≥0 be a filtered probability space.A stochastic process {X(t) : t ∈ [0, T]} is termed F t -adapted if, at any time t, the past and present are known, making X(t) measurable with respect to the σ-algebra F t .
Definition 4 (Multiplicative Noise [55]).A stochastic process where the noise intensity depends on the system's current state is represented mathematically as dX(t) = a(X(t), t)dW(t), where a(X(t), t) is a function of the state variable X(t) and time t, and dW(t) is a Wiener process increment.
Definition 5 (Wiener Process [56]).A Wiener process W(t) models Brownian motion, providing a continuous stochastic process with independent increments such that It exhibits long-range dependence and self-similarity, with the Hurst parameter controlling the degree of smoothness and roughness of the paths.

Chronological Progression of Numerical Methods for FSPDEs
In 2016 and 2018, respectively, Liu et al. [35,37] proposed a spectral approach utilizing Fourier spectral methods for solving one-dimensional spatial FSPDEs influenced by specific additive noises.These FSPDEs are formulated as where (−∆) α , with 1/2 < α ≤ 1, denotes the fractional Laplacian, and ∂ 2 tx W(t, x) represents the mixed second-order derivative of the Brownian sheet.In ref. [35], periodic boundary conditions are considered which enforce spatial periodicity, ensuring the continuity of function values and gradients at the boundaries.On the contrary, in [37], homogeneous Dirichlet boundary conditions are specified, which imply a zero derivative of the solution with respect to time at both boundary points.
In their research, the construction of the space fractional derivative relies on the Laplacian's eigenvalues and eigenfunctions under specific boundary conditions.The space-time noise is approximated by employing piecewise constant functions in the time direction and suitable approximations in the spatial direction.Subsequently, the resulting FSPDE approximations are tackled using Fourier spectral methods.Moreover, for the linear problem, precise error estimates in the L 2 -norm are derived, clarifying the relationships between the error bounds and the fractional powers.Furthermore, they offer a holistic approach to FSPDEs, allowing for the simultaneous consideration of all spatial points, thus enhancing the computational efficiency and providing analytical insights into stochastic processes.To tackle the nonlinear problem, they proposed employing fast Fourier transforms (FFTs), offering flexibility to adapt these algorithms for solving other FSPDEs with multiplicative noises.Additionally, they ensure a high accuracy, especially when solutions exhibit smoothness properties aligned with the FFT's spectral representation.
Applying Fourier spectral methods to solve one-dimensional FSPDEs driven by special additive noises presents notable advantages in terms of accuracy and efficiency.These methods excel in accurately approximating solutions, particularly in linear cases with smooth solutions, leveraging their capability to efficiently capture high-frequency components.In addition, using the FFT to address the nonlinear aspects of FSPDEs brings several benefits, particularly in terms of the computational efficiency and accuracy.FFT-based methods effectively manage the complexities of nonlinear FSPDEs by transforming nonlinear terms into the frequency domain, enabling easier manipulation of convolutions and products and reducing the computational complexity.Hence, the FFT proves adept at accurately capturing the dynamics of nonlinear FSPDEs, offering efficient and precise numerical solutions.Nevertheless, challenges arise when employing Fourier spectral methods for FSPDEs, such as managing certain boundary conditions, dealing with aliasing effects resulting from finite truncation of Fourier series, and addressing grid dependency issues, particularly in non-periodic domains.Additionally, nonlinear equations may introduce spectral pollution or cause convergence loss.Despite these challenges, careful consideration and management of boundary conditions, aliasing, and nonlinearities maintain Fourier spectral methods as potent tools for accurately and efficiently tackling one-dimensional FSPDEs driven by special additive noises.
In 2017, Li et al. [36] proposed a Galerkin finite element technique, classified as a spectral method, to tackle the stochastic space-time fractional wave equation driven by additive space-time white noise.That is, where (−∆) β denotes the fractional Laplacian derived from the spectral decomposition of the Dirichlet Laplacian, and W(t, x) signifies infinite-dimensional Brownian motion.This equation describes the decay of waves with a power-law effect influenced by white Gaussian noise, which is the focal point of investigation in this study.Initially, Galerkin finite element approximations serve as the primary numerical approach.This method involves discretizing the problem domain and representing the solution through finite linear combinations of basis functions, allowing for approximation across discrete elements.An essential aspect of this approach is discretizing the space-time additive noise, leading to a modeling discrepancy.Consequently, a regularized version of the stochastic space-time fractional wave equation emerges, suitable for numerical analysis.These methods offer adaptability to complex geometries and boundary conditions, provide accurate solutions with refined meshes, and achieve convergence through variational principles.Then, the regularity properties of this equation are analyzed to clarify its mathematical characteristics, stability attributes, and convergence behavior, offering crucial insights into the solution's behavior.
A nuanced understanding of both the advantages and limitations of Galerkin finite element methods is essential for their effective application in addressing the studied stochastic space-time fractional wave equation.Galerkin finite element methods, which fall under the category of spectral methods, come with computational overhead, especially for extensive or irregular domains, requiring meticulous mesh-generation procedures.Numerical stability concerns may arise, particularly in the presence of steep gradients or high-frequency oscillations, necessitating careful selection of basis functions and mesh refinement strategies.Additionally, scalability to higher dimensions may be limited due to the increased computational complexities associated with the higher dimensionality.
In 2017, Hosseini and Asgari [44] developed an efficient numerical technique for solving stochastic time fractional stiff PDEs.The focal equations are the stochastic time frac-tional Burgers equation and the stochastic time fractional Kuramoto-Sivashinsky equation, both subject to additive noise.
whereas the stochastic time fractional Kuramoto-Sivashinsky equation is expressed as: where σ denotes the noise amplitude, W(t) stands for a Wiener process on L 2 (R) (the space of real-valued square integrable functions on R), and ∂ α 0,t [•] signifies the Caputo fractional derivative of order α.
They employed polynomial chaos expansion to discretize the random variable, leading to a coupled time fractional deterministic system of PDEs.This method offers several advantages for addressing stochastic time fractional stiff PDEs.By combining polynomial chaos expansion and Fourier spectral methods, both the random variable and spatial domain are effectively and accurately discretized, respectively.Incorporating an exponential integrator method effectively handles stiffness in the resulting semi-discrete system, thereby enhancing numerical stability.Additionally, the use of a predictor-corrector method in the temporal domain further improves the solution accuracy.Numerical experiments demonstrate the method's effectiveness in solving stochastic time fractional Burgers and Kuramoto-Sivashinsky equations, highlighting its practicality and reliability.However, a notable drawback is the absence of a discussion regarding the convergence and stability of the proposed algorithm, potentially impacting its robustness and applicability to broader scenarios.Further investigation into these aspects would strengthen confidence in the method's performance and expand its utility to a wider range of problems.
In 2018, Darehmiraki [46] proposed an efficient numerical method, which falls under the finite difference and meshless categories, for solving a linear SPDE with fractional temporal orders.This FSPDE is formulated in a separable Hilbert space H and takes the following form: where ∇ represents the gradient operator, σ > 0, and α ∈ (0, 1].The positive constants v > 0 and γ > 0 quantify advection and diffusion processes, respectively.The functions f and g are known, ensuring the existence of a unique solution for the model.The term Ẇ denotes a space-time white noise with x ∈ R d , assumed to be adapted with respect to a filtered probability space Ω, F , P, {F t } t≥0 , where F is complete, and the filtration {F t , t ≥ 0} is right continuous.Ẇ(x, t) is a generalized process with covariance specified by a zero-mean Gaussian noise with the following covariance structure: where both (possibly generalized) functions ς and Λ are assumed to be nonnegative and nonnegative definite.Meshless methods based on radial basis functions (RBFs) offer a promising avenue for efficiently solving PDEs.In this approach, the time fractional derivative is approximated using a scheme of order O(τ 2−α ) through spline approximation, where 0 < α < 1, while spatial derivatives are approximated using the Kansa approach.This method demonstrates its efficiency and effectiveness through numerical examples, showcasing its prowess in solving FSPDEs.RBF-based meshless methods provide an attractive solution for numerically addressing PDEs using collocation schemes, eliminating the need for traditional mesh-based methods and their associated integrals over elements.RBFs, which depend solely on the radial distance, enable efficient function approximation and surface reconstruction from scattered data across various dimensions.Kansa's method, which directly employs RBFs, exemplifies this mesh-free approach, often utilizing the multiquadric (MQ) RBF.A notable advantage of RBFs lies in their adaptability to high-dimensional geometries without mesh generation, making them ideal for problems with intricate or irregular boundaries.Nevertheless, the accuracy and convergence of RBF-based approximations depend on factors such as RBF selection, collocation point distribution, and the interpolation strategy, necessitating meticulous analysis and consideration in their application.Despite their advantages, meshless methods based on RBFs also have limitations and considerations.One crucial factor is the dependency of the approximation accuracy on several parameters, including the choice of RBF and the distribution of collocation points.Selecting appropriate parameters can be challenging and may require careful tuning to achieve optimal results.Additionally, the computational cost of solving PDEs using RBFs can be significant, especially for large-scale problems, as the direct use of RBFs involves solving dense linear systems.Furthermore, while meshless methods eliminate the need for mesh generation, they may introduce challenges in enforcing boundary conditions and handling irregular geometries, potentially affecting the accuracy and robustness of the solution.
In 2018, Gunzburger et al. [47] applied a methodology based on the finite-difference framework.This methodology focuses on investigating the stability and convergence of numerical approximations for a time-fractional SPDE utilizing the backward Euler convolution quadrature method in time.The equation under consideration is where Ω ⊂ R d , with d ∈ {1, 2, 3}, represents a bounded region with a Lipschitz boundary ∂Ω.f (x, t) denotes a given deterministic source function, u 0 (x) is a given deterministic initial condition, and Ẇ(x, t) is a space-time white noise (the time derivative of a cylindrical Wiener process in L 2 (Ω)).The Laplacian operator is represented by ∆ : D(∆) → L 2 (Ω), defined on the domain By utilizing a discrete analogue of the inverse Laplace transform, this paper derives an integral representation of the numerical solution, which serves as the basis for proving the sharp convergence rate of the numerical approximation where Ψ n is the approximate solution at the n-th time step, τ is the time step size, α ∈ (0, 2/d), and d denotes the spatial dimension.This advancement allows for a refined error analysis, yielding convergence rates in line with the theoretical predictions.However, it is essential to note that the analysis is limited to the specific growth properties of the Laplacian operator's eigenvalues, which may restrict its applicability to more general abstract operators, such as semilinear problems and multiplicative space-time white noise.
The backward Euler convolution quadrature method, categorized as a finite difference method, serves as a numerical approach for discretizing and approximating the time evolution of differential equations, especially in time-dependent scenarios.It is commonly utilized to handle equations exhibiting stiffness or non-stiffness, including those with time delays or fractional derivatives.
In this method, time integration occurs by discretizing the convolution integral arising from representing time derivatives through convolution quadrature.Unlike conventional time-stepping techniques such as explicit or implicit schemes, which directly advance the solution in time, convolution quadrature approximates the time convolution integral using numerical quadrature methods.The backward Euler variant of convolution quadrature discretizes the integral backward in time, employing past solution values to approximate the current one.This approach proves particularly effective for stiff problems; it is adept at providing stable and accurate numerical solutions by adequately capturing inherent memory effects in time-dependent equations.
While the backward Euler convolution quadrature method offers benefits, it also presents drawbacks when applied to solve SPDEs.On the positive side, it boasts unconditional stability, ensuring resilience in handling the stiff terms and complex boundary conditions typical in SPDEs.Additionally, the method yields accurate solutions, especially for problems featuring memory effects or fractional derivatives.Its simplicity facilitates implementation, requiring minimal memory storage and thereby enhancing the efficiency of large-scale simulations.However, drawbacks include potential temporal discretization errors, particularly with large time steps or rapidly changing stochastic processes, which can compromise the solution accuracy.Moreover, computational costs may arise, especially for stiff SPDEs, and addressing non-standard boundary conditions can prove challenging.The method's limited adaptivity may also result in inefficiencies in capturing varying solution dynamics or time scales during simulations.Despite these limitations, careful parameter tuning and implementation can effectively harness the backward Euler convolution quadrature method for solving SPDEs.
In 2018, Zou [48] introduced a Galerkin finite element approach for the time-fractional stochastic heat equation driven by multiplicative noise.This equation is motivated by examining heat transfer in porous media with thermal memory and incorporating random effects.The focus is on determining the spatial and temporal regularity properties of the mild solution under specific conditions.This approach involves utilizing the conventional Galerkin finite element method in the spatial domain, falling under the category of finite difference methods, alongside the Gorenflo-Mainardi-Moretti-Paradisi (GMMP) scheme in the temporal domain, categorized as a finite difference method.They establish convergence error estimates for both the semi-discrete and fully discrete schemes.
The main equation under consideration is the time-fractional stochastic heat equation perturbed by multiplicative noise on the finite time interval [0, T]: subject to the initial condition: and the boundary conditions: where u(x, t) is a random function, D is a bounded subset of R d (where d = 1, 2, 3) with the boundary ∂Ω, and the operator ∆ represents the Laplacian acting on R d .The coefficient σ is a real-valued continuous function, and W(t) denotes the Wiener process given on a filtered probability space Ω, F , P, {F t } t≥0 .
The GMMP scheme, initially proposed in the early 2000s [59], holds significance as a numerical method tailored for solving fractional differential equations.Developed by Gorenflo, Mainardi, Moretti, and Paradisi, this scheme adeptly handles Caputo-type fractional derivatives ubiquitous in various applications of fractional calculus.Its methodology revolves around discretizing the fractional derivative term using difference approximations, often employing a convolution quadrature approach.While esteemed for its accuracy and adaptability to diverse FDEs with varying initial and boundary conditions, the GMMP scheme is not without its limitations.Its computational complexity can pose challenges, particularly for large-scale or high-dimensional systems, and parameter sensitivity necessitates careful consideration.Nonetheless, with meticulous implementation, the GMMP scheme remains a robust tool, resulting in stable numerical solutions even for stiff fractional differential equations and significantly advancing the field of fractional calculus methodologies.
In 2019, Zou [41] introduced another set of numerical methods falling under the category of finite difference methods for solving time-fractional SPDEs that model anomalous diffusion in porous media with random effects and thermal memory.The primary focus is on addressing a time-fractional stochastic Cauchy problem, which is formulated as follows: where is the Caputo fractional derivative (0 < α < 1), −△ is the Laplacian with zero Dirichlet boundary condition, I 1−α t [•] is the fractional integral operator, Ẇ(t) is white noise representing random effects, and σ is globally Lipschitz continuous.
In this investigation, the numerical methodology involves employing a conventional finite element approximation in the spatial domain, along with temporal discretization using finite difference approximation.Robust convergence error estimates are derived for both semi-discrete and fully discrete schemes within a semigroup framework.These analyses rely on error estimates for corresponding deterministic equations and the optimal regularity properties of mild solutions to time-fractional SPDEs.The proposed approach provides a flexible and accurate means of solving time-fractional SPDEs.This strategy enables efficient discretization of spatial derivatives in complex geometries and irregular domains while accurately approximating time derivatives using finite difference schemes.Nevertheless, challenges related to stability, convergence, and numerical dissipation may arise, particularly for stiff problems or when employing explicit time-stepping schemes.Achieving accurate and efficient solutions requires careful balancing of spatial and temporal discretization errors and the integration of adaptivity strategies.
In 2019, Ding and Nieto [49] studied the existence and uniqueness of mild solutions for fractional stochastic evolution equations.The focus is particularly on a system of FSPDEs represented by: The methodology employed integrates a numerical strategy combining the subspace decomposition approach with the waveform relaxation technique, known as the frequency decomposition waveform relaxation method (FDWRM).
The FDWRM, categorized under the finite difference framework, is an iterative numerical method utilized for solving systems of ordinary or partial differential equations.It divides the time domain into multiple frequency bands, treating each band separately.This decomposition simplifies the problem into distinct components, corresponding to specific frequency elements, for more efficient resolution.By iteratively solving the system across various frequency bands, the FDWRM captures the system's behavior across different time scales, particularly benefiting problems with oscillatory or multiscale characteristics.This frequency-specific approach enhances the efficiency by exploiting the problem's frequency traits, with applications in areas such as circuit simulation, signal processing, and computational electromagnetics.This paper aims to validate and demonstrate the applicability of the FDWRM through a comprehensive examination of its convergence properties.It covers definitions related to fractional calculus and the integral operator featuring the Mittag-Leffler function in the kernel.Discussions encompass the existence and uniqueness of mild solutions, the introduction of frequency subspace decomposition, and the convergence analysis of the FDWRM.By addressing fractional stochastic evolution differential equations involving almost sectorial operators, this paper introduces a methodological framework that not only establishes theoretical findings but also provides a numerical methodology for practical implementation.Leveraging the FDWRM, the proposed methodology demonstrates effectiveness in managing large-scale nonlinear systems and lends itself to parallelization.
In 2020, Mirzaee and Samadyar [50] presented a semi-discrete numerical scheme for solving time-fractional stochastic advection-diffusion equations.These equations are given by: ∂ α 0,t u(x, t) = β∂ x u(x, t) + γ∂ 2 x u(x, t) + f (x, t) + σ Ẇ(t), where 0 < α ≤ 1 denotes the fractional order of the time derivative and σ represents the amplitude of the stochastic term.This model offers a comprehensive understanding of the intricate relationship between fractional dynamics, deterministic influences, and stochastic fluctuations in advection-diffusion processes.
The methodology combines finite difference schemes and radial basis function (RBF) interpolation, categorized as finite difference and meshless methods, respectively.Initially, the time-fractional stochastic advection-diffusion equation is transformed into elliptic stochastic differential equations through a finite difference scheme.Subsequently, a meshless method utilizing RBFs approximates the resulting equation.The solution is achieved by discretizing the domain in the time direction with the finite difference method and approximating the unknown function in the space direction using generalized inverse multiquadric RBFs.Notably, direct simulation of noise terms at collocation points in each time step provides a significant advantage.The stability and convergence of the scheme are thoroughly established.
The proposed method offers several advantages for solving FPDEs.By coupling a finite difference method with RBF interpolation, it provides a semi-discrete approach that effectively captures both the temporal dynamics and spatial variations.This combination allows for a flexible and efficient numerical solution, particularly suitable for problems with complex geometries and non-uniform spatial distributions.Additionally, the method leverages the strengths of both techniques, exploiting the accuracy of finite difference methods in time discretization while benefiting from the flexibility and interpolation capabilities of RBFs in spatial approximation.However, there are also potential drawbacks to consider.The computational cost of RBF interpolation may be significant, especially for high-dimensional problems or large datasets, potentially limiting the scalability of the method.Additionally, the accuracy of the solution may depend on the choice of RBF and its parameters, requiring careful tuning and validation.Despite these challenges, the proposed method offers a promising approach for solving FSPDEs, particularly in scenarios where traditional numerical methods face limitations.
In 2020, Sweilam et al. [42] carried out a comprehensive numerical investigation into a stochastic fractional advection-diffusion model with multiplicative noise represents the Caputo fractional time derivative, and I 1−α t [•] denotes the integral of fractional order.The second-order space operator, K = −∆ + ∇, incorporates the Laplacian ∆ in one dimension and the gradient ∇.The operator K : D(K) ⊂ H → H operates within a Hilbert space H, where D(K) represents the domain of the linear operator K. Additionally, the boundary condition is set to zero, and W(t) signifies a Q-Wiener process.
Utilizing the finite difference framework, this study effectively tackles the primary equation of interest, a time fractional stochastic advection-diffusion equation, by employing the Galerkin finite element method in the spatial domain and the finite difference method in the temporal domain, with the Caputo sense employed for the fractional derivative.Notably, the inclusion of error estimations through the Galerkin finite element method enhances the rigor of the analysis, providing valuable insights into the accuracy and reliability of the numerical results.By obtaining the mild solution in terms of Mittag-Leffler's function and employing Laplace transform techniques for both physical and semi-discrete cases, the study demonstrates a deep understanding of fractional calculus principles.
The Galerkin finite element method and the standard finite element method are numerical techniques commonly used to approximate solutions to FSPDEs across diverse scientific and engineering applications.While both methods aim to provide accurate approximations, they differ primarily in their variational formulation and choice of weight functions, which play a crucial role in capturing the fractional aspects of the equations.In the standard finite element method for FSPDEs, the variational formulation involves choosing weight functions v independently of the basis functions u i , resulting in the weak form a(u, v) = L(v), where a is a bilinear form representing the differential equation, and L is a linear functional.Conversely, in the Galerkin finite element method for FSPDEs, the weight functions v are selected to be the same as the basis functions u i , yielding the weak form a(u, u i ) = L(u i ).Mathematically, this distinction reflects the choice of weight functions in the variational formulation, which can impact the accuracy and efficiency of the numerical approximation.
Furthermore, the solution process and resulting system of equations may also differ between the two methods.In the standard finite element method for FSPDEs, the solution process typically involves assembling a global stiffness matrix and solving the resulting linear system, while the Galerkin finite element method may result in a system of equations with a different structure due to the choice of weight functions.This structural difference can affect the computational cost and stability of the solution process.Despite these differences, both methods aim to offer accurate approximations to the solution of FSPDEs, with the Galerkin finite element method often exhibiting improved consistency and stability properties owing to its choice of weight functions that are tailored to the underlying fractional dynamics.Ultimately, the selection between the two methods depends on factors such as the specific problem at hand, computational resources, and the user's expertise in handling FSPDEs.
In 2020, Babaei et al. [38] explored a stochastic heat equation of fractional order with additive noise.They employed the spectral collocation method with sixth-kind Chebyshev polynomials, falling under the spectral framework, to address this problem.The stochastic heat equation considered in their study is represented by subject to the following initial and boundary conditions: where (x, t) ∈ Ω × [0, T], Ω := [0, l], and α ∈ (0, 1).Constants µ, ϑ, and λ are real, while ∂Ω represents the boundary of Ω.The source term functions f (x, t), φ(x, t), and η(x) represent the stochastic processes defined on Ω, F , P, {F t } t≥0 , with u(x, t) denoting an unknown stochastic function.Furthermore, ∂ α 0,t [•] refers to Caputo's fractional derivative.This method involves discretizing the problem domain and approximating the solution using Chebyshev polynomials, facilitating the conversion of the original differential equation into a system of linear algebraic equations.Through the solution of this system, the numerical solution to the primary problem is effectively obtained.The rigorous establishment of convergence in this spectral collocation technique ensures the accuracy and reliability of the obtained numerical results, further validated through various test problems undertaken to assess the proposed scheme's efficacy and applicability.
While the spectral collocation method with sixth-kind Chebyshev polynomials presents several advantages for solving FSPDEs, it also entails certain limitations.Despite offering a high accuracy, particularly in approximating smooth solutions, and efficiently handling non-uniform grids and boundary conditions, the method's sensitivity to the choice of basis function can lead to spectral pollution and accuracy losses if not carefully selected.Moreover, challenges such as the occurrence of the Gibbs phenomenon in problems with discontinuities or singularities and potential difficulties in handling highly oscillatory solutions or those with steep gradients underscore the need for meticulous attention to numerical stability and precision.Thus, while the spectral collocation method offers significant benefits, optimal results hinge on thoughtful consideration of the basis function and problem-specific characteristics.
In 2021, Mirzaee et al. [45] proposed a numerical solution approach for the timefractional stochastic advection-diffusion equation.The method involves employing the finite difference method and spline approximation under the categories of finite difference methods, respectively.The equation under consideration takes the general form: along with the initial and boundary conditions.Herein, β, γ, and σ denote real constants.The unknown stochastic process to be approximated is denoted as u(x, t), with the Caputo partial derivative of order α represented by ∂ α 0,t u(x, t).The proposed methodology involves transforming the equation into a system of second-order differential equations with appropriate boundary conditions.Subsequently, a suitable numerical method, such as the backward differentiation formula, is applied to solve the resulting system.Error analysis is conducted under mild conditions, disregarding error terms of O(∆t 2 ) in the system.To assess the algorithm's features, including accuracy, efficiency, and reliability, various test problems are considered.
Combining the finite difference method and spline approximation offers several advantages and disadvantages.On the positive side, the finite difference method provides a straightforward and intuitive approach for discretizing differential equations, particularly in domains with regular grids, offering simplicity in implementation and computational efficiency.Spline approximation, on the other hand, allows for flexible and accurate interpolation of functions between grid points, enabling the representation of complex geometries and smooth variations in the solution.Moreover, spline methods can handle irregularly spaced data points, providing robustness in handling non-uniform grids.However, this combined approach may suffer from limitations such as an increased computational cost, especially when dealing with high-dimensional problems or fine discretizations.Additionally, the accuracy of spline interpolation may degrade in regions with sharp gradients or discontinuities, necessitating careful consideration of mesh refinement strategies.Despite these challenges, the combination of finite difference and spline approximation methods remains a versatile and widely used approach in numerical analysis and scientific computing.
In 2021, as a continuation of their research, Mirzaee et al. [51] employed finite difference and meshless frameworks for the 2D stochastic time-fractional sine-Gordon equation where γ i (for i = 1, 2, 3, 4) are known positive constants, and the functions f (x, t) and p(x) are known.The sought-after solution is denoted by u(x, t).The function sin(u(x, t)) represents the nonlinear sine function of the dependent variable.
This equation holds significant importance in various domains such as condensed matter physics, nonlinear optics, and mathematical biology.It serves as a mathematical model for elucidating the behavior of waves or solitons in two-dimensional media under the influence of both deterministic (sine-Gordon) and stochastic (noise) effects.Understanding the dynamics governed by this equation is pivotal for a wide array of applications ranging from the modeling of complex systems to the advancement of cutting-edge technologies.
Their proposed method offers several advantages in tackling the 2D stochastic timefractional sine-Gordon equation.By combining RBFs and finite difference methods, it efficiently addresses the challenges associated with a high dimensionality and irregular domains, providing a versatile approach for modeling complex physical phenomena.The finite difference technique helps mitigate issues related to dimensionality, while RBF-based interpolation proves effective for problems defined on irregular domains.Additionally, this method provides a valuable framework for studying wave dynamics in diverse fields such as condensed matter physics and nonlinear optics, where the interplay of deterministic and stochastic effects is essential.However, the lack of specificity regarding the equation's application context and the absence of convergence analysis are notable limitations that require further investigation.
In 2022, Loghman et al. [43] examined the nonlinear vibration behavior of a microbeam with a fractional viscoelastic core and subjected to random excitation.The formulation of the nonlinear FSPDE is based on Von-Karman's nonlinear strain and the Kelvin-Voigt fractional viscoelastic model, namely, Following this, the micro-beam encounters random excitation, represented by Gaussian white noise for its amplitude, and is described by where ρ signifies the material density, A refers to the micro-beam's cross-sectional area, I represents the moment of inertia, l indicates the characteristic length of the micro-beam, and c 0 denotes any damping coefficient or external damping applied to the system.Their approach involves discretizing the stochastic equation governing the beam's motion using the finite difference method and the Galerkin finite element method within the finite difference frameworks.The white noise excitation undergoes numerical discretization, and a comparison is made between the statistical linearization method and the spectral representation method.
The proposed method showcases adaptability in accommodating the various boundary conditions of the micro-beam.Particularly, the influence of the fractional-order derivative on the micro-beam's statistical properties varies notably for 0 < α < 0.6 and 0.6 < α < 1.In this context, the methodology proposed by Logman et al. aligns with prior research highlighting the advantages of combining the finite difference method with the Galerkin finite element method.This approach has been demonstrated to offer improved versatility and accuracy in addressing different types of partial differential equations.Sweilam et al. [42] further illustrate the effectiveness of this combined approach in dealing with a stochastic fractional advection-diffusion model with multiplicative noise.By harnessing the complementary strengths of both numerical techniques, their methodology aims to provide robust and efficient solutions capable of capturing the complex dynamics inherent in stochastic systems.
In 2023, Noupelah et al. [52] studied the strong convergence of a numerical approximation method for a general time-fractional SPDE defined on a bounded domain Λ ⊂ R d (d = 1, 2, 3) with a smooth boundary.Such a model is pertinent for modeling random phenomena affecting particle transport in media with thermal memory, that is where ∂ α 0,t [•] denotes the Caputo time derivative with order 0.5 ≤ α ≤ 1, where H represents a Hurst parameter satisfying 0.5 ≤ H ≤ 1, and A stands for an unbounded (not necessarily self-adjoint) operator that is assumed to produce a semigroup S(t) = e −tA .Their paper specifies the deterministic mappings F, G, and Φ, which operate from the Hilbert space H to itself.The initial data X 0 are random, and the noise terms W(t) and B H (t) are H-valued Q-Wiener and fractional Q 1 -Brownian motions, respectively.These motions are defined in a filtered probability space Ω, F , P, {F t } t≥0 , and the covariances operators Q : H → H and Q 1 : H → H are positive and linear self-adjoint operators.
The authors establish both the existence and uniqueness of solutions and employ conventional finite elements for spatial discretization within the finite difference framework.They also utilize the generalized Euler time discretization (GETD) method for temporal discretization, which falls under the finite difference category.The manuscript rigorously demonstrates the temporal and spatial convergence of the fully discrete scheme, revealing convergence rates that depend on the regularity of the initial data, the exponent of the fractional derivative, and the Hurst parameter H.
The GETD method is a numerical technique used to solve stiff ODEs or PDEs with stiff components.It proves particularly effective for problems where stiffness arises from fast oscillations or rapidly changing dynamics.This method combines exponential time differencing with a generalized form of time-stepping to efficiently handle stiffness.Unlike traditional explicit or implicit methods, which may struggle with stiff problems due to either restrictive time step sizes or high computational costs, the GETD method strikes a balance between accuracy and efficiency.By employing a combination of exponential terms with varying coefficients, the GETD method can adaptively adjust the time integration scheme to capture the dynamics of both stiff and non-stiff components of the problem.This allows for larger time step sizes compared to traditional methods, resulting in significant computational savings while maintaining accuracy.However, combining the standard finite element method in the spatial domain with the GETD method in the temporal domain for solving time FSPDEs on a bounded domain presents challenges.While the method efficiently handles stiff temporal components, delicate stability considerations and careful parameter tuning are necessary to maintain accuracy.Implementing this combined approach involves complexities in integrating spatial and temporal discretizations, especially when accommodating non-trivial boundary conditions.Additionally, the adaptivity of the GETD method in adjusting time step sizes may lead to unexpected behavior, while inherent numerical dissipation and dispersion introduce further complexities.
In 2024, Heydari et al. [39] introduced a numerical solution method within the spectral method classification for the stochastic fractional two-dimensional Sobolev equation where u is the unknown stochastic process, Their approach employs a collocation technique that relies on discrete Chebyshev polynomials as the set of basis functions.Moreover, the collocation method employs a finite expansion of discrete Chebyshev polynomials to approximate the solution.The expansion coefficients, treated as unknowns, are determined by substituting the approximation into the stochastic fractional problem.Operational matrices related to conventional and stochastic integrals, as well as fractional and ordinary derivatives of these polynomials, are extracted and utilized.This process results in a system of linear algebraic equations, which, when solved, yields the expansion coefficients and, consequently, the solution of the original stochastic fractional problem.
The proposed method offers numerous benefits for solving the stochastic fractional two-dimensional Sobolev equation.By employing a collocation approach with discrete Chebyshev polynomials, the method establishes a systematic framework for discretizing and solving the problem, capitalizing on the efficiency and accuracy of polynomial approximations.Additionally, utilizing operational matrices simplifies the conversion of the stochastic fractional problem into a set of linear algebraic equations, streamlining the computational process.Moreover, the method's flexibility enables the inclusion of various boundary conditions and source terms, broadening its applicability to different problem scenarios.However, relying on polynomial approximations may result in accuracy and convergence challenges, especially in areas of rapid variation or near discontinuities in the solution.Furthermore, the computational cost associated with solving the resulting set of linear equations could be prohibitive for problems with large spatial domains or high-dimensional stochastic inputs.Despite these drawbacks, the proposed method offers a promising approach to tackling the stochastic fractional two-dimensional Sobolev equation, striking a balance between computational efficiency and numerical accuracy.
In 2024, Nasiri et al. [40] focused on numerically resolving a time-fractional stochastic backward parabolic equation.Moreover, incorporating stochastic elements into real-life models often gives rise to stochastic backward inverse problems, highlighting the need for robust numerical methodologies.In the field of computational methods for solving differential equations, dealing with backward problems is a big challenge, especially when we can only guess the starting information from observations made at the end [60,61].When we add randomness to real-world models, it often leads to tricky backward problems, making it important to have strong numerical methods.The authors employ a unique approach tailored to address the challenges arising from the stochastic characteristics of the problem.The central equation under examination is the time-fractional stochastic parabolic equation ∂ α 0,t u(x, t) − ∇ • (D(x)∇u(x, t)) + λ(x)u(x, t) = f (x, t) + σ(x)W(t), (x, t) ∈ Ω × (0, T), subject to the boundary and initial conditions u(x, t) = 0, u(x, 0) = ϕ(x), where the diffusion coefficient D is a known continuous function in Ω ⊂ R 2 , while λ and σ are known functions in L 2 (Ω).This equation governs the sought quasi-solution in the study.When the initial function ϕ is unknown, we encounter an inverse backward problem.
To uniquely solve this problem, an additional condition is necessary.That is u(x, T) = φ(x), where ϕ and φ are random functions.This study explores quasi-solutions for fractional stochastic backward equations, showcasing the adaptability of deterministic techniques in this scenario.The theoretical framework establishes the existence and uniqueness of these quasi-solutions, introducing a novel numerical method grounded in 2D Chebyshev wavelets and classified under the SPM-ST approach.By employing this method, their paper addresses the inherent challenges posed by ill posedness in the system of linear equations through the application of the Levenberg-Marquardt regularization technique.A comprehensive convergence analysis is conducted to underscore the reliability and efficacy of the proposed numerical algorithm in practical applications.
While this approach offers several advantages, it also presents notable challenges.On the positive side, it introduces a novel methodology for tackling the stochastic nature of the problem and obtaining quasi-solutions to fractional stochastic backward equations.The numerical technique effectively approximates the solution and is supported by the convergence analysis, boosting confidence in its reliability.However, challenges arise, particularly in managing the ill-posed system of linear equations and the potential computational complexity introduced by the regularization technique.Additionally, while the methodology demonstrates adaptability from deterministic to stochastic backward equations, practical implementation may require careful consideration of computational resources and parameter tuning.Nonetheless, this developed methodology offers a promising solution for real-world inverse problems involving time-fractional stochastic backward equations, with potential applications across various fields requiring robust and efficient numerical solutions.

Conclusions
This article provided an extensive review of the latest advancements in numerical methods and practical implementations in the field of fractional stochastic partial differential equations.By integrating fractional calculus, stochastic processes, and differential equations, FSPDEs offer a powerful framework for modeling complex dynamical systems characterized by memory and randomness.This article systematically introduced crucial foundational concepts and definitions for understanding FSPDEs, followed by a comprehensive overview of the diverse numerical methods and analytical techniques developed to address these equations.Notably, the significant expansion in numerical methods, such as spectral and finite element methods, highlights their potential for innovative applications across various disciplines.Adopting a chronological organization based on publication year enhanced the coherence and readability of this paper.While a comparative analysis of these methods was not within the scope of this review, the systematic survey offered valuable insights into their application across a spectrum of FSPDE scenarios.By providing researchers and practitioners with a comprehensive understanding of existing numerical methods, this paper encourages the pursuit of new and creative strategies to tackle the challenges posed by fractional stochastic partial differential equations.Ultimately, the aim was to stimulate innovation in developing numerical methods for FSPDEs, thereby contributing to advancements in the field.

Table 1 .
Strengths and weaknesses of numerical methods including error analysis.Can be significant, especially with large numbers of nodes Weaknesses: Limited applicability, prone to Gibbs phenomena, computational cost Truncation Error: Can be minimized by selecting appropriate interpolation schemes Discretization Error: Can be controlled through grid refinement Round-off Error: Generally low due to high-precision arithmetic Method: Finite Difference Methods Strengths: Simple implementation, efficient for irregular geometries, suitable for timedependent problems Weaknesses: Lower accuracy (linear convergence), difficulty with boundary conditions, limited scalability Truncation Error: Prone to truncation errors, especially with high-order schemes Discretization Error: Can be reduced through finer mesh resolutions Round-off Error: Moderate to high due to accumulation of errors Method: Meshless Methods Strengths: Highly flexible, no grid generation needed, suitable for complex deformations Weaknesses: Reduced accuracy, sensitivity to node distribution, increased computational cost Truncation Error: Can be controlled through adaptive refinement strategies Discretization Error: Highly dependent on node distribution and density Round-off Error: (t 2 ) − B H (t 1 ) are Gaussian distributed with mean zero and variance |t 2 − t 1 | 2H .•