Analyzing the continuity of the mild solution in finite element analysis of semilinear stochastic subdiffusion problems

1.   Introduction

Fractional calculus has gained increasing attention due to its potential applications in various fields of science and engineering ^[19]. However, the study of fractional calculus has been predominantly focused on deterministic equations, using either deterministic or probabilistic methods. This approach limits the modeling of real-world phenomena where the propagation speed can be finite, as heat flow can be disrupted by material response. In contrast, the classical heat equation assumes infinite speed of heat flow. Recent studies have shown that materials with thermal memory can exhibit finite heat flow speed ^[2]. This is due to the convolution term in the definition of fractional derivatives and integrals, which implies that the nearer past has a stronger influence on the present. Additionally, if the internal energy of the material is affected by past random effects, it can be modeled as fractionally integrated additive noise, represented as ${_{0}{\rm I}_{t}^{\gamma}}\dot{W}(t)$ using the classical Wiener process. Here, ${_{0}{\rm I}_{t}^{\gamma}} u$ (or ${_{RL}{\rm D}_{0, t}^{-\gamma}} u$) denotes the Riemann-Liouville fractional integral of order $\gamma$ of the function $u$ defined by

In this paper, we consider the time fractional semilinear stochastic partial differential equation driven by fractionally integrated additive noise,

where $0 < \alpha < 1$, $0 \leq \gamma \leq 1$, $A = -\Delta$ is a self-adjoint, positive definite, not necessarily bounded operator on the Hilbert space with domain $D(A) = H^2(D)\cap H_0^1(D)$, where $D \subset \mathbb{R}^{d}, d = 1, 2, 3$ denotes a bounded convex polygonal domain. Here, $\dot{W}(t) = \frac{{\rm d} W(t)}{{\rm d} t}$ denotes the white noise, and the time fractional derivative $_{C}{\rm D}_{0, t}^{\alpha} u$ with order $\alpha \in (0, 1)$ is defined as follows ^[12],

where $\Gamma(\cdot)$ denotes the Gamma function.

In recent years, the numerical solution of fractional partial differential equations has become a major focus of research. This is because analytical expressions for such equations are generally difficult to obtain, prompting researchers to explore numerical methods instead. The presence of singular convolution kernels in fractional operators further complicates the solving process. To tackle this challenge, a plethora of remarkable mathematical techniques and innovative approaches have emerged in recent years. These techniques not only play a crucial role in theoretical investigations but also demonstrate great potential in practical applications. For example, a novel technique employing double reduction order and a newly constructed nonlinear compact difference operator has been developed to simulate nonlocal problems on graded meshes ^[13]. In a separate study, researchers ^[15] have devised a conservative, positivity-preserving, nonlinear finite volume scheme suitable for multi-term nonlocal Nagumo-type equations using distorted meshes. Additionally, ^[16] proposed a positivity-preserving finite volume scheme tailored for subdiffusion equations on nonconforming quadrilateral distorted meshes with hanging nodes. Furthermore, ^[20] addresses the numerical solution of the three-dimensional nonlocal evolution equation with a weakly singular kernel.

Many researchers are also dedicated to studying techniques for solving stochastic partial differential equations. Interested readers are encouraged to explore ^[5,6,14,17] for further works in this area. In prior research, our focus revolved around weak convergence analysis of the L1 scheme for a stochastic subdiffusion problem, as well as leveraging the spectral method for strong approximation of stochastic semilinear subdiffusion and superdiffusion equations driven by fractionally integrated additive noise ^[8,9]. In this article, we delve deeper into the analysis of non-smooth data in this model. Notably, our model problem (1.1), in contrast to the formulation studied in ^[6,14], exhibits greater generality and necessitates the involvement of three distinct Mittag-Leffler solution operators (namely $E(t)$, $\overline{E}(t)$, and $\widetilde{E}(t)$, as detailed in Section 2) due to the presence of time fractional derivative and fractionally integrated additive noise.

Upon applying the Riemann-Liouville derivative operator $_{RL}{\rm D}_{0, t}^{1-\alpha} : = ({_{0}{\rm I}_{t}^{\alpha}})'$ on both sides of (1.1), it is then formally equivalent to a semilinear fractional Volterra type evolution equation

