Numerical Analysis of Local Discontinuous Galerkin Method for the Time-Fractional Fourth-Order Equation with Initial Singularity

: In this paper, efﬁcient methods seeking the numerical solution of a time-fractional fourth-order differential equation with Caputo’s derivative are derived. The solution of such a problem has a weak singularity near the initial time t = 0. The Caputo time-fractional derivative with derivative order α ∈ ( 0,1 ) is discretized by the well-known L1 formula on nonuniform meshes; for the spatial derivative, the local discontinuous Galerkin (LDG) ﬁnite element method is used. Based on the discrete fractional Gronwall’s inequality, we prove the stability of the proposed scheme and the optimal error estimate for the solution, i.e., ( 2 − α ) -order accurate in time and ( k + 1 ) -order accurate in space, when piece-wise polynomials of degree at most k are used. Moreover, a second-order and nonuniform time-stepping scheme is developed for the fractional model. The scheme uses the L 2-1 σ formula for the time fractional derivative and the LDG method for the space approximation. The stability and temporal optimal second-order convergence of the scheme are also shown. Finally, some numerical experiments are presented to conﬁrm the theoretical results.


Introduction
In this paper, we consider the local discontinuous Galerkin (LDG) finite element method for the following time-fractional fourth-order problem with periodic boundary condition C D α 0,t u + c 1 u x + c 2 u xx + c 3 u xxx + c 4 u xxxx = f (x, t), (x, t) ∈ D = (0, L) × (0, T], u(x, 0) = u 0 (x), x ∈ Ω = (0, L), where c 1 , c 2 , c 3 and c 4 > 0 are arbitrary constants. Without loss of generality, we assume that c 1 > 0 and c 3 > 0; however, we do not require the sign of c 2 to be positive or negative.
The source term f (x, t) and the initial value u 0 (x) are given functions. The term C D α 0,t u represents the Caputo fractional derivative of order α (0 < α < 1) with respect to t, which is [1,2], in which the operator ∂ s denotes the partial derivative with respect to s, and Γ(·) is the usual gamma function. Time-fractional partial differential equations with fourth-order spatial derivatives have been widely used in various fields, such as bridge slabs, airplane wings, floor system and window glass (e.g., [2][3][4][5]). On the assumption that the analytical solution is sufficiently smooth, many numerical methods have been devised for this kind of problem. In [6], a fully discrete LDG scheme was proposed to solve the time-fractional fourth-order equation and was proven to be stable and convergent with order O(τ −α h k+1 + τ 2−α + τ −α/2 h k+1/2 + h k+1 ),

Preliminaries
Let us start by presenting the notations, definitions and projections used in this paper.

Tessellation and Function Space
Let T h = {I j = (x j− 1 2 , x j+ 1 2 )} N j=1 be the partition of Ω, where x 1 2 = 0 and x N+ 1 2 = L are the two boundary endpoints. For each cell I j , the cell center and cell length are denoted by x j = (x j− 1 2 + x j+ 1 2 )/2 and h j = x j+ 1 2 − x j− 1 2 , respectively. We use h = max j h j to denote the length of the largest cell. Throughout this paper, it is assumed that T h is a quasi-uniform mesh; namely, there exists a fixed positive constant ρ independent of h such that ρh ≤ h j ≤ h for any j = 1, . . . , N, when h goes to zero. Define the finite element space V h = v ∈ L 2 (Ω) : v| I j ∈ P k (I j ), ∀j = 1, . . . , N , where P k (I j ) denotes the space of polynomials of degree no more than k on the cell I j . Note that the functions in this space are allowed to have discontinuities across element interfaces. As usual, we use u − . In particular, if l = 0, we use an unmarked norm · to represent the usual L 2 norm on Ω.
To end this subsection, we list some inverse properties of the finite element space V h . For any v h ∈ V h , there exists a positive constants µ independent of v h and h, such that Here and below, Γ h is the union of all cell boundary points, and for any w ∈ H 1 (T h ), the L 2 norm on Γ h is defined by

