Galerkin Finite Element Approximation of a Stochastic Semilinear Fractional Wave Equation Driven by Fractionally Integrated Additive Noise

: We investigate the application of the Galerkin ﬁnite element method to approximate a stochastic semilinear space–time fractional wave equation. The equation is driven by integrated additive noise, and the time fractional order α ∈ ( 1,2 ) . The existence of a unique solution of the problem is proved by using the Banach ﬁxed point theorem, and the spatial and temporal regularities of the solution are established. The noise is approximated with the piecewise constant function in time in order to obtain a stochastic regularized semilinear space–time wave equation which is then approximated using the Galerkin ﬁnite element method. The optimal error estimates are proved based on the various smoothing properties of the Mittag–Lefﬂer functions. Numerical examples are provided to demonstrate the consistency between the theoretical ﬁndings and the obtained numerical results.


Introduction
Consider the following stochastic semilinear space-time fractional wave equation driven by fractionally integrated additive noise, with 1 < α < 2, 1 2 where D is a bounded domain in R d , d = 1, 2, 3 with smooth boundary ∂D, and C 0 D α t u(t) and R 0 D −γ impact the internal energy [4]. For physical systems, stochastic perturbations arise from many natural sources, which cannot always be ignored. Therefore, it is necessary to include them in the corresponding deterministic model. It is not possible to find the analytic solution of the space-time fractional Equation (1). Therefore, one needs to introduce and analyze some efficient numerical methods for solving (1). Li et al. [5] considered the Galerkin finite element method of (1) for the linear case with the additive Gaussian noise, that is, f = 0 and γ = 0, and obtain the error estimates. In [6], the authors studied the Galerkin finite element method for approximating the semilinear stochastic time-tempered fractional wave equations with multiplicative Gaussian noise and additive fractional Gaussian noise, but they only established error estimates for α ∈ ( 3 2 , 2). Extensive theoretical results exist for the stochastic subdiffusion problem with α ∈ (0, 1), as seen in works such as [7][8][9][10][11][12], alongside corresponding numerical approximations in works including [13][14][15][16][17]. Regarding the theoretical and numerical findings for the stochastic wave equation, we recommend exploring references such as [18][19][20][21]. For theoretical advancements in fractional-order nonlinear differential equations, recent works such as [22][23][24][25][26][27] and their references provide a comprehensive overview.
In this paper, our focus lies on the application of the Galerkin finite element method to solve (1). Firstly, we establish the existence of a unique solution for (1) using the Banach fixed point theorem. Additionally, we analyze the spatial and temporal regularities of the solution. To approximate the noise, we employ a piecewise constant function in time, resulting in a stochastic regularized equation. This equation is then tackled using the Galerkin finite element method. We provide corresponding error estimates, utilizing the various smoothing properties exhibited by the Mittag-Leffler functions. We extend the error estimates in [5] from the linear case of (1) with Gaussian additive noise to the semilinear case with the more general integrated additive noise. We also extend the error estimates of [6] for the stochastic semilinear time fractional wave equation from α ∈ ( 3 2 , 2) to α ∈ (1, 2).
The paper is organized as follows. In Section 2, we provide some preliminaries and notations. In Section 3, we focus on the continuous problem and establish the existence, uniqueness, and regularity results for the problem (1). In Section 4, we discuss the approximation of the noise and obtain an error estimate for the regularized stochastic semilinear fractional superdiffusion problem. In Section 5, we consider the finite element approximation of the regularized problem and derive optimal error estimates. Finally, in Section 6, we present numerical experiments that validate our theoretical findings.
Throughout this paper, we denote C as a generic constant that is independent of the step size τ and the space step size h, which could be different at different occurrences. Additionally, we always assume > 0 is a small positive constant.