so the existence and uniqueness of a mild solution $u$ can be proved according to the literature methods ^[1] analogously, even only under some assumptions, via a standard Banach fixed point argument.

The main contributions of the paper are the following:

(i) The paper introduces finite element analysis for semilinear stochastic subdiffusion problems driven by fractionally integrated additive noise. It explores the smoothness of the solution and employs complex integration techniques to approximate the error of the solution operator under non-smooth data.

(ii) The paper establishes the continuity conditions of the mild solution, providing insights into the behavior and regularity of the solution when dealing with non-smooth data. We accurately demonstrate the impact of the fractional parameters $\alpha$ and $\gamma$ on the convergence rates through numerical examples, offering valuable insights into the sensitivity and dependence of the solution on these parameters.

The remaining sections of this paper are structured as follows. In the upcoming section, we lay out the framework for the stochastic partial differential equation (SPDE) and establish significant smoothing properties using Mittag-Leffler functions. Subsequently, we derive essential error estimates for deterministic subdiffusion under non-smooth data, followed by obtaining optimal convergence estimates for the finite element method in the presence of non-smooth data.

2.   Preliminaries

It is known that if $A = -\Delta$ with homogenous Dirichlet boundary conditions, one has $A\varphi_k = \lambda_k \varphi_k, \; k \in \mathbb{N},$ where $0 < \lambda_1 \leq \lambda_2 \leq \cdots \leq \cdots \leq \lambda_k \leq \cdots$, $\mathop {\lim }\limits_{k \to \infty } {\lambda _k} = \infty$, and the eigenvectors $\{\varphi _k\}_{k = 1}^{\infty}$ form an orthonormal basis for $H$. Let $H = L_2(D)$ be a separable Hilbert space with inner product $(\cdot, \cdot)$ and norm $\|\cdot\|$. Let $(\Omega, \mathscr{F}, \mathbb P, \{\mathscr{F}_{t}\}_{t\geq 0})$ be a filtered probability space, with Bochner spaces $L_p(\Omega; H) = L_p((\Omega, \mathscr{F}, \mathbb P); H)$. Let $\mathbb{E}$ denote the expectation (with respect to $\mathbb{P}$). We recall an abstract framework to describe the noise $W(t)$ in the model (1.1) more precisely. A Wiener process $W(t)$ with a covariance operator $Q$ may be characterized by the Fourier type series as follows:

where $Q$ is a bounded, linear, self-adjoint, positive definite operator on $H$, with the pairs of eigenvalue and eigenfunction $\{(\mu_k, \varphi_k)\}_{k = 1}^{\infty}$. The $\{ {\beta _k}(t)\} _{k = 1}^\infty$ is a sequence of independently and identically distributed standard Brownian motions.

For any $\nu\in \mathbb{R}$, we introduce the space ${{\dot H}^\nu }(D) = D({A^{\frac{\nu }{2}}})$ with norm $|v|_\nu ^2 = \|{A^{\frac{\nu }{2}}}v\|{^2} = \sum\limits_{k = 1}^\infty {\lambda _k^\nu } {(v, {\varphi _k})^2}$, where $\{\varphi _k\}_{k = 1}^{\infty}$ is the orthonormal basis in $H$.

Let $\mathcal{L} = \mathcal{L}(H)$ denote the space of all bounded linear operators on $H$ and let $\mathcal{L}_{2}^{0} = HS(Q^{\frac{1}{2}}(H), H)$ be the space of Hilbert-Schmidt operators from $Q^\frac{1}{2}(H)$ to $H$, i.e.,

furnished with the norm $\big\|T\big\|_{\mathcal{L}_2^0 }^2 = {\sum\limits_{k = 1}^\infty {\big\|T {Q^{\frac{1}{2}}}{\varphi _k}\big\|} ^2},$ thus $\|T\|_{\mathcal{L}_2^0 } = \|T{Q^{\frac{1}{2}}}\|_{HS} < \infty$, for $T \in \mathcal{L}_2^0$.

We also need to recall the Burkholder-Davis-Gundy inequality ^[11], for $p \geq 2$,

