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

Abstract

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 difficult 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 efficiency and robustness of the presented method.

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

$$D_{0,\mathrm{t}}^{\alpha }u(\mathrm{x},\mathrm{t})=\left(\mu +\vartheta \dot{B}\left(\mathrm{t}\right)\right)u_{\mathrm{x}\mathrm{x}}(\mathrm{x},\mathrm{t})+\lambda u_{\mathrm{x}}(\mathrm{x},\mathrm{t})+f(\mathrm{x},\mathrm{t}),$$

(1)

where $(\mathrm{x},\mathrm{t})\in \mathcal{L}\times \mathcal{I}$, with the boundary and initial conditions

$$\begin{array}{c}\hfill u(\mathrm{x},\mathrm{t})=\phi (\mathrm{x},\mathrm{t}),\mathrm{x}\in \partial \mathcal{L},\mathrm{t}\in \mathcal{I},\end{array}$$

(2)

$$\begin{array}{c}\hfill u(\mathrm{x},0)=\eta \left(\mathrm{x}\right),\mathrm{x}\in \mathcal{L},\end{array}$$

(3)

where $\alpha \in (0,1)$, $\mu$, $\vartheta$ and $\lambda$ are real constants, $\mathcal{I}:=[0,T]$, $\mathcal{L}:=[0,l]$ and $\partial \mathcal{L}$ is the boundary of $\mathcal{L}$. Also, $\dot{B}\left(\mathrm{t}\right):=\frac{\mathrm{d}B\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}$ denotes a time white noise where $B\left(\mathrm{t}\right)$, $\mathrm{t}\in \mathcal{I}$ is the Brownian motion adapted to a filtration $F_B={\left\{F_{\mathrm{t}}\right\}}_{\mathrm{t}\in \mathcal{I}}$ in a probability space $({\mathsf{\Omega }}_B,F_B,{\mathbb{P}}_B)$ [36]. Moreover, the source term $f(\mathrm{x},\mathrm{t})$, $\phi (\mathrm{x},\mathrm{t})$ and $\eta \left(\mathrm{x}\right)$ are some stochastic processes defined on $({\mathsf{\Omega }}_B,F_B,{\mathbb{P}}_B)$ and $u(\mathrm{x},\mathrm{t})$ is an unknown stochastic function to be found. Moreover, the operator $D_{0,\mathrm{t}}^{\alpha }[·]$ denotes Caputo fractional derivative defined as:

$$\begin{array}{c}\hfill D_{0,\mathrm{t}}^{\alpha }u(x,\mathrm{t})=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\mathrm{t}}\frac{1}{{(\mathrm{t}-\xi )}^{\alpha }}\frac{\partial u}{\partial \xi }(\mathrm{x},\xi )\mathrm{d}\xi ,\alpha \in (0,1),\end{array}$$

(4)

and $\mathsf{\Gamma }(·)$ 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 $\alpha =1$. The authors have proposed a wavelet Galerkin method to find the solution to this equation. When $\vartheta =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. Recently, a few authors applied the SKCPs to solve some types of differential equations [44,45,46].

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.

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

Definition 1.

The shifted SKCPs on $[0,l]$ are defined by

$${\mathcal{J}}_m\left(x\right)={\widehat{\mathcal{J}}}_m((2/l)x-1),m=0,1,2,\dots ,$$

where ([43])

$${\widehat{\mathcal{J}}}_m\left(x\right)=\prod _{i=0}^{\lfloor \frac{m}{2}\rfloor -1}\frac{2i+{(-1)}^{m+1}+4}{-5-(2i+2\left[\frac{m}{2}\right]+{(-1)}^{m+1})}{\mathcal{E}}_m\left(x\right),$$

(5)

and

$${\mathcal{E}}_m\left(x\right)=\sum _{\tau =0}^{\lfloor \frac{m}{2}\rfloor }\prod _{\kappa =0}^{\lfloor \frac{m}{2}\rfloor -(\tau +1)}\left({\displaystyle \frac{{(-1)}^m-2(\kappa +\lfloor \frac{m}{2}\rfloor )-5}{{(-1)}^{m+1}+2(\kappa +2)}}\right)\frac{\left(\lfloor \frac{m}{2}\rfloor \right)!}{\tau !(\tau -\lfloor \frac{m}{2}\rfloor )!}{x}^{m-2\tau }.$$

(6)

The explicit form of shifted SKCPs as follows: [45]

$${\mathcal{J}}_m\left(x\right)=\sum _{r=0}^m{\overline{\theta }}_{r,m}{(x/l)}^r,$$

(7)

where