Lemma 1 ([28]
). (Itõ isometry property) Let ψ : [0, T] × Ω → H be a strongly measurable mapping such that t 0 E ψ(s) 2 ds < ∞. Let B(t) denote a real-valued standard Brownian motion. Then, the following isometry equality holds for t ∈ (0, T]: To represent the solution of (1) in the integral form, we utilize the Laplace transform technique to write down the solution representation in terms of the Mittag-Leffler functions. The Mittag-Leffler functions are defined in [28], and we use them to express the solution in a compact form.
The following Lemma is related to the bounds of the Mittag-Leffler functions.

Existence, Uniqueness, and Regularity Results
This section is dedicated to studying the existence, uniqueness, and regularity results of the mild solution of the stochastic semilinear space-time fractional superdiffusion model (1).

Assumption 1.
There is a positive constant C such that the nonlinear function f : and Assumption 2. The sequence {σ k (t)} with its derivative is uniformly bounded by µ k and γ k , respectively, i.e., where the series ∑ ∞ k=1 µ k and ∑ ∞ k=1 γ k are convergent.
Proof. Set C [0, T]; L 2 (Ω; H) λ , λ > 0, as the set of functions in C [0, T]; L 2 (Ω; H) with the following weighted norm u 2 For any fixed λ > 0 this norm is the same as the standard norm on C [0, T]; L 2 (Ω; H) . We can therefore define a nonlinear map T : For any λ > 0, the function u ∈ C [0, T]; L 2 (Ω; H) is a solution of (11) if and only if u is a fixed point of the map T : In order to apply the Banach fixed point theorem it suffices to show that for an appropriately chosen λ > 0, T is a contraction mapping. We first show that T u ∈ C [0, T]; L 2 (Ω; H) for any u ∈ C [0, T]; L 2 (Ω; H) . By (15) and the Cauchy-Schwarz inequality we obtain with Based on the smoothing properties of the solution operators and by Assumption 1, it follows that By the smoothing property of the operatorĒ α,β,γ , the isometry property of Brownian motion, we have with 0 ≤ r ≤ κ that, Note that v 1 , v 2 ∈ L 2 (Ω; H) and u ∈ C [0, T]; L 2 (Ω; H) ; we then obtain sup t∈[0,T] E T u 2 < ∞, which implies that T u ∈ C [0, T]; L 2 (Ω; H) . Now we consider the contraction property. For any given two functions u 1 and u 2 in C [0, T]; L 2 (Ω; H) λ , it follows from (15) with the smoothing property and boundedness ofĒ α,β,γ with p = 0 and Note that With α ∈ (1, 2), choose sufficiently large λ, we obtain Hence, T (u 1 ) − T (u 2 ) λ ≤ δ u 1 − u 2 λ . Therefore, by the Banach contraction mapping theorem, there exists a unique fixed point u of the nonlinear map T which is the mild solution given by (15) of the problem (1).
Assume that Assumptions 1-3 hold. Let v 1 ∈ L 2 (Ω;Ḣ q ), v 2 ∈ L 2 (Ω;Ḣ p ) with p, q ∈ [0, 2β]. Then, the following regularity results hold for the solution u of (11) with r ∈ [0, κ] and 0 ≤ p ≤ r ≤ 2β, 0 ≤ q ≤ r ≤ 2β, Proof. From the definition of the mild solution (11) and t ∈ (0, T] with 0 ≤ p, q ≤ r ≤ 2β, it follows that, with r ∈ [0, κ], For I 1 and I 2 , by the regularity Lemma in [5], we have and For I 3 using the smoothing property of the operatorĒ α,β (t) and the Assumption 1, we have For I 4 , by isometry property of the Brownian motion and Assumption 2 and the smoothing property of the operatorĒ α,β,γ (t), we arrive at, with 0 ≤ r ≤ κ, which implies that Hence the proof of the theorem is completed.
Proof. The proof is similar to the existence and uniqueness theorem; therefore, we will only indicate the changes in that proof. Set C [0, T]; L 2 (Ω;Ḣ 2β ) λ , t > 0 as the set of functions in C [0, T]; L 2 (Ω;Ḣ 2β ) with the following weighted norm: For the proof, it is now enough to show that the map T : λ is a contraction. We first show that T u ∈ C [0, T]; L 2 (Ω;Ḣ 2β ) for any u ∈ C [0, T]; L 2 (Ω;Ḣ 2β ) . By (15) and the Cauchy-Schwarz inequality, we obtain with u ∈ C [0, T]; L 2 (Ω;Ḣ 2β ) , By the smoothing properties of E α,β andẼ α,β with p = q, and using the Assumption 1, it follows that For the integral E| t 0Ē α,β,γ (t − s)dW(s)| 2 2β , a use of the isometry property and Assumptions 3 and 4 and the smoothing property of the operatorĒ α,β,γ , yields, with To resolve the integral 2Ē α,β,γ (s) 2 ds < ∞, it is enough to choose r = 2β, which means that k = r = 2β since 0 ≤ r ≤ κ. Hence, we need to restrict 2γ > 1 in order to obtain κ = 2β by Assumption 3. With such choices of κ and r and by noting that 1 2 < γ ≤ 1, we arrive at We ) . Next, we look at the contraction property of the mapping T . For any given two functions u 1 and u 2 in C [0, T]; L 2 (Ω;Ḣ 2β ) λ , it follows from (15) that Based on the same argument of the existence and uniqueness theorem proof, the rest of the proof follows and this concludes the proof.