for strongly measurable functions $\phi: [0, T] \rightarrow \mathcal{L}_2^0$. It's important to emphasize that when $p = 2$, it is the Ito isometry.

By using time fractional Duhamel's principle and Laplace transform, we can obtain the mild solution of (1.1) as follows:

where

The two parameter function of the Mittag-Leffler type plays a very important role in the fractional calculus ^[5]. We recall the following important properties of the Mittag-Leffler function essential in our analysis.

Lemma 2.1. Let $0 < \alpha < 2$ and $\beta \in \mathbb{R}$ be arbitrary and $\frac{{\pi \alpha }}{2} < \mu < \min (\pi, \alpha \pi)$, then there exists a constant $C = C(\alpha, \beta, \mu)$ such that

where the notation ${\mathbb{Z}^ - }$ denotes the set of nonpositive integers, i.e., ${\mathbb{Z}^ - } = \{0, \; -1, \; -2, \ldots \}$.

Throughout we always make the following standing assumption on the fractional orders $\alpha$ and $\gamma$, which is sufficient to ensure the well-posedness of Eq (1.1) (see ^[2,5]).

Assumption 2.1. Let $0 < \alpha < 1$, $0 \leq \gamma \leq 1$ with $\alpha + \gamma > \frac{1}{2}$.

On the Assumption 2.1, for $\beta \in (0, \kappa]$, the regularity of noise can be characterized as follows (see Lemma A.1 of ^[5]):

where (with $\varepsilon > 0$ small)

Now we state the smoothing property of the operators $E(t)$, $\overline{E}(t)$, and $\widetilde{E}(t)$.

Lemma 2.2. There exists $C$ such that for $t > 0$, we have

Proof. We just prove (2.13). The other conclusions can be obtained by using the similar method. By Lemma 2.1, we have

This completes the proof of the lemma. □

Remark 2.1. Overall, the corresponding conclusions are the smoothing properties of the heat semigroup, when $\alpha \rightarrow 1$ and $\gamma \rightarrow 0$, see ^[17].

In the qualitative theory of nonlinear PDEs, the subsequent Gronwall type inequalities play a very important role in error estimate.

Lemma 2.3. ^[3] Let $T > 0$, $N \in \mathbb{N}$, $k = \frac{T}{N}$, and $t_n = nk$ for $0\leq n\leq N$. If $\zeta_1, ..., \zeta_N \geq 0$ satisfy for some $M_0, \; M_1 \geq 0$, and $\mu, \; \nu > 0$, the inequality

then there exists a constant $M_2 = M_2(\mu, \nu, M_1, T)$ such that ${\zeta _n} \le {M_0}{M_2}(1 + t_n^{ - 1 + \mu }), \; 1 \le n \le N$.

To mimic Assumption 2.14 as stated in book ^[7], let $\mathcal{P}_T$ be the $\sigma$-field of predictable stochastic processes, and $\mathcal{B}(S)$ be the Borel $\sigma$-field of $S$. Given the available options, we shall make some reasonable assumptions about the nonlinear parts.

Assumption 2.2. The mapping $f:[0, T] \times \Omega \times H \rightarrow \dot{H}^{-1}$, $(t, \omega, h) \rightarrow f(t, \omega, h)$ is $\mathcal{P}_T \times \mathcal{B}(H)/\mathcal{B}(\dot{H}^{-1})$ measurable. When $\delta \in (0, \frac{1}{2})$, there exists a constant $C$ that satisfies the following expression:

for all $h \in H$, $0 \leq t_1 < t_2 \leq T$, $\omega \in \Omega$.

3.   Non-smooth data analysis for stochastic problem

In this section, we formulate the Galerkin finite element methods for spatial discretization in combination with time discretization based on an exponential Euler type method for approximation of (1.1). Let $\mathcal{T}_h$ be a regular shaped quasi-uniform triangulation of the domain $D$, and let $S_h \subset H_0^1(D)$ be the space of continuous piecewise linear functions on the triangulation $\mathcal{T}_h$. We define the $L_{2}$-projection ${P_h}:H \to {S_h}$ by

and the Ritz projection ${R_h}:H_0^1 \to {S_h}$ by

