Analyzing the continuity of the mild solution in ﬁnite element analysis of semilinear stochastic subdi ﬀ usion problems

: This paper aimed to further introduce the ﬁnite element analysis of non-smooth data for semilinear stochastic subdi ﬀ usion problems driven by fractionally integrated additive noise. The mild solution of this stochastic model consisted of three di ﬀ erent Mittag-Le ﬄ er functions. We analyzed the smoothness of the solution and utilized complex integration to approximate the error of the solution operator under non-smooth data. Consequently, optimal convergence estimates were obtained, and we also obtained the continuity conditions of the mild solution. Finally, the inﬂuence of the fractional parameters α and γ on the convergence rates were accurately demonstrated through numerical examples.


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 I γ t Ẇ(t) using the classical Wiener process.Here, 0 I γ t u (or RL D −γ 0,t u) denotes the Riemann-Liouville fractional integral of order γ of the function u defined by 0 I γ t u(t) ≡ RL D −γ 0,t u(t) = 1 Γ(γ) t 0 (t − σ) γ−1 u(σ)dσ.
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 positivitypreserving finite volume scheme tailored for subdiffusion equations on nonconforming quadrilateral distorted meshes with hanging nodes.Furthermore, [20] addresses the numerical solution of the threedimensional 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), E(t), and 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 D 1−α 0,t := ( 0 I α t ) 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 α and γ 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.

It is known that if
and the eigenvectors {ϕ k } ∞ k=1 form an orthonormal basis for H. Let H = L 2 (D) be a separable Hilbert space with inner product (•, •) and norm • .Let (Ω, F , P, {F t } t≥0 ) be a filtered probability space, with Bochner spaces L p (Ω; H) = L p ((Ω, F , P); H).Let E denote the expectation (with respect to 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 k=1 is a sequence of independently and identically distributed standard Brownian motions.
For any ν ∈ R, we introduce the space Ḣν (D) = D(A k=1 is the orthonormal basis in H. Let L = L(H) denote the space of all bounded linear operators on H and let We also need to recall the Burkholder-Davis-Gundy inequality [11], for p ≥ 2, for strongly measurable functions φ : [0, T ] → L 0 2 .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.
Throughout we always make the following standing assumption on the fractional orders α and γ, which is sufficient to ensure the well-posedness of Eq (1.1) (see [2,5]).
Lemma 2.2.There exists C such that for t > 0, we have (2.10) (2.11) (2.12) 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.
In the qualitative theory of nonlinear PDEs, the subsequent Gronwall type inequalities play a very important role in error estimate.
N , and t n = nk for 0 ≤ n ≤ N. If ζ 1 , ..., ζ N ≥ 0 satisfy for some M 0 , M 1 ≥ 0, and µ, ν > 0, the inequality To mimic Assumption 2.14 as stated in book [7], let P T be the σ-field of predictable stochastic processes, and B(S ) be the Borel σ-field of S .Given the available options, we shall make some reasonable assumptions about the nonlinear parts. 1  2 ), there exists a constant C that satisfies the following expression: )

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 T h be a regular shaped quasi-uniform triangulation of the domain D, and let S h ⊂ H 1 0 (D) be the space of continuous piecewise linear functions on the triangulation T h .We define the L 2 -projection P h : H → S h by and the Ritz projection R h : where a(u, χ) = (∇u, ∇χ) is the associated bilinear form.Note that by interpreting the righthand side of (3.1) as a duality pairing between Ḣ1 (D) and Ḣ−1 (D), one may extend P h to be a bounded operator from Ḣ−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) ∈ S h such that where the discrete Laplacian A h is defined by Naturally, we present the discrete analogues of operators E h (t), E h (t), and E h (t) as follows Hence, the fully discrete approximation of (1.1) is given by with initial value U 0 h = P h u 0 .Let us introduce and prove some Lemmas that will play an important role later on.Lemma 3.2.[10] Let 0 ≤ ω ≤ µ ≤ 2. For α ∈ (0, 1), there exists a constant C such that We now turn to the non-smooth data error estimates of the approximations to E(t)g, g ∈ H in the semi-discrete case.Lemma 3.3.Let 0 ≤ s ≤ 1 and 0 ≤ r ≤ 2 with r + s ≤ 2. For g ∈ H, there holds Proof.In the case s = 0, by the inverse Laplace transform, for any given g ∈ H, we have where Let us first show (3.7).For any fixed g ∈ H, we have, by (3.8) and (3.9), By [4], p. 820, we get (z  For I, with z = δe iφ , we have δ = t −1 , where t −1 < π τ for sufficiently small τ, For II, with z = re ±iθ , we have r ≥ δ, δ = 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 E(t).More details can be found in Lemma 4.4 of [5].
Lemma 3.4.Let 0 ≤ s ≤ 1 and 0 ≤ r ≤ 2 with r + s ≤ 2. For g ∈ H, there holds Based on the previous discussion, we are ready to prove the error estimates for the fully discrete approximation.
Under the assumptions of smoothness of the solution operator and the nonlinear term, we can obtain the following conclusions.Theorem 3.2.Let p ∈ [2, ∞) be given such that Assumptions 2.13-2.17hold, then the unique mild solution u to (2.3) is continuous with respect to the norm • L p (Ω; Ḣs ) .
Proof.According to the definition of continuity, our goal is to show that lim with either t 1 or t 2 fixed.Therefore, with the expression of a mild solution u, we divide u (t 2 ) − u (t 1 ) L p( Ω; Ḣs ) into five parts, and only need to show that each of the parts approaches zero when t 2 → t 1 .By utilizing the triangle inequality, we obtain u(t 1 ) − u(t 2 )
For the expression M 1 , we apply (2.11) from Lemma 2.2.One yields 1 ), and the validity of lim For M 2 , we insert a node t 3 and split it into three parts.
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 1+γ−s 2 α > 0. For M 3 , we prove that the following property holds first.One obtains (1 + t 4 λ j ) 4 (v, ϕ j ) 2 , and then we have u(σ) L p (Ω; Ḣr ) ).Since 1+r−s 2 α > 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 − s 2 )α + γ − 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 − s 2 )α + γ − 1 2 > 0. Thus, the theorem is proved.

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 α ∈ (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 = Furthermore, we assume that W(t) has the following Fourier series expansion: Thus, the semi-discrete solution U N m,k,h satisfies where u 0,k = (u 0 , ϕ k ), F k (•) = (F(•), ϕ k ), and β k (σ) for k = 1, . . ., 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 ∆t.To approximate the exact solution, we used the fully discrete solution with a small time step size of ∆t = 2 −10 and a spatial step size of ∆x = 2 −10 .The theoretical rate of convergence is αν for ν ∈ (0, 1), which approaches α as ν → 1.
To conduct the experiments, we fixed a small spatial step size ∆x and considered a sequence of moderate time step sizes ∆t i = 2 −i , i = 2, . . ., 5. We performed M = 100 simulations for each time step size ∆t i .In each simulation ω j , j = 1, 2, . . ., M, we generated N h independent Brownian motions β k (t), k = 1, 2, . . ., 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 α and γ.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 α and γ.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 α and γ, 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.

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.