Approximation of Fractionally Integrated Noise
Let 0 = t 1 < t 2 < · · · < t N < t N+1 = T be the discretization of [0, T] and t = T N be the time step size. The noise dβ k (s) ds can be approximated by using Euler method, 1), and N (0, 1) is the normally distributed random variable with mean 0 and variance 1. Assume that σ n k (s) is some approximation of σ k (s). To be able to obtain an approximation of here, χ i (t) is the characteristic function for the ith time step length [t i , t t i+1 ], i = 0, 1, 2, · · · , N − 1 and σ n k is some approximations of σ k . The following is the regularized stochastic space-time fractional superdiffusion problem. Let u n be an approximation of u defined by The solution of (32) takes the following form:

Assumption 5 ([5]
). Suppose that the coefficients σ n k (t) are generated in such a way that To regularize the noise dW n (s) ds , we need the following regularity assumption.
By following similar proofs as in Theorems 1 and 2, we can establish the following theorems for the approximate solution u n (t).
Then the following regularity result for the solution u n of Equation (33) holds with r ∈ [0, κ] and 0 ≤ q, p ≤ r ≤ 2β, (34) Suppose that Assumptions 1-6 hold. Let u and u n be the solutions of Equations (1) and (32), respectively. We have, for any given > 0, Proof. Subtracting (33) from (11), we obtain where By the definitions of dW and dW n , we now rewrite G 2 as We first estimate E G 1 2 . From the form of G 1 , using the smoothing property of the operatorĒ α,β (t − s) and Assumption 1, we arrive with 1 < α < 2 at For the estimate of E G 21 2 , using the Ito isometry property and Assumption 6, we obtain Note that, for α + γ < 3 2 , a use of the boundedness property of the Mittag-Lefler function yields Furthermore, for 1 < α + γ < 3, by using the asymptotic property of the Mittag-Lefler function, we have Thus, we now arrive at We now estimate G 22 . We first denote dβ m (s) and replace the variable s withs in the second term of G 22 . Using the orthogonality property of φ k , k = 1, 2, · · · , we obtain Thus, a use of the Cauchy-Schwarz inequality yields For I 1 , using the mean value theorem and the Assumption 5, we arrive at Now, following the same estimates as in (41) For I 2 , we note by the Mittage-Leffler function property that hence, Now we estimate s s (t − τ) α+γ−2 dτ for the different α and γ. We shall show that, with Case 1. We now consider the case α + γ < 3 2 . Ifs < s, then with 0 < < 1 2 , it implies that Since a θ − b θ ≤ (a − b) θ , for a > b > 0 and 0 < θ < 1, then for α + γ < 3 2 , and this implies that Similarly, we may show that for s <s, with 0 < < 1 2 , Therefore, we obtain, for α + γ < 3 2 , Case 2. Next, consider the case 3 2 ≤ α + γ < 3. Ifs < s then we obtain, Similarly, for s <s, it follows that Therefore, for 3 2 Note that Thus, we derive the following estimate for I 2 . For α + γ < 3 2 , For 3 2 ≤ α + γ < 3, Together with these estimates we obtain the following results.
An application of the Gronwall's Lemma completes the rest of the proof.