The Semi-Discrete LDG Scheme
As usual, we first introduce some auxiliary variables approximating various order derivatives of the solution and rewrite Equation (1) into a first-order system, Then, the semi-discrete LDG scheme is as follows: ∀t ∈ (0, T], find u h , q h , p h , r h ∈ V h , such that for any v h , ρ h , φ h , ψ h ∈ V h and j = 1, . . . , N, it holds that Notice that the periodic boundary conditions are considered, i.e., ζ − 1 2 Denote, by (v, w) = ∑ N j=1 I j vwdx, the inner product in L 2 (Ω). Summing up the variational formulations (5) over j = 1, 2, . . . , N, we can write the semi-discrete LDG scheme (5a)-(5d) in the global form: Here, Using the definitions of the above operators, the following lemmas can be obtained, and the proof is straightforward (refer to [23] where C µ is a positive constant that is independent of h but may depend on the inverse constant µ.
Lemma 4 presents the important relationships between the primal and auxiliary variables, which are crucial to the stability analysis.

Projection and Interpolation Property
In what follows, we define the elliptic projection. For any function u, q = u x , p = q x , and r = p x , the elliptic projection is the unique solution In addition, since u h is determined as an additive constant in the elliptic problem with periodic boundary conditions, to ensure that (12a)-(12d) is well-defined, we assume [24] (u − U h , 1) = 0.
By referring to Lemma 4.2 in [23], it can be found that the elliptic projection defined above exists uniquely and satisfies the following approximation properties.
Let U h , Q h , P h , R h ∈ V h be the elliptic projection (12), we have where C is a constant depending on the regularity of u but is independent of h.

Nonuniform L1-LDG Scheme
In this section, we propose a fully discrete numerical scheme to solve the timefractional fourth-order Equation (1), hereafter referred to as the nonuniform L1-LDG scheme, which discretizes the Caputo time-fractional derivative using the L1 formula on nonuniform meshes and the LDG method to discretize the spatial derivative.

The Fully Discrete Numerical Scheme and Its Stability
For a given finite time T > 0, denote t n = T(n/M) r , and let n = 0, 1, . . . , M be the mesh points, where r ≥ 1. Let τ n = t n − t n−1 , n = 1, . . . , M be the time mesh sizes. If r = 1, then the mesh is uniform.

Lemma 6 ([22]
). Assume that ∂ l t u(x, t) ≤ Ct l−α for l = 0, 1, 2. Then, Lemma 7 ([19]). Assume that u(x, ·) ∈ C 2 ((0, T]) and ∂ l t u(x, t) ≤ Ct l−α for l = 0, 1, 2, and then Let u n h , q n h , p n h , r n h ∈ V h be the approximation of u(x, t n ), q(x, t n ), p(x, t n ), r(x, t n ), respectively. Then, the fully discrete nonuniform L1-LDG scheme for problem (1) is as follows: find u n h , q n h , p n h , r n h ∈ V h such that, for any Now, we turn to the stability analysis of scheme (18). We first introduce the following discrete fractional Gronwall inequality.

Lemma 9 ([26]
). Let the functions u n = u(x, t n ) be in L 2 (Ω) for n = 0, 1, . . . , M. Then, one has the following inequality Theorem 1. If the graded mesh satisfies the maximum time-step condition , then, for n = 1, . . . , M − 1, the solution u n h of the fully discrete nonuniform L1-LDG scheme (18) satisfies Proof. Choosing the test function in (18a) as v h = u n h and using Lemma 2, we obtain It follows from Lemma 1 that By the Cauchy-Schwarz inequality and (18c), one has Applying Lemma 1, the equality (18d) and Lemma 4, we find Substituting (22)- (24) into (21), we arrive at Therefore, if we take = 4 3 c 4 and use the Cauchy-Schwarz inequality again, then By using the Young's inequality together with Lemma 9, we obtain the estimate Therefore, applying Lemma 8 with v n = u n h , φ n = 0, ψ n = f n , λ 0 = . This completes the proof.

Error Estimate of the Nonuniform L1-LDG Scheme
We are now ready to show the optimal error estimate of scheme (18). Assume that the solution u of time-fractional fourth-order problem (1) satisfies Theorem 2. Let u n be the exact solution of Equation (1) that satisfies the smoothness assumption (28), and u n h be the numerical solution of the nonuniform L1-LDG scheme (18). Then, for n = 1, 2, . . . , M, the following estimate holds where C is a positive constant independent of M and h.
Proof. For any t > 0, denote (e n u , e n q , e n p , e n r ) = (u n − u n h , q n − q n h , p n − p n h , r n − r n h ).
Let (U n h , Q n h , P n h , R n h ) be the elliptic projection defined in (12) at time t = t n . Then, we divide the error e n ς in the form for ς = u, q, p, r, where In order to obtain the error equation of the fully discrete numerical scheme, we need to present the weak formulation of (4) at t n , which is, where ϕ n = (u n , q n , p n ) and v, ρ, φ, ψ are test functions. Then, we can find the error equation by subtracting (18) from (33) that, for any v h , ρ h , φ h , ψ h ∈ V h and n = 1, . . . , M, where e n ϕ = (e n u , e n q , e n p ). According to the definition (12) of elliptic projection and (33b)-(33d), we obtain Therefore, we have from (34a) and (35a) that Here, ξ n ϕ = (ξ n u , ξ n q , ξ n p ) and η n ϕ = (η n u , η n q , η n p ). By using (8), (35c) and (35d), the above equation can be further written as Substituting (31) into (34b)-(34d) and using (35b)-(35d), we obtain Taking v h = ξ n u in (37a), we obtain the following identity From Lemma 2, we arrive at It has been shown in ( [23], Lemma 9) that By the similar argument to prove inequality (11) in Lemma 3 (refer ( [27], Lemma 3.3) for similar analysis), we can find Then, combine (40)-(41) and Lemma 4 to obtain where we have used the interpolating property (14) in the first inequality.
From Theorem 2, it can be concluded that the optimal order of convergence (i.e., (2 − α)th-order accurate in time and (k + 1)-order accurate in space) for the solution can be obtained if we use nonuniform L1 formula in the temporal direction and the LDG method in space. However, the numerical solution generated by scheme (18) will be limited to being (2 − α)-order accurate in time even if the solution is sufficiently smooth. Therefore, developing high-order numerical algorithms for the time-fractional fourth-order problem (1) is also indispensable and will be studied in the next section.

Nonuniform L2-1 σ -LDG Scheme
In the section, based on the LDG method in the spatial direction and L2-1 σ approximation in the time direction, we propose a fully discrete numerical scheme (called the nonuniform L2-1 σ -LDG scheme for brevity) with high temporal accuracy to solve the time-fractional fourth-order Equation (1).

Moreover, it was proven in [25] that
where π A is a positive constant. Let q = u x , p = q x , r = p x , and then the weak form of the time-fractional fourth-order Equation (1) at t n+σ is formulated as (r n,σ , ρ) = −H − (p n,σ , ρ), (57b) where v, ρ, φ, ψ are test functions; ϕ n,σ = (u n,σ , q n,σ , p n,σ ); the bilinear operators H ± and L are defined in (6) and (9), respectively; H(ϕ n,σ ; v) has been given in (8); and The LDG method introduced in Section 2 is used for spatial discretization. Then, the fully discrete nonuniform L2-1 σ -LDG approximation scheme for (1) reads as: find u n,σ h , Below, we study the stability of the nonuniform L2-1 σ -LDG scheme (58). The following lemmas play a key role in proving the stability for the nonuniform meshes.

Theorem 3. If the graded mesh satisfies the maximum time-step condition
Proof. Taking the test function v h = u n,σ h in (58a) and applying Lemma 2, we obtain

Numerical Examples
The purpose of this section is to numerically validate the efficiency of schemes (18) and (58) for solving the time-fractional fourth-order Equation (1) with initial singularity. All the algorithms were implemented using MATLAB R2016a, and were run in a 3.10 GHz PC with 16 GB RAM and a Windows 10 operating system. Example 1. Consider the problem (1) with c 1 = c 3 = c 4 = 1, c 2 = −1, Ω = (0, 1), T = 1, u 0 (x) = 0, and the periodic boundary condition in use. In this case, the source term The analytical solution is given by u(x, t) = (t α + t 3 ) sin(2πx). This solution displays a weak singularity at t = 0.
To solve Example 1, we apply the nonuniform L1-LDG scheme (18) with r = (2 − α)/α in computation. The L 2 -norm errors and temporal convergence orders of the numerical solution u n h for α = 0.4, 0.6, 0.8 and different time-steps are listed in Table 1. The convergence orders of α = 0.4 and 0.6 are close to (2 − α), which is consistent with the theoretical prediction in Theorem 2. However, the accuracy of α = 0.8 is slightly lower. In Tables 2 and 3, for fixed M = 4000 and r = (2 − α)/α, we observe that the spatial convergence order for (18) is (k + 1), which is in agreement with the theoretical analysis. The numerical solutions of the scheme (18) for different α are given in Figures 1-3. Figure 4 depicts the L 2 -norm errors versus N between the numerical solution and the exact solution for different α at t = 1. The graphs show good agreement between the two solutions.       Example 2. The purpose of this example is to investigate the accuracy and efficiency of the proposed nonuniform L2-1 σ -LDG method (57). For simplicity, the equation in Example 1 is still regarded as a test problem, but it is solved by the scheme (57). The L 2 -norm errors at time t = 1 and convergence orders in the temporal direction with different α and r are shown in Tables 4-6. The orders of convergence displayed indicate that the order of convergence is O(M − min{rα,2} ), which coincides with Theorem 4. From Tables 5 and 6, we can also see that the grading parameter r ≥ 2/α yields the temporal optimal second-order accuracy. Then, we refine the spatial step size with a fixed temporal step size M = 2000. The L 2 -norm errors at time t = 1 and convergence orders in the spatial direction are shown in Table 7. The results imply that the algorithm (57) has an accuracy of O(h k+1 ) in space.

Concluding Remarks
In this paper, we studied the numerical algorithms for the time-fractional fourthorder equation with an initial singularity. In the temporal direction, two types of finite difference schemes were proposed and analyzed, including the nonuniform L1 scheme and nonuniform L2-1 σ scheme. In the spatial direction, the LDG method was utilized. Detailed proofs of L 2 stability and optimal error estimates for the schemes were derived using the discrete fractional Gronwall-type inequalities. Finally, some numerical examples were presented to verify the theoretical results. Data Availability Statement: All the data were computed using our algorithms.