$${\overline{\theta }}_{r,m}=\left\{\begin{array}{ll}\frac{2^{2r-m}}{(2r+1)!}\sum _{i=\lfloor \frac{r+1}{2}\rfloor }^{\frac{m}{2}}{\displaystyle \frac{{(-1)}^{\frac{m}{2}+i+r}(2i+r+1)!}{(2i-r)!}},& meven,\\ \frac{2^{2r-m+1}}{(m+1)(2r+1)!}\sum _{i=\lfloor \frac{r}{2}\rfloor }^{\frac{m-1}{2}}{\displaystyle \frac{{(-1)}^{\frac{m+1}{2}+i+r}(i+1)(2i+r+2)!}{(2i-r+1)!}},& modd.\end{array}\right.$$

Theorem 1.

([46]) Suppose ${\mathrm{L}}_{\mathtt{W}}^2\left(\tilde{\mathsf{\Lambda }}\right)$ is the square integrable function space according to the weight $\mathtt{W}(x,t)={(2x-1)}^2{(2t-1)}^2\sqrt{x-x^2}\sqrt{t-t^2}$. Let $g(x,t)\in L_{\mathtt{W}}^2\left(\tilde{\mathsf{\Lambda }}\right)$ is considered with ${∥\frac{{\partial }^{6}g(x,t)}{\partial x^3\partial t^3}∥}_2\le \varsigma$ for some constant $\varsigma >0$, satisfies the expansion $g(x,t)={\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^{\infty }}{c}_{i,j}{\mathcal{J}}_i\left(x\right){\mathcal{J}}_j\left(t\right)$. If

$${\mathtt{G}}_{N,M}(x,t)=\sum _{i=0}^N\sum _{j=0}^M{c}_{i,j}{\mathcal{J}}_i\left(x\right){\mathcal{J}}_j\left(t\right),$$

(8)

is an approximation of $g(x,t)$, then

$$|g(x,t)-{\mathtt{G}}_{N,M}(x,t)|<\frac{\varsigma }{2^{N+M}},$$

$$|\frac{\partial g}{\partial x}(x,t)-\frac{\partial {\mathtt{G}}_{N,M}}{\partial x}(x,t)|<\xi \frac{N}{2^{N+M-2}},$$

$$|\frac{{\partial }^{2}g}{\partial x^2}(x,t)-\frac{{\partial }^2{\mathtt{G}}_{N,M}}{\partial x^2}(x,t)|<\varrho \frac{N^3}{2^{N+M-8}},$$

where ξ and ϱ are two positive constants.

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

$$u(\mathrm{x},\mathrm{t})\simeq {\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})=\sum _{i=0}^N\sum _{j=0}^M{\delta }_{i,j}{\mathcal{J}}_i\left(\mathrm{x}\right){\overline{\mathcal{J}}}_j\left(\mathrm{t}\right)=\mathbf{J}{\left(\mathrm{x}\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(\mathrm{t}\right),$$

(9)

where

$$\begin{array}{c}\hfill \mathbf{J}\left(\mathrm{x}\right)={[{\mathcal{J}}_0\left(\mathrm{x}\right),\dots ,{\mathcal{J}}_i\left(\mathrm{x}\right),\dots ,{\mathcal{J}}_N\left(\mathrm{x}\right)]}^T,\end{array}$$

(10)

$$\begin{array}{c}\hfill \overline{\mathbf{J}}\left(\mathrm{t}\right)={[{\overline{\mathcal{J}}}_0\left(\mathrm{t}\right),\dots ,{\overline{\mathcal{J}}}_j\left(\mathrm{t}\right),\dots ,{\overline{\mathcal{J}}}_M\left(\mathrm{t}\right)]}^T,\end{array}$$

(11)

in which ${\mathcal{J}}_i\left(\mathrm{x}\right)={\widehat{\mathcal{J}}}_i((2/l)\mathrm{x}-1)$ on the interval $\mathcal{L}$ and ${\overline{\mathcal{J}}}_j\left(\mathrm{t}\right)={\widehat{\mathcal{J}}}_i((2/T)\mathrm{t}-1)$ on the interval $\mathcal{I}$. Moreover

$$\mathbf{C}={\left(\begin{array}{ccc}{\delta }_{0,0}& \cdots & {\delta }_{0,M}\\ ⋮& \ddots & ⋮\\ {\delta }_{N,0}& \cdots & {\delta }_{N,M}\end{array}\right)}_{(N+1)\times (M+1)},$$

is an unknown coefficients matrix.

Theorem 2.

Let $\overline{\mathbf{J}}\left(t\right)$ is the shifted SKCPs vector as (11), then

$$D_{0,t}^{\alpha }\overline{\mathbf{J}}\left(t\right)={\mathsf{\Phi }}^{\alpha }\left(t\right),$$

(12)

where ${\mathsf{\Phi }}^{\alpha }\left(t\right)$ is Caputo’s fractional derivative of the vector $\overline{\mathbf{J}}\left(t\right)$ and is defined as

$${\mathsf{\Phi }}^{\alpha }\left(t\right)={\left[0,\sum _{r=1}^1{\psi }_{r,1}^{\alpha }\left(t\right),\dots ,\sum _{r=1}^j{\psi }_{r,j}^{\alpha }\left(t\right),\dots ,\sum _{r=1}^M{\psi }_{r,M}^{\alpha }\left(t\right)\right]}^T,$$

(13)

where

$${\psi }_{r,j}^{\alpha }\left(t\right)=\frac{\mathsf{\Gamma }(r+1)}{T^r\mathsf{\Gamma }(r+1-\alpha )}{\overline{\theta }}_{r,j}{t}^{r-\alpha }.$$

Proof. 

Due to the analytic form (7), we have

$$D_{0,\mathrm{t}}^{\alpha }{\overline{\mathcal{J}}}_0\left(\mathrm{t}\right)={\overline{\theta }}_{0,0}{D}_{0,\mathrm{t}}^{\alpha }\left(1\right)=0,$$

(14)

Also, we know that [7]

$$D_{0,\mathrm{t}}^{\alpha }{\mathrm{t}}^r=\frac{\mathsf{\Gamma }(r+1)}{\mathsf{\Gamma }(r+1-\alpha )}{\mathrm{t}}^{r-\alpha }.$$

(15)

for $r\ge 1$. So, for $j=1,\dots ,M$, we get

$$\begin{array}{c}\hfill D_{0,\mathrm{t}}^{\alpha }{\overline{\mathcal{J}}}_j\left(\mathrm{t}\right)=\sum _{r=0}^j{\overline{\theta }}_{r,j}{D}_{0,\mathrm{t}}^{\alpha }{(\mathrm{t}/T)}^r=\sum _{r=1}^j{\psi }_{r,j}^{\alpha }\left(\mathrm{t}\right)\end{array}$$

in which ${\psi }_{r,j}^{\alpha }\left(\mathrm{t}\right)=\frac{\mathsf{\Gamma }(r+1)}{T^r\mathsf{\Gamma }(r+1-\alpha )}{\overline{\theta }}_{r,j}{\mathrm{t}}^{r-\alpha }$. □

According to Equations (1) and (9) and by applying Theorem 2, we have

$$\mathbf{J}{\left(\mathrm{x}\right)}^T\mathbf{C}{\mathsf{\Phi }}^{\alpha }\left(\mathrm{t}\right)=\left(\mu +\vartheta \dot{B}\left(\mathrm{t}\right)\right){\mathbf{J}}_{\mathrm{x}\mathrm{x}}{\left(\mathrm{x}\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(\mathrm{t}\right)+\lambda {\mathbf{J}}_{\mathrm{x}}{\left(\mathrm{x}\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(\mathrm{t}\right)+f(\mathrm{x},\mathrm{t}),$$

(16)

where

$${\mathbf{J}}_{\mathrm{x}}\left(\mathrm{x}\right)={[{\mathcal{J}}_0^{{}^{\prime }}\left(\mathrm{x}\right),\cdots ,{\mathcal{J}}_i^{{}^{\prime }}\left(\mathrm{x}\right),\cdots ,{\mathcal{J}}_N^{{}^{\prime }}\left(\mathrm{x}\right)]}^T,$$

$${\mathbf{J}}_{\mathrm{x}\mathrm{x}}\left(\mathrm{x}\right)={[{\mathcal{J}}_0^{{}^{″}}\left(\mathrm{x}\right),\cdots ,{\mathcal{J}}_i^{{}^{″}}\left(\mathrm{x}\right),\cdots ,{\mathcal{J}}_N^{{}^{″}}\left(\mathrm{x}\right)]}^T,$$

and from the conditions (2) and (3) and Equation (9), we have

$$\begin{array}{c}\hfill \mathbf{J}{\left(0\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(\mathrm{t}\right)=\phi (0,\mathrm{t}),\end{array}$$

(17)

$$\begin{array}{c}\hfill \mathbf{J}{\left(l\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(\mathrm{t}\right)=\phi (l,\mathrm{t}),\end{array}$$

(18)

$$\begin{array}{c}\hfill \mathbf{J}{\left(\mathrm{x}\right)}^T\mathbf{C}\overline{\mathbf{J}}\left(0\right)=\eta \left(\mathrm{x}\right).\end{array}$$

(19)

Let ${\mathtt{x}}_0=0,{\mathtt{x}}_N=l,$ and ${\mathtt{x}}_1,\dots ,{\mathtt{x}}_{N-1},$ are the roots of ${\mathcal{J}}_{N-1}\left(\mathrm{x}\right)$. Also, suppose ${\mathtt{t}}_j,j=1,\dots ,M,$ are roots of ${\overline{\mathcal{J}}}_M\left(\mathrm{t}\right)$. By considering these collocation nodes, we define

$$\begin{array}{ccc}\hfill \mathsf{\Lambda }& =& {\left[\mathbf{J}\left({\mathtt{x}}_1\right),\dots ,\mathbf{J}\left({\mathtt{x}}_i\right),\dots ,\mathbf{J}\left({\mathtt{x}}_{N-1}\right)\right]}^T,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\mathsf{\Lambda }}_{\mathrm{x}}& =& {\left[{\mathbf{J}}_{\mathrm{x}}\left({\mathtt{x}}_1\right),\dots ,{\mathbf{J}}_{\mathrm{x}}\left({\mathtt{x}}_i\right),\dots ,{\mathbf{J}}_{\mathrm{x}}\left({\mathtt{x}}_{N-1}\right)\right]}^T,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\mathsf{\Lambda }}_{\mathrm{x}\mathrm{x}}& =& {\left[{\mathbf{J}}_{\mathrm{x}\mathrm{x}}\left({\mathtt{x}}_1\right),\dots ,{\mathbf{J}}_{\mathrm{x}\mathrm{x}}\left({\mathtt{x}}_i\right),\dots ,{\mathbf{J}}_{\mathrm{x}\mathrm{x}}\left({\mathtt{x}}_{N-1}\right)\right]}^T,\hfill \end{array}$$

(22)

where the matrices $\mathsf{\Lambda }$, ${\mathsf{\Lambda }}_{\mathrm{x}}$ and ${\mathsf{\Lambda }}_{\mathrm{x}\mathrm{x}}$ are of the order $(N-1)\times (N+1)$ and

$$\begin{array}{ccc}\hfill \mathsf{\Psi }& =& {\left[\overline{\mathbf{J}}\left({\mathtt{t}}_1\right),\dots ,\overline{\mathbf{J}}\left({\mathtt{t}}_j\right),\dots ,\overline{\mathbf{J}}\left({\mathtt{t}}_M\right)\right]}_{(M+1)\times M},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}^{\alpha }& =& {\left[{\mathsf{\Phi }}^{\alpha }\left({\mathtt{t}}_1\right),\dots ,{\mathsf{\Phi }}^{\alpha }\left({\mathtt{t}}_j\right),\dots ,{\mathsf{\Phi }}^{\alpha }\left({\mathtt{t}}_M\right)\right]}_{(M+1)\times M}.\hfill \end{array}$$

(24)

By evaluating (16) at $(N-1)\times M$ collocation points $({\mathtt{x}}_i,{\mathtt{t}}_j)$ for $i=1,\dots ,N-1$ and $j=1,\dots ,M$, we have

$$\mathsf{\Lambda }\mathbf{C}{\mathsf{\Psi }}^{\alpha }={\mathsf{\Lambda }}_{xx}\mathbf{C}\mathsf{\Psi }\mathcal{B}+\lambda {\mathsf{\Lambda }}_x\mathbf{C}\mathsf{\Psi }+\mathcal{F},$$

(25)

where

$$\mathcal{B}=diag\left(\mu +\vartheta {\mathrm{b}}_1,\cdots ,\mu +\vartheta {\mathrm{b}}_j,\cdots ,\mu +\vartheta {\mathrm{b}}_M\right),$$

in which ${\mathrm{b}}_j=B\left({\mathtt{t}}_j\right)-B\left({\mathtt{t}}_{j-1}\right)$, ${\mathtt{t}}_0=0$ and

$$\mathcal{F}={\left[f_{i,j}\right]}_{(N-1)\times M},f_{i,j}=f({\mathtt{x}}_i,{\mathtt{t}}_j),i=1,\dots ,N-1,j=1,\dots ,M.$$

Also, by evaluating (17) and (18) at collocation points ${\mathtt{t}}_j$ and (19) at collocation points ${\mathtt{x}}_i$, we get

$$\begin{array}{c}\hfill \mathbf{J}{\left(0\right)}^T\mathbf{C}\mathsf{\Psi }={\mathrm{Y}}_0,\end{array}$$

(26)

$$\begin{array}{c}\hfill \mathbf{J}{\left(l\right)}^T\mathbf{C}\mathsf{\Psi }={\mathrm{Y}}_l,\end{array}$$

(27)

$$\begin{array}{c}\hfill \overline{\mathsf{\Lambda }}\mathbf{C}\overline{\mathbf{J}}\left(0\right)=\overline{\mathrm{Y}},\end{array}$$

(28)

where

$${\mathrm{Y}}_0={[\phi (0,{\mathtt{t}}_1),\dots ,\phi (0,{\mathtt{t}}_j),\dots ,\phi (0,{\mathtt{t}}_M)]}^T,{\mathrm{Y}}_l={[\phi (l,{\mathtt{t}}_1),\dots ,\phi (l,{\mathtt{t}}_j),\dots ,\phi (l,{\mathtt{t}}_M)]}^T,$$

$$\overline{\mathsf{\Lambda }}={\left[\mathbf{J}\left({\mathtt{x}}_0\right),\dots ,\mathbf{J}\left({\mathtt{x}}_i\right),\dots ,\mathbf{J}\left({\mathtt{x}}_N\right)\right]}^T,\overline{\mathrm{Y}}={[\eta \left({\mathtt{x}}_0\right),\dots ,\eta \left({\mathtt{x}}_i\right),\dots ,\eta \left({\mathtt{x}}_N\right)]}^T.$$

Using the Kronecker product, Equation (25) transforms to

$$\begin{array}{c}\hfill \mathcal{A}\mathcal{X}={\mathcal{T}}_{\mathrm{vec}},\end{array}$$

(29)

where

$$\mathcal{A}={{\mathsf{\Psi }}^{\alpha }}^T\otimes \mathsf{\Lambda }-{\left(\mathsf{\Psi }\mathcal{B}\right)}^T\otimes {\mathsf{\Lambda }}_{xx}-\lambda {\mathsf{\Psi }}^T\otimes {\mathsf{\Lambda }}_x,$$

and $\mathcal{X}=\mathrm{vec}\left(\mathbf{C}\right),{\mathcal{T}}_{\mathrm{vec}}=\mathrm{vec}\left(\mathcal{F}\right).$ Also, Equations (26)–(28) are equivalent to

$$\begin{array}{c}\hfill \overline{\mathtt{E}}\mathcal{X}=\overline{\mathrm{Y}},{\mathtt{E}}_0\mathcal{X}={\mathrm{Y}}_0,{\mathtt{E}}_l\mathcal{X}={\mathrm{Y}}_l,\end{array}$$

(30)

where

$$\overline{\mathtt{E}}=\overline{\mathbf{J}}{\left(0\right)}^T\otimes \overline{\mathsf{\Lambda }},{\mathtt{E}}_0={\mathsf{\Psi }}^T\otimes \mathbf{J}{\left(0\right)}^T,{\mathtt{E}}_l={\mathsf{\Psi }}^T\otimes \mathbf{J}{\left(l\right)}^T.$$

Thus, from Equations (29) and (30), we obtain a system of linear equations $\mathbf{A}\mathcal{X}=\mathbf{B}$ in which

$$\mathbf{A}={\left[{\mathcal{A}}^T,{\overline{\mathtt{E}}}^T,{\mathtt{E}}_0^T,{\mathtt{E}}_l^T\right]}^T,\mathbf{B}={\left[{\mathcal{T}}_{\mathrm{vec}}^T,{\overline{\mathrm{Y}}}^T,{\mathrm{Y}}_0^T,{\mathrm{Y}}_l^T\right]}^T.$$

Solving this system leads to an estimation ${\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})$ for the solution of (1)–(3), in the form (9).

4. Convergence Analysis

In the following, we examine the convergence of the approximate solution expressed in the form (9) for the problem (1)–(3).

Theorem 3.

Let ${\mathtt{U}}_{N,M}(x,t)$ is the approximate solution obtained by the procedure presented in Section 3 and $u(x,t)$ is the exact solution of (1)–(3). Consider the residual error ${\mathcal{R}}_{N,M}(x,t)$ of this numerical solution. Then, $\mathbb{E}\parallel {\mathcal{R}}_{N,M}{(x,t)\parallel }_{\infty }$ tends to zero, when $N\to \infty$ and $M\to \infty$.

Proof. 

Suppose ${\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})$, for $(\mathrm{x},\mathrm{t})\in \mathcal{L}\times \mathcal{I}$, satisfies the equation

$$\begin{array}{c}\hfill D_{0,\mathrm{t}}^{\alpha }{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})=\left(\mu +\vartheta \dot{B}\left(\mathrm{t}\right)\right)\frac{{\partial }^2{\mathtt{U}}_{N,M}}{\partial {\mathrm{x}}^2}(x,t)\\ \hfill +\lambda \frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{x}}(\mathrm{x},\mathrm{t})+f(\mathrm{x},\mathrm{t})+{\mathcal{R}}_{N,M}(\mathrm{x},\mathrm{t}),\end{array}$$

(31)

where ${\mathcal{R}}_{N,M}(\mathrm{x},\mathrm{t})$ is the residual function. Now, from Equations (1) and (31), we get

$$\begin{array}{cc}\hfill \mathbb{E}\parallel {\mathcal{R}}_{N,M}{(\mathrm{x},\mathrm{t})\parallel }_{\infty }\le & \mathbb{E}\parallel D_{0,\mathrm{t}}^{\alpha }\left(u(\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})\right){\parallel }_{\infty }\\ & \hfill +\mathbb{E}\left(\parallel \mu +\vartheta \dot{B}{\left(\mathrm{t}\right)\parallel }_{\infty }\parallel u_{\mathrm{x}\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{{\partial }^2{\mathtt{U}}_{N,M}}{\partial {\mathrm{x}}^2}(\mathrm{x},\mathrm{t}){\parallel }_{\infty }\right)\\ & \hfill +\left|\lambda \right|\mathbb{E}\parallel u_{\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{x}}(\mathrm{x},\mathrm{t}){\parallel }_{\infty }.\end{array}$$

(32)

By using Theorem 1, we have

$$\begin{array}{c}\hfill \parallel u_{\mathrm{t}}(\mathrm{x},\mathrm{t})-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{t}}(\mathrm{x},\mathrm{t}){\parallel }_{\infty }=\underset{(\mathrm{x},\mathrm{t})\in \mathcal{L}\times \mathcal{I}}{sup}|u_{\mathrm{t}}(\mathrm{x},\mathrm{t})-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{t}}(\mathrm{x},\mathrm{t})|\\ \hfill <\frac{{\theta }_{1}M}{2^{N+M-2}},\end{array}$$

where ${\theta }_1$ is a positive constant, thus

$$\begin{array}{c}\hfill \mathbb{E}\parallel D_{0,\mathrm{t}}^{\alpha }\left(u(\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})\right){\parallel }_{\infty }\le \mathbb{E}\left({\int }_0^{\mathrm{t}}\frac{{\parallel (\mathrm{t}-\tau )}^{-\alpha }{\parallel }_{\infty }}{\mathsf{\Gamma }(1-\alpha )}\parallel u_{\tau }(\mathrm{x},\tau )-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \tau }(\mathrm{x},\tau ){\parallel }_{\infty }\mathrm{d}\tau \right)\\ \hfill <\frac{{\theta }_{1}M}{\mathsf{\Gamma }(1-\alpha )2^{N+M-2}}\mathbb{E}\left({\int }_0^{\mathrm{t}}{\parallel {(\mathrm{t}-\tau )}^{-\alpha }\parallel }_{\infty }\mathrm{d}\tau \right).\end{array}$$