where $a(u, \chi) = (\nabla u, \nabla \chi)$ is the associated bilinear form.

Note that by interpreting the righthand side of (3.1) as a duality pairing between ${{\dot H}^1}(D)$ and ${{\dot H}^{-1}}(D)$, one may extend $P_h$ to be a bounded operator from ${{\dot H}^{-1}}(D)$ to $S_{h}$. It is well-known that the operators $P_h$ and $R_h$ have the following approximation properties ^[18].

Lemma 3.1. The operators $P_h$ and $R_h$ satisfy

The semi-discrete Galerkin FEM scheme for (1.1) is to find $u_h(t)\in S_h$ such that

where the discrete Laplacian $A_h$ is defined by

Naturally, we present the discrete analogues of operators $E_h(t)$, $\overline{E}_h(t)$, and $\widetilde{E}_h(t)$ as follows

Now, let $0 < {t_0} < {t_1} < \cdots < {t_M} = T$ be a uniform partion of time interval $[0, T]$, with time step $\Delta t = t_{m+1}-t_{m}, \; m = 0, 1, \cdots, M-1$. Hence, the fully discrete approximation of (1.1) is given by

with initial value $U_h^0 = P_hu_0$.

Let us introduce and prove some Lemmas that will play an important role later on.

Lemma 3.2. ^[10] Let $0 \le \omega \le \mu \le 2$. For $\alpha \in (0, 1)$, there exists a constant $C$ such that

We now turn to the non-smooth data error estimates of the approximations to $\overline{E}(t) g, \, g \in H$ in the semi-discrete case.

Lemma 3.3. Let $0\leq s\leq 1$ and $0\leq r \leq 2$ with $r+s \leq 2$. For $g\in H$, there holds

Proof. In the case $s = 0$, by the inverse Laplace transform, for any given $g \in H$, we have

where $\Gamma_{\theta, \delta}^{\tau} = \{ z \in \Gamma_{\theta, \delta}: |\Im z | \leq \frac{\pi}{\tau} \}$, and $\Gamma_{\theta, \delta} = \{ z \in\mathbb{C}: z = r e^{\pm i \theta}, r \geq \delta \} \cup \{ z \in \mathbb{C}: z = \delta e^{i \phi}, |\phi| \leq \theta \}, \; \frac{\pi}{2} < \theta < \pi$, $\frac{\pi}{\tau} > \delta$.

Let us first show (3.7). For any fixed $g \in H$, we have, by (3.8) and (3.9),

By ^[4], p. 820, we get $\big \| \big ((z^{\alpha} +A)^{-1} - (z^{\alpha} +A_{h})^{-1}P_{h} \big) g \big \| \leq C h^2 \| g\|, \; \forall z \in \Gamma_{\theta, \delta}$. Hence, we have

For $I$, with $z = \delta e^{i \phi},$ we have $\delta = t^{-1}$, where $t^{-1} < \frac{\pi}{\tau}$ for sufficiently small $\tau$,

For $II$, with $z = r e^{\pm i \theta}$, we have $r \geq \delta, \, \delta = t^{-1}$,

where $c, \; C$ is some suitable positive constant.

Meanwhile, by (2.12) and the triangle inequality,

Similarly, for $s = 1$, there holds

Now the desired assertion follows by interpolation, which completes the proof. □

The next lemma gives an error estimate on $\widetilde{E}(t)$. More details can be found in Lemma 4.4 of ^[5].

Lemma 3.4. Let $0\leq s\leq 1$ and $0\leq r \leq 2$ with $r+s \leq 2$. For $g\in H$, there holds

Based on the previous discussion, we are ready to prove the error estimates for the fully discrete approximation.

Theorem 3.1. For $0\leq r \leq 2$ and ${r\alpha+2\gamma > 1}$, by the Assumptions 2.1 and 2.2, with $\nu \in [0, \beta)$, $\beta \in (0, \kappa]$, $\alpha + \gamma \in (\frac{1}{2}, 1)$ and $\alpha (2 - \kappa + \beta - \nu) + 2\gamma -1 \in (0, 1)$, then there holds