Finite Element Approximation and Error Analysis
Let D be the spatial domain and let T h be a shape regular and quasi-uniform triangulation of the domain D with spatial discretization parameter h = max K∈T h h K , where h K is the diameter of K. Let V h ⊂Ḣ β , 1 2 < β ≤ 1 be the piecewise linear finite element space with respect to the triangulation T h , that is, Let P h : L 2 (D) → V h and R h :Ḣ β → V h be the L 2 projection and fractional Ritz projection and are the eigenpairs of the discrete Laplacian, that is, forms an orthonormal basis of V h ⊂ H. Further, we introduce the following fractional discrete For χ ∈ V h , the discrete norm can be defined by |χ| 2 where v h 1 = P h v 1 , v h 2 = P h v 2 are chosen as L 2 projections of the initial functions v 1 , v 2 ∈ H. As it is in the continuous case, the solution of (57) takes the form where We have the following smoothing properties:

Lemma 8 ([5]). (Inverse Estimate in V h ) For any > s, there exists a constant C independent of h such that
We now consider the error estimates.
Suppose that Assumptions 1-6 hold. Let u n and u h n be the solutions of (32) and (57), respectively. Let v 1 , v 2 ∈ L 2 (Ω;Ḣ q ) with 0 ≤ q ≤ 2β. Then, there exists a positive constant C such that for any > 0, with r ∈ [0, κ] and 0 ≤ max(q, β) ≤ r ≤ 2β, Proof. Introducingũ h n (t) ∈ V h as a solution of an intermediate discrete system We split the error u h n (t) − u n (t) := (u h n (t) −ũ h n (t)) + (ũ h n (t) − u n (t)) := ζ(t) + η(t). Again using P h u n we split η(t), From Lemma 6 it follows that, with r ∈ [β, 2β], which means that To estimate θ, note that θ satisfies the following equation and hence, the representation of solution θ is written as Choose p = 0 and p = β separately, from Lemma 6 with γ = 0 and Lemma 7, it follows that for q = − 2β + p and 0 < < 2β that Now an application of regularity shows where we used the fact that β ds < ∞ since 0 < α 2β < 2 and 0 ≤ α( r−q β ) < 2. We now combine (63), (66), and (67) to arrive at an estimate for η as, with p = 0, β and 0 ≤ p ≤ r ≤ 2β, Now to estimate ζ, note that and therefore we now write ζ(t) in the integral form as Again, choose p = 0, β. From Lemma 6 with γ = 0 and Lemma 7, it follows for q = p and for 1 < α < 2, that Combining (68) and (71) it follows for p = 0 and β, and 0 ≤ p ≤ r ≤ 2β that An application of the Gronwall's Lemma completes the rest of the proof.
The main result of this paper is obtained by combining Theorems 6 and 7.

Remark 1.
In particular, when the noise is the trace class noise, i.e., In this case we have η n which are consistent with the results for the stochastic heat equation.