Since $0<\tau <\mathrm{t}\le T$, hence, we get

$$\begin{array}{c}\hfill \mathbb{E}\parallel D_{0,\mathrm{t}}^{\alpha }\left(u(\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})\right){\parallel }_{\infty }<\frac{{\theta }_1{T}^{1-\alpha }M}{\mathsf{\Gamma }(1-\alpha )2^{N+M-2}}.\end{array}$$

(33)

Also, from Theorem 1, we have

$$\begin{array}{c}\hfill \parallel u_{\mathrm{x}\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{{\partial }^2{\mathtt{U}}_{N,M}}{\partial {\mathrm{x}}^2}(\mathrm{x},\mathrm{t}){\parallel }_{\infty }=\underset{(\mathrm{x},\mathrm{t})\in \mathcal{L}\times \mathcal{I}}{sup}|u_{\mathrm{x}\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{{\partial }^2{\mathtt{U}}_{N,M}}{\partial {\mathrm{x}}^2}(\mathrm{x},\mathrm{t})|\\ \hfill <\frac{{\theta }_2{N}^3}{2^{N+M-8}},\end{array}$$

(34)

$$\begin{array}{c}\hfill \parallel u_{\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{x}}(\mathrm{x},\mathrm{t}){\parallel }_{\infty }=\underset{(\mathrm{x},\mathrm{t})\in \mathcal{L}\times \mathcal{I}}{sup}|u_{\mathrm{x}}(\mathrm{x},\mathrm{t})-\frac{\partial {\mathtt{U}}_{N,M}}{\partial \mathrm{x}}(\mathrm{x},\mathrm{t})|\\ \hfill <\frac{{\theta }_{3}N}{2^{N+M-2}},\end{array}$$

(35)

where ${\theta }_2$ and ${\theta }_3$ are positive constants. Let $\overline{\gamma }={\parallel \dot{B}\left(\mathrm{t}\right)\parallel }_{\infty }$, then, from the relations (32)–(35), it can be concluded that

$$\begin{array}{cc}\hfill \mathbb{E}\parallel {\mathcal{R}}_{N,M}{(\mathrm{x},\mathrm{t})\parallel }_{\infty }& <\frac{{\theta }_1{T}^{1-\alpha }M}{\mathsf{\Gamma }(1-\alpha )2^{N+M-2}}+\left(\right|\mu |+\overline{\gamma }\left|\vartheta \right|)\frac{{\theta }_2{N}^3}{2^{N+M-8}}+\left|\lambda \right|\frac{{\theta }_{3}N}{2^{N+M-2}}\hfill \\ & <\frac{{\theta }_1{T}^{1-\alpha }M}{\mathsf{\Gamma }(1-\alpha )2^{N+M-8}}+\left(\right|\mu |+\overline{\gamma }\left|\vartheta \right|)\frac{{\theta }_2{N}^3}{2^{N+M-8}}+\left|\lambda \right|\frac{{\theta }_3{N}^3}{2^{N+M-8}}\hfill \\ & <\widehat{\theta }\frac{M+2N^3}{2^{N+M-8}},\hfill \end{array}$$

(36)

where

$$\widehat{\theta }=max\{\frac{{\theta }_1{T}^{1-\alpha }}{\mathsf{\Gamma }(1-\alpha )},\left(\left|\mu \right|+\overline{\gamma }\left|\vartheta \right|\right){\theta }_2,\left|\lambda \right|{\theta }_3\}.$$

Moreover, for $\mathrm{x}\in \partial \mathcal{L}$ and $\mathrm{t}\in \mathcal{I}$, ${\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})$ satisfies the following equation

$$\begin{array}{cc}\hfill \mathbb{E}\parallel {\mathcal{R}}_{N,M}{(\mathrm{x},\mathrm{t})\parallel }_{\infty }& =\mathbb{E}{\parallel \phi (\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})\parallel }_{\infty }\hfill \\ & =\mathbb{E}\parallel u(\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}{(\mathrm{x},\mathrm{t})\parallel }_{\infty }\hfill \\ & =\underset{(\mathrm{x},\mathrm{t})\in \partial \mathcal{L}\times \mathcal{I}}{sup}|u(\mathrm{x},\mathrm{t})-{\mathtt{U}}_{N,M}(\mathrm{x},\mathrm{t})|<\frac{{\theta }_4}{2^{N+M}},\hfill \end{array}$$