Proof. Due to (1.1), it is formally then equivalent to a nonlinear fractional Volterra type evolution equation. We can obtain $\mathop {\sup}\limits_{t \in [0, T]} \|u(t)\|{_{{L_p}(\Omega; {{\dot H}^\nu })}} \le C$. Subtracting (3.6) from (2.3) and by taking norms, one obtains

Here, we note that $I_1$, $I_2$, and $I_3$ correspond to the spatial finite element discretization error, while $I_4$ corresponds to the temporal error.

The estimate of $I_1$ is a consequence of Lemma 3.2. For $\omega = 0$ and $\mu = 2-r$, we get

For $I_2$, by using (3.7) and $r > 0$, we have

In the case $r = 0$, similarly, we split $I_2$ into two terms,

For $I_{21}$, noting $\| \overline E (\sigma)\| \leq C \sigma^{\alpha-1}$ and $\| \overline E_h (\sigma)\| \leq C \sigma^{\alpha-1}$, then

For $I_{22}$, by (3.7) with $r = 0$, we derive

where $l_h = \ln\frac{t_m}{h^{2-r}}$.

Now, we estimate $I_3$. Using (3.10), by Itô's formula, for trace class and ${r\alpha+2\gamma > 1}$, we obtain

For $I_4$, by Lemma 2.2, the idea is similar to ^[9] with different coefficient condition $\alpha (2 - \kappa + \beta - \nu) + 2\gamma -1 > 0$, so we get $I_4 \leq \Delta t^{\alpha\nu}$. Finally, combining the Gronwall's inequality of Lemma 2.3, we complete the proof of Theorem 3.1. □

Under the assumptions of smoothness of the solution operator and the nonlinear term, we can obtain the following conclusions.

Theorem 3.2. Let $p \in [2, \infty)$ be given such that Assumptions 2.13–2.17 hold, then the unique mild solution $u$ to (2.3) is continuous with respect to the norm $\|\cdot\|_{L_{p}(\Omega; \dot{H}^s)}$.

Proof. According to the definition of continuity, our goal is to show that

with either $t_1$ or $t_2$ fixed.

Therefore, with the expression of a mild solution $u$, we divide $\left\|u\left(t_2\right)-u\left(t_1\right)\right\|_{L_p\left(\Omega; \dot{H}^s\right)}$ into five parts, and only need to show that each of the parts approaches zero when $t_2 \rightarrow t_1$. By utilizing the triangle inequality, we obtain

For the expression $M_1$, we apply (2.11) from Lemma 2.2. One yields

and the validity of $\mathop {\lim }\limits_{{t_2} - {t_1}} {M_1} = 0$ here is obvious when $\nu \geq 0$.

For $M_2$, we insert a node $t_3$ and split it into three parts.

For $M_{21}$, we apply (2.12) from Lemma 2.2, and Assumption 2.2. We obtain

Due to $\delta \in (0, \frac{1}{2})$, therefore $\frac{1-s}{2}\alpha + \delta > 0$, and thus when $t_2 \rightarrow t_1$, it is valid.

For $M_{22}$,

which is formally identical to $M_{21}$. When $t_2 \rightarrow t_3 \in [t_1, t_2]$ holds, the conclusion is valid.

For $M_{23}$, with the help of (2.12), one gets

To ensure the continuity property holds, it is sufficient to only have the integration in the last inequality be well-defined and satisfy the conditions of $\frac{1+\gamma-s}{2}\alpha > 0$.

For $M_{3}$, we prove that the following property holds first. One obtains

and then we have

Since ${\frac{{1 + r - s}}{2}\alpha } > 0$ holds here, the conclusion is valid.

For $M_4$, by applying (2.14) and the Ito formula, we can obtain

In order for the conclusion to hold, it is necessary for the conditions of ${(1 - \frac{s}{2})\alpha + \gamma - \frac{1}{2}} > 0$ to be satisfied here.

For $M_5$, according to the requirement, we first prove that the following property holds.

then we have

The last term of the inequality only needs to satisfy the condition of ${(1 - \frac{s}{2})\alpha + \gamma - \frac{1}{2}} > 0$. Thus, the theorem is proved. □

4.   Numerical implementation