Numerical Simulations
In this section, we will explore three numerical examples of the stochastic semilinear fractional wave equation. For simplicity, we will focus on the Laplacian operator, that is, β = 1 in Equation (1). Our goal is to approximate Equation (1) with various functions f (u) and examine the experimentally determined orders of convergence in time. We consider both the cases of trace class and white noises. Specifically, we choose the following functions: By comparing the experimentally determined orders of convergence with the theoretical findings in Theorem 8, we observe consistent results, as expected. All the numerical simulations in this paper are performed on an Acer Aspire 5 Laptop.
To complete this, let us first introduce the numerical method for solving (8). Consider, with α ∈ (1, 2), where f (r), r ∈ R, u 0 (x) and u 1 (x) are given smooth functions. Here, with γ ∈ [0, 1], where β m (t), m = 1, 2, . . . are the Brownian motions. Here, e m (x) = √ 2 sin mπx denotes the eigenfunctions of the operator A = − ∂ 2 ∂x 2 with D(A) = H 1 0 (0, 1) ∩ H 2 (0, 1). Further, let (γ m , e m ), m = 1, 2, . . . be the eigenpairs of the covariance operator Q of the stochastic process W(t), that is, We shall consider two cases in our numerical simulations. Case 1: The white noise case, e.g., γ m = m −β with β = 0, which implies that where tr(Q) denotes the trace of the operator Q. Case 2: The trace class case, e.g., γ m = m −β with β > 1, which implies that The numerical methods for solving stochastic time fractional partial differential equations are similar to the numerical methods for solving deterministic time fractional partial differential equations. The only difference is that we have the extra term g in the stochastic case and we need to consider how to approximate g.
Let 0 < t 0 < t 1 < · · · < t N = T be a partition of the time interval [0, T] and τ the time step size. Let 0 = x 0 < x 1 < · · · < x M = 1 be a partition of the space interval [0, 1] and h the space step size.
We next consider how to calculate g n , which is more complicated than F n . Approximating the Riemann-Liouville fractional integral by the Lubich first-order convolution quadrature formula and truncating the noise term to M − 1 terms, we obtain the lth component of g n by, with l = 1, 2, . . . , M − 1, where w (−γ) j , j = 0, 1, 2, . . . , n are generated by the Lubich first-order method, with γ ∈ [0, 1], To solve (91), we first need to generate M − 1 Brownian motions β H m (t), m = 1, 2, . . . , M − 1. This can be performed by using MathWorks MATLAB function fbm1d.m [29], which gives the value of the fractional Brownian motion with the Hurst parameter H ∈ (0, 1) at any fixed time T. Let Nre f = 2 7 and T = 1 and let dtre f = T/Nre f denote the reference time step size. Let 0 = t 0 < t 1 < . . . Since we do not know the exact solution of the system, we shall use the reference time step size dtre f and the space step size h = 2 −7 to calculate the reference solution vre f . The spacial discretization is based on the linear finite element method.
We then choose kappa = 2 5 , 2 4 , 2 3 , 2 2 and consider the different time step size τ = dtre f * kappa to obtain the approximate solutions V n at t n = nτ.
Let us discuss how to calculate the lth component of g n in MATLAB. Denote The lth component of the vector g n satisfies g n (l) = w γ * dWdt, l = 1, 2, . . . , M − 1.
Finally, we shall consider how to calculate the L 2 projections P h u 0 and P h u 1 of u 0 and u 1 , respectively. Here, we only consider the case P h u 0 . The calculation of P h u 1 is similar. Assume that By the definition of P h , we obtain Hence, u 0 can be calculated by Example 1. Consider the following stochastic time fractional PDE (Partial Differential Equation) , where f (t, x) = x 2 (1 − x) 2 e t − (2 − 12x + 12x 2 )e t and the initial value u 0 (x) = x 2 (1 − x) 2 , u 1 (x) = x and g(t, x) is defined by (78).
Let v(t, x) = u(t, x) − u 0 (x) − tu 1 (x) and transform the system (93)-(95) of u into the system of v. We shall consider the approximation of v at T = 1. We choose the space step size h = 2 −6 and the time step size dtre f = 2 −7 to obtain the reference solution vref. To observe the time convergence orders, we consider the different time step sizes τ = kappa * dtre f with kappa = [2 5 , 2 4 , 2 3 , 2 2 ] to obtain the approximate solution V. We choose M1 = 20 simulations to calculate the following L2 error at T = 1 with the different time step sizes By Theorem 8, the convergence order should be In Table 1, we consider the case of trace class noise, where γ m = m −2 for m = 1, 2, . . . . We observe that the experimentally determined time convergence orders are consistent with our theoretical convergence orders, as indicated in the numbers in the brackets. We have included the CPU time in seconds for running 20 simulations in each experiment. The CPU times exhibit similarity across the other tables; hence, we have decided not to include them in subsequent tables.
In Table 2, we consider the case of white noise, where γ m = 1 for m = 1, 2, . . . . We observe that the experimentally determined time convergence orders are slightly lower than the orders in the trace class noise case, as we expected.
The reference line in Figure 1 is determined by four points (log 2(∆t j ), p log 2(∆t j )), j = 1, 2, 3, 4 and the blue line in Figure 1 is determined by the four points (log 2(∆t j ), log 2(e j )), j = 1, 2, 3, 4. If these two lines are parallel, then we may conclude that the experimentally determined convergence order is almost p. We use the similar approach to obtain other figures below.
In Figure 2, we plot the experimentally determined orders of convergence with γ = 0.6 and α = 1.1 as shown in Table 2 for the white noise. We observe that the convergence order is almost O(τ) in the figure, where the reference line represents the order O(τ).  Table 1.  Table 2.
where f (u) = sin(u) and the initial values u 0 ( We use the same notations as in Example 1. In Table 3, we consider the case of trace class noise, where γ m = m −2 for m = 1, 2, . . . . We observe that the experimentally determined time convergence orders are consistent with our theoretical convergence orders, as indicated in the numbers in the brackets.
In Table 4, we consider the case of white noise, where γ m = 1 for m = 1, 2, . . . . We observe that the experimentally determined time convergence orders are slightly lower than the orders in the trace class noise case, as expected.
In Figure 3, we plot the experimentally determined orders of convergence with γ = 0.6 and α = 1.6 for the trace class noise as shown in Table 3. The expected convergence order is O(τ min{1,α+γ−1/2} ) = O(τ). The reference line in the figure represents the order O(τ), which is consistent with our observations.
We use the same notations as in Example 1. In Table 5, we consider the trace class noise, i.e., γ m = m −2 , m = 1, 2, . . . , and observe that the experimentally determined time convergence orders are consistent with our theoretical convergence orders. The numbers in the brackets denote the theoretical convergence orders.
In Table 6, we consider the white noise, i.e., γ m = 1, m = 1, 2, . . . , and observe that the experimentally determined time convergence orders are slightly less than the orders in the trace class noise case, as expected.
In Figure 5, we plot the experimentally determined orders of convergence with γ = 0.6 and α = 1.6 in Table 5 for the trace class noise. The expected convergence order is O(τ min{1,α+γ−1/2} ) = O(τ). We indeed observe this in the figure where the reference line is for the order O(τ).
In Figure 6, we plot the experimentally determined orders of convergence with γ = 0.6 and α = 1.6 in Table 6 for the white noise. We observe that the convergence order is almost O(τ) in the figure where the reference line is for the order O(τ).  Table 5.  Table 6.

Conclusions
In this paper, we explore a numerical approach to approximate the stochastic semilinear space-time fractional wave equation. We establish the existence of a unique solution for this equation by using the Banach fixed point theorem, assuming that the nonlinear function satisfies the global Lipschitz condition. To obtain the stochastic regularized problem, we approximate the noise using a piecewise constant function in time. The finite element method is then employed to approximate the stochastic regularized problem. Furthermore, we propose a natural extension of this work, which involves considering more general nonlinear functions, such as the Allen-Cahn equation, as well as exploring different types of noise, including fractional noise with a Hurst parameter H ∈ (0, 1). Additionally, it might be very interesting to investigate more advanced fractal-fractional derivatives [30] in our stochastic space-time fractional wave equation.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.