(37)

and for $\mathrm{x}\in \mathcal{L}$, we have

$$\begin{array}{c}\hfill \mathbb{E}\parallel {\mathcal{R}}_{N,M}{(\mathrm{x},0)\parallel }_{\infty }=\mathbb{E}{\parallel \eta \left(\mathrm{x}\right)-{\mathtt{U}}_{N,M}(\mathrm{x},0)\parallel }_{\infty }\\ \hfill =\mathbb{E}\parallel u(\mathrm{x},0)-{\mathtt{U}}_{N,M}{(\mathrm{x},0)\parallel }_{\infty }\\ \hfill =\underset{x\in \mathcal{L}}{sup}|u(\mathrm{x},0)-{\mathtt{U}}_{N,M}(\mathrm{x},0)|<\frac{{\theta }_5}{2^{N+M}}.\end{array}$$

(38)

Therefore, from Equations (36)–(38), we can see that $\mathbb{E}\parallel {\mathcal{R}}_{N,M}{(\mathrm{x},\mathrm{t})\parallel }_{\infty }$ tends to zero, when N, $M\to \infty$. □

5. 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\left(\mathrm{t}\right)$, we employ the approach described in [36]. Consider a discretization of $B\left(\mathrm{t}\right)$. We set ${\mathtt{t}}_0=0$ and let ${\mathtt{t}}_j,j=1,\dots ,M$, are the considered collocation points, where ${\mathtt{t}}_i<{\mathtt{t}}_j$ for $i<j$. Also, let $B_j=B\left({\mathtt{t}}_j\right)$ and

