A Collocation Approach for Solving Time-Fractional Stochastic Heat Equation Driven by an Additive Noise

: A spectral collocation approach is constructed to solve a class of time-fractional stochastic heat equations (TFSHEs) driven by Brownian motion. Stochastic differential equations with additive noise have an important role in explaining some symmetry phenomena such as symmetry breaking in molecular vibrations. Finding the exact solution of such equations is difﬁcult in many cases. Thus, a collocation method based on sixth-kind Chebyshev polynomials (SKCPs) is introduced to assess their numerical solutions. This collocation approach reduces the considered problem to a system of linear algebraic equations. The convergence and error analysis of the suggested scheme are investigated. In the end, numerical results and the order of convergence are evaluated for some numerical test problems to illustrate the efﬁciency and robustness of the presented method.


Introduction
Many models in physics, chemistry, and engineering reveal stochastic effects and are introduced as stochastic partial differential equations (SPDEs) [1,2]. Some phenomena in various fields such as population dynamics [3], motions of ions in crystals [4], optimal pricing in economics [5] and thermal noise [6] show stochastic behaviors. Fractional stochastic partial differential equation (FSPDE) is an example of these equations that have attracted more attention recently.
In recent decades, investigations have shown that fractional calculus provides some new ways for a better understanding of behaviors of real-world phenomena. Fractional-order operators give helpful tools for modeling inherited memory characteristics of real applications. Scientists proposed models for numerous phenomena in engineering, fluid mechanics, physics [7][8][9][10][11], finance [12,13], geomagnetic [14] and hydrology [15] based on fractional differential and integral equations. Non-Markovian anomalous diffusion in materials with memory, such as, viscoelastic substances is an example of these applications [16], in which the mean square displacement of particles grows faster or slower than in the case of normal diffusion.
In many applications, it is more realistic to represent the mathematical model of the problem in a non-deterministic state. In other words, some kinds of randomness and uncertainty are considered in the mathematical formulation of the problem. Hence, stochastic functional equations have arisen in many situations and numerous problems in different fields of science are modeled as fractional stochastic differential or integral equations [17][18][19]. Many theoretical investigations on the fractional stochastic differential equations have been made by researchers in the literature. Liu et al. studied some properties of fractional stochastic heat equations [20]. Ralchenko and Shevchenko [21] surveyed the existence and uniqueness of mild solution for a special type of stochastic heat equations of fractional order. Roozbahani et al. [22] proved the unique solvability of a class of SPDEs. Moghaddam et al. [23] proved the existence and uniqueness of solution for some delay stochastic differential equations of fractional order. Moreover, Mishura et al. [24] investigated mild and weak solutions for a SPDE with second order elliptic operator in divergence form. Since the exact solutions of these equations are scarcely known, researchers have examined several numerical algorithms to solve them. Finite difference schemes [25,26], finite element approaches [27][28][29], wavelets Galerkin method [30], B-spline collocation method [31,32], hat function operational matrix method [33], mean-square dissipative method [34] and operational matrix of Chebyshev wavelets [35] are a number of these schemes.
In the present work, we consider the following TFSHE where (x, t) ∈ L × I, with the boundary and initial conditions where α ∈ (0, 1), µ, ϑ and λ are real constants, I := [0, T], L := [0, l] and ∂L is the boundary of L. Also,Ḃ(t) := dB(t) dt denotes a time white noise where B(t), t ∈ I is the Brownian motion adapted to a filtration F B = {F t } t∈I in a probability space (Ω B , F B , P B ) [36]. Moreover, the source term f (x, t), ϕ(x, t) and η(x) are some stochastic processes defined on (Ω B , F B , P B ) and u(x, t) is an unknown stochastic function to be found. Moreover, the operator D α 0,t [·] denotes Caputo fractional derivative defined as: and Γ(·) represents the Gamma function. Equation (1) is a FSPDE driven by additive noise that takes into account both memory and environmental noise effects. Many physical and engineering models are built based on these types of stochastic equations. Fractional stochastic heat equations [20,[37][38][39], stochastic Burgers equation [40] and stochastic coupled fractional Ginzburg-Landau equation [41] are some examples of these applications. The problem (1)-(3) has been considered in [30], in the case α = 1. The authors have proposed a wavelet Galerkin method to find the solution to this equation. When ϑ = 0, Equation (1) reduces to an advection-dispersion equation of fractional order describing the transport of passive tracers in a porous medium in groundwater hydrology [42].
Many numerical schemes with Chebyshev polynomials basis functions are established in literature to solve various types of problems. Masjed-Jamei in [43] introduced a class of symmetric orthogonal polynomials. The six various types of Chebyshev polynomials are special cases of this basic class. To our experience, the approaches based on the SKCPs expansions result very accurate numerical estimations. Hence, we motivated to employ this kind of Chebyshev polynomials for solving TFSHEs.
The structure of this work is organized as follows. In Section 2, the basic concepts of the SKCPs theory are described. In Section 3, the collocation scheme based on the SKCPs is applied. The convergence of the numerical procedure is considered in Section 4. The accuracy of the proposed approach is analyzed in Section 5 by three numerical test problems. In the end, the main concluding remarks are presented in Section 6.

The Shifted SKCPs and Their Properties
In this section, some necessary preliminaries and relevant properties of the shifted SKCPs utilized in the next sections, are reviewed. and The explicit form of shifted SKCPs as follows: [45] J whereθ is an approximation of g(x, t), then where ξ and are two positive constants.

The SKCPs-Collocation Approach
In the following, we describe a numerical technique to solve problem (1)- (3). For this reason, we consider the numerical solution of (1) as follows where , is an unknown coefficients matrix.
is the shifted SKCPs vector as (11), then where Φ α (t) is Caputo's fractional derivative of the vectorJ(t) and is defined as where Proof. Due to the analytic form (7), we have Also, we know that [7] D α 0,t t r = for r ≥ 1. So, for j = 1, . . . , M, we get According to Equations (1) and (9) and by applying Theorem 2, we have where and from the conditions (2) and (3) and Equation (9), we have Let x 0 = 0, x N = l, and x 1 , . . . , x N−1 , are the roots of J N−1 (x). Also, suppose t j , j = 1, . . . , M, are roots ofJ M (t). By considering these collocation nodes, we define where the matrices Λ, Λ x and Λ xx are of the order (N − 1) × (N + 1) and By evaluating (16) Also, by evaluating (17) and (18) at collocation points t j and (19) at collocation points x i , we get where Using the Kronecker product, Equation (25) transforms to where and X = vec(C), T vec = vec(F ). Also, Equations (26)-(28) are equivalent tō whereĒ Thus, from Equations (29) and (30), we obtain a system of linear equations AX = B in which Solving this system leads to an estimation U N,M (x, t) for the solution of (1)-(3), in the form (9).

Convergence Analysis
In the following, we examine the convergence of the approximate solution expressed in the form (9) for the problem (1)-(3).
where R N,M (x, t) is the residual function. Now, from Equations (1) and (31), we get By using Theorem 1, we have where θ 1 is a positive constant, thus Also, from Theorem 1, we have where θ 2 and θ 3 are positive constants. Letγ = Ḃ (t) ∞ , then, from the relations (32)- (35), it can be concluded that Moreover, for x ∈ ∂L and t ∈ I, U N,M (x, t) satisfies the following equation and for x ∈ L, we have Therefore, from Equations (36)-(38), we can see that E R N,M (x, t) ∞ tends to zero, when N, M → ∞.

Applications and Results
We assess the applicability of our proposed approach to solve some stochastic heat equations of fractional order.
To simulate the Brownian motion B(t), we employ the approach described in [36]. Consider a discretization of B(t). We set t 0 = 0 and let t j , j = 1, . . . , M, are the considered collocation points, where t i < t j for i < j. Also, let B j = B(t j ) and From the definition of Brownian motion B(t) on (Ω B , ℱ B , P B ), we know that ℬ(0) = 0 with the probability 1. ℬ(τ) − ℬ(r) ∼ √ τ − rN (0, 1), for 0 ≤ r < τ ≤ T, where N (0, 1) is a normally distributed random variable with zero mean and unit variance. Also, ℬ(τ 2 ) − ℬ(τ 1 ) and ℬ(ν 2 ) − ℬ(ν 1 ) are independent for 0 ≤ τ 1 < τ 2 < ν 1 < ν 2 ≤ T. Thus, we let B 0 = t 0 with the probability 1, and where each dB j is an independent random variable of the form ∆ j N (0, 1). Throughout the section, unless stated otherwise, we assume that T = 1, l = 1 and N = M. Also, we evaluate the numerical solution u(x, t) alongP discretized paths and finally, the average of the results over these paths is considered. The L ∞ -norm error is evaluated using the following definition: where U N (ξ i , τ j ) and u(x, t), are computed by the exact and numerical solutions defined in (9) at the collocation points x = ξ i and t = τ j , respectively. The convergence order is defined by the following formula: where E N i ∞ denotes the L ∞ -norm error for N i (i = 1, 2) collocation points. The numerical computations are performed on a personal computer using a 1.70 GHz processor and the codes are written in Matlab software.

Example 1.
Consider the time-fractional stochastic equation subject to the conditions: where α ∈ (0, 1), B(t) is a Brownian motion and u(x, t) = αt 2 x 5 exp x 2 is the exact solution to the above problem. Now, we evaluate u(x, t) alongP = 80 discretized Brownian paths. To approximateḂ(t), we use the discretized scheme described at the beginning of this section. Figure 1 displays the exact and approximate solution with α = 0.9 and Figure 2 shows the exact solution and its estimations for different values of α, when N = 12. These figures confirm that the resulted numerical solutions have good compatibility with the exact solution. Table 1 displays the l ∞ -norm errors and convergence orders for α = 0.25, 0.75 and several values of N. Also, Figure 3 show the behaviour of the absolute error of u(x, t) for different values of N, when α = 0.5. Table 1 and Figure 3 confirm the accuracy of the obtained numerical approximations.
The numerical solution is evaluated alongP = 100 discretized Brownian paths. Table 2 displays the l ∞ -norm errors and order of convergence for several values of α and N. This table shows the high accuracy of the introduced scheme. Also, Figure 4 displays the exact and numerical solution of u(x, t) when α = 0.5, N = 10 and Figure 5 indicates the absolute error together with the contour plot for N = 16. It can be seen that the numerical solution is in well agreement with the exact solution.
The numerical solutions are evaluated alongP = 50 discretized Brownian paths. Figure 6 shows the numerical solution at t = 1 for different values of α, when ϑ = 0.5 and N = 10. Figure 7 displays the estimation of u(x, t) when ϑ = 0.15, 0.2, N = 8 and α = 1. The results are compared with wavelets Galerkin (WG) method [30]. This figure confirms that the present method gives more smooth solution than the numerical scheme in [30]. Also, Figures 8 and 9 indicate the approximate solutions and the contour plots for several values of ϑ when α = 0.45. The results confirm that the employed approach is very efficient.

Conclusions
According to numerous applications of FSPDEs, a new numerical scheme was introduced to solve a class of stochastic heat equations of fractional order with additive noise subject to suitable conditions. This numerical method was based on a collocation approach with the SKCPs basis functions. The convergence of the proposed method was proved. Three illustrative examples were investigated to authenticate the efficiency of the discussed approach. The obtained numerical results approved the accuracy of this method.