We conducted several numerical experiments in this section to validate our previous theoretical findings. To illustrate the theoretical results obtained in Theorem 3.1 for $\alpha \in (0, 1)$, we provide two numerical examples. Specifically, we set $T = 0.1$, $D = (0, 1)$, and $F(u) = \sin(u)$.

To explain the computer implementation of the full-discrete method (3.6), we make the assumption that the covariance operator $Q$ possesses the same eigenfunctions as $A$, such that $Qv = \sum\limits_{k = 1}^\infty {{\mu _k}(v, {\varphi _k}){\varphi _k}}$. Furthermore, we assume that $W(t)$ has the following Fourier series expansion: $W(t) = \sum_{k = 1}^{\infty} \mu_{k}^{1/2} \varphi_{k} \beta_{k}(t)$. Thus, the semi-discrete solution ${U_{m, k, h}^N}$ satisfies

where ${u_{0, k}} = ({u_0}, {\varphi _k})$, ${F_k}(\cdot) = (F(\cdot), {\varphi _k})$, and $\beta_k(\sigma)$ for $k = 1, \dots, N$ are mutually independent standard Brownian motions.

In our experiments, we investigated the dependence of the error estimates in Theorem 3.1 on the time step size $\Delta t$. To approximate the exact solution, we used the fully discrete solution with a small time step size of $\Delta t = 2^{-10}$ and a spatial step size of $\Delta x = 2^{-10}$. The theoretical rate of convergence is $\alpha\nu$ for $\nu \in (0, 1)$, which approaches $\alpha$ as $\nu \rightarrow 1$.

To conduct the experiments, we fixed a small spatial step size $\Delta x$ and considered a sequence of moderate time step sizes $\Delta t_i = 2^{-i}, \; i = 2, \dots, 5$. We performed $M = 100$ simulations for each time step size $\Delta t_i$. In each simulation $\omega_j$, $j = 1, 2, \dots, M$, we generated $N_h$ independent Brownian motions $\beta_k(t)$, $k = 1, 2, \dots, N_h$. The initial data was set to $u_0 = 1$, and we defined the nonlinear operator $F(u) = \sin(u)$.

In Figure 1, we investigated the effects of different parameters $\alpha$ and $\gamma$. We observed that the numerical results are consistent with the theoretical results stated in Theorem 3.1. The numerical values vary slightly due to the limited range of $\alpha$ and $\gamma$. As the mesh size is refined, the numerical experiments support the theoretical findings in Theorem 3.1.

In the following analysis, we focus on the spatial convergence. To this end, we consider a fixed number of space steps $M = 200$ and final time $T = 1$, and we obtain the reference solution at $N = 480$. Next, we compute the numerical results for different combinations of fractional orders $\alpha$ and $\gamma$, with trace class noise (with $m = 2$). Table 1 illustrates the results, which show that an $O(h^2)$ convergence order is observed for all combinations. We observed that the initial convergence rate changes slowly, but as the mesh is refined, the results are in excellent agreement with the theoretical predictions from Theorem 3.1.

5.   Conclusions

In this study, we proposed a numerical method for solving stochastic semilinear subdiffusion equations driven by fractionally integrated additive noise. The temporal discretization relies on Mittag-Leffler integrators, while the spatial discretization is based on the finite element method. The effectiveness of the proposed method was demonstrated through illustrative examples that provided support for the theoretical analysis. In the next phase, we aim to enhance efficiency by incorporating the stochastic fractional system model. We also plan to preserve important physical properties and structures, such as positivity preservation, maximum principle, long time behavior, and investigate singular solutions. This includes studying the uniform L1 long time behavior of time discretization for time-fractional partial differential equations with non-smooth data. Additionally, we will develop a finite volume scheme for two-dimensional time-fractional stochastic Fokker-Planck equations on distorted meshes, which preserves the maximum principle.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the NSFC under Grant No. 12147115, the University Natural Science Research Project of Anhui Province under Grant No. 2023AH050254, the Fundamental Research Program of Shanxi Province under Grant No. 202103021224317, the Lvliang High Tech Research and Development Program under Grant No. 2023GXYF14, and the present job was also partially supported by the School Research Project of AUFE, Grant No. ACKYC21049.

Conflict of interest

The authors declare no conflict of interest.