$${\Delta }_j={\mathtt{t}}_j-{\mathtt{t}}_{j-1},j=1,\dots ,M.$$

(39)

From the definition of Brownian motion $B\left(t\right)$ on $({\mathsf{\Omega }}_B,{\mathcal{F}}_B,{\mathbb{P}}_B)$, we know that $\mathcal{B}\left(0\right)=0$ with the probability 1. $\mathcal{B}\left(\tau \right)-\mathcal{B}\left(r\right)\sim \sqrt{\tau -r}\mathcal{N}(0,1)$, for $0\le r<\tau \le T$, where $\mathcal{N}(0,1)$ is a normally distributed random variable with zero mean and unit variance. Also, $\mathcal{B}\left({\tau }_2\right)-\mathcal{B}\left({\tau }_1\right)$ and $\mathcal{B}\left({\nu }_2\right)-\mathcal{B}\left({\nu }_1\right)$ are independent for $0\le {\tau }_1<{\tau }_2<{\nu }_1<{\nu }_2\le T$. Thus, we let $B_0={\mathtt{t}}_0$ with the probability 1, and

$$B_j=B_{j-1}+\mathrm{d}{B}_j,j=1,\dots ,M,$$

(40)

where each $\mathrm{d}{B}_j$ is an independent random variable of the form $\sqrt{{\Delta }_j}\mathcal{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(\mathrm{x},\mathrm{t})$ along $\overline{\mathrm{P}}$ discretized paths and finally, the average of the results over these paths is considered.

The $L_{\infty }$-norm error is evaluated using the following definition:

$$\parallel {\mathsf{E}}_N{\parallel }_{\infty }=\underset{1\le i,j\le N}{max}|u({\xi }_i,{\tau }_j)-{\mathtt{U}}_N({\xi }_i,{\tau }_j)|,$$

(41)

where ${\mathtt{U}}_N({\xi }_i,{\tau }_j)$ and $u(\mathrm{x},\mathrm{t})$, are computed by the exact and numerical solutions defined in (9) at the collocation points $\mathrm{x}={\xi }_i$ and $\mathrm{t}={\tau }_j$, respectively. The convergence order is defined by the following formula:

$$\mathtt{CO}={log}_{\frac{N_1}{N_2}}\frac{\parallel {\mathsf{E}}_{N_1}{\parallel }_{\infty }}{\parallel {\mathsf{E}}_{N_2}{\parallel }_{\infty }},$$

(42)

where $\parallel {\mathsf{E}}_{N_i}{\parallel }_{\infty }$ denotes the $L_{\infty }$-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

$$\begin{array}{c}\hfill D_{0,t}^{\alpha }u(x,t)=\dot{B}\left(t\right)u_{xx}(x,t)+u_x(x,t)+f(x,t),\end{array}$$

subject to the conditions:

$$\begin{array}{c}\hfill u(x,0)=0,\\ \hfill u(0,t)=0,u(1,t)=\alpha exp\left(1\right)t^2,\end{array}$$

where $\alpha \in (0,1)$, $B\left(t\right)$ is a Brownian motion and

$$\begin{array}{c}\hfill f(x,t)=\frac{\alpha \mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }(3-\alpha )}{t}^{2-\alpha }{x}^{5}exp\left(x^2\right)-\alpha t^2{x}^{4}exp\left(x^2\right)(5+2x^2)\\ \hfill -2\alpha \dot{B}\left(t\right)t^2{x}^{3}exp\left(x^2\right)(10+11x^2+2x^4).\end{array}$$

$u(x,t)=\alpha t^2{x}^{5}exp\left(x^2\right)$ is the exact solution to the above problem.

Now, we evaluate $u(\mathrm{x},\mathrm{t})$ along $\overline{\mathrm{P}}=80$ discretized Brownian paths. To approximate $\dot{B}\left(t\right)$, we use the discretized scheme described at the beginning of this section. Figure 1 displays the exact and approximate solution with $\alpha =0.9$ and Figure 2 shows the exact solution and its estimations for different values of $\alpha$, when $N=12$. These figures confirm that the resulted numerical solutions have good compatibility with the exact solution.

Table 1. Example 1: The $l_{\infty }$-norm errors and convergence orders.

N	$\mathit{\alpha }=0.25$		$\mathit{\alpha }=0.75$
	$\parallel {\mathsf{E}}_{\mathit{N}}{\parallel }_{{}_{\infty }}$	CO		$\parallel {\mathsf{E}}_{\mathit{N}}{\parallel }_{{}_{\infty }}$	CO
6	$1.6754\times {10}^{-2}$	$-$		$1.1450\times {10}^{-2}$	$-$
9	$1.7591\times {10}^{-4}$	$11.2375$		$1.8956\times {10}^{-3}$	$4.4354$
12	$2.5055\times {10}^{-7}$	$22.7823$		$9.5343\times {10}^{-6}$	$18.3967$
15	$1.0902\times {10}^{-9}$	$24.3666$		$8.4121\times {10}^{-8}$	$21.1988$

displays the $l_{\infty }$-norm errors and convergence orders for $\alpha =0.25,0.75$ and several values of N. Also, Figure 3 show the behaviour of the absolute error of $u(\mathrm{x},\mathrm{t})$ for different values of N, when $\alpha =0.5$. Table 1 and Figure 3 confirm the accuracy of the obtained numerical approximations.

Example 2.

Suppose the time-fractional stochastic equation

$$\begin{array}{c}\hfill D_{0,t}^{\alpha }u(x,t)=\left(\pi +\dot{B}\left(t\right)\right)u_{xx}(x,t)-2u_x(x,t)+f(x,t),\end{array}$$

where $\alpha \in (0,1)$, $B\left(t\right)$ is a Brownian motion, and

$$\begin{array}{c}f(x,t)=\frac{3}{\mathsf{\Gamma }(2-\alpha )}{t}^{1-\alpha }\left(x^2-\frac{2tx}{\alpha -2}+\frac{2t^2}{(\alpha -2)(\alpha -3)}\right)sin\left(\pi x\right)\\ -(\pi +\dot{B}\left(t\right))(x+t)[6sin\left(\pi x\right)+6\pi (x+t)cos\left(\pi x\right)\\ -{\pi }^2{(x+t)}^{2}sin\left(\pi x\right)]+2{(x+t)}^2\left(\pi (x+t)cos\left(\pi x\right)+3sin\left(\pi x\right)\right).\end{array}$$

With these assumptions, the exact solution is $u(x,t)={(x+t)}^{3}sin\left(\pi x\right)$.

The numerical solution is evaluated along $\overline{\mathrm{P}}=100$ discretized Brownian paths.

Table 2. Example 2: The $l_{\infty }$-norm errors and convergence order.

N	$\mathit{\alpha }=0.25$		$\mathit{\alpha }=0.75$
	$\parallel {\mathsf{E}}_{\mathit{N}}{\parallel }_{{}_{\infty }}$	CO		$\parallel {\mathsf{E}}_{\mathit{N}}{\parallel }_{{}_{\infty }}$	CO
6	$7.6688\times {10}^{-3}$	$-$		$6.2094\times {10}^{-3}$	$-$
9	$2.8269\times {10}^{-4}$	$8.1401$		$1.3207\times {10}^{-4}$	$9.4964$
12	$1.7723\times {10}^{-7}$	$25.6347$		$1.6552\times {10}^{-8}$	$23.2270$
15	$9.9296\times {10}^{-11}$	$33.5528$		$7.3368\times {10}^{-11}$	$34.6026$

displays the $l_{\infty }$-norm errors and order of convergence for several values of $\alpha$ and N. This table shows the high accuracy of the introduced scheme. Also, Figure 4 displays the exact and numerical solution of $u(\mathrm{x},\mathrm{t})$ when $\alpha =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.

Example 3.

Let

$$\begin{array}{c}\hfill D_{0,t}^{\alpha }u(x,t)=\left(\frac{1}{{\pi }^2}+\vartheta \dot{B}\left(t\right)\right)u_{xx}(x,t),\end{array}$$

subject to:

$$\begin{array}{c}\hfill u(0,t)=u(1,t)=0,\\ \hfill u(x,0)=sin\left(\pi x\right),\end{array}$$

where $\alpha \in (0,1)$ and $B\left(t\right)$ is a Brownian motion.

The numerical solutions are evaluated along $\overline{\mathrm{P}}=50$ discretized Brownian paths. Figure 6 shows the numerical solution at $t=1$ for different values of $\alpha$, when $\vartheta =0.5$ and $N=10$. Figure 7 displays the estimation of $u(\mathrm{x},\mathrm{t})$ when $\vartheta =0.15,0.2$, $N=8$ and $\alpha =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, Figure 8 and Figure 9 indicate the approximate solutions and the contour plots for several values of $\vartheta$ when $\alpha =0.45$. The results confirm that the employed approach is very efficient.

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