Error Analysis of the Nonuniform Alikhanov Scheme for the Fourth-Order Fractional Diffusion-Wave Equation

: This paper considers the numerical approximation to the fourth-order fractional diffusion-wave equation. Using a separation of variables, we can construct the exact solution for such a problem and then analyze its regularity. The obtained regularity result indicates that the solution behaves as a weak singularity at the initial time. Using the order reduction method, the fourth-order fractional diffusion-wave equation can be rewritten as a coupled system of low order, which is approximated by the nonuniform Alikhanov scheme in time and the finite difference method in space. Furthermore, the H 2 -norm stability result is obtained. With the help of this result and a priori bounds of the solution, an α -robust error estimate with optimal convergence order is derived. In order to further verify the accuracy of our theoretical analysis, some numerical results are provided.


Introduction
As of the last few decades, fractional partial differential equations (FPDEs) have gained widespread acceptance in both physics and engineering [1,2] to describe a variety of phenomena such as the universal response of electromagnetic fields [3], the denoising of images [4], etc.The exact solution to FPDEs can be constructed using tools such as Laplace transforms, Fourier transforms, and Mellin transforms; see [5,6].In reality, only a few FPDEs can obtain the exact solution.Thus, the numerical simulation of FPDEs has been greatly developed.A number of numerical methods have recently been developed for solving time-fractional diffusion equations with a weakly singular solution.The regularity result for this problem was presented by Stynes et al. [7], and the L1 scheme was constructed on graded meshes.For the nonlinear time-fractional diffusion equation, Liao et al. [8] proposed the nonuniform Alikhanov scheme and constructed the corresponding discrete fractional Gronwall inequality.A nonuniform L2 scheme was designed by Kopteka [9] and the local error of the formulated scheme was estimated.By using the Müntz polynomial, Hou et al. [10] proposed a method for solving the time-fractional diffusion equation using the fractional spectral method.
The purpose of this paper is to develop a numerical method that is efficient in solving the following fourth-order fractional diffusion-wave equation: u(0, x) = φ(x) and u t (0, x) = φ(x) for x ∈ Ω, (1b) where 1 < α < 2, κ > 0, c ≥ 0, Ω = [x l , x r ] × [y l , y r ] ⊂ R 2 , φ(x) ∈ C ∞ ( Ω), and φ(x) ∈ C ∞ ( Ω).Here, D α t is the standard Caputo fractional derivative defined by As a matter of fact, Model (1) can be used to describe a wide variety of phenomena in physics, chemistry, engineering, and medicine, such as surfactant and liquid delivery into the lung [11], ice formation [12], forming grooves on flat surfaces [13], removing noise from medical magnetic resonance images [14], and so on.The main advantage of Model (1) is that it is able to be used to describe the evolution processes between diffusion and wave propagation, which is one of its most significant advantages.A number of numerical methods have recently been developed for solving the fractional diffusion-wave equation [15,16].Despite this, it is relatively rare to simulate Equation ( 1) with a weakly singular solution numerically.A combination of the L1 scheme and the compact difference scheme is applied in the temporal and spatial directions by Hu and Zhang [17] to resolve Problem (1) by means of order reduction methods.As described in [18], Li and Vong solved Problem (1) by applying the parametric quintic spline in spatial direction and the weighted shifted Grünwald Letnikov scheme in the temporal direction.Bhrawy et al. [19] proposed an efficient and accurate spectral method to simulate Equation (1).To solve Equation (1) with nonlinear terms, Huang et al. [20] used piecewise rectangular formulas and the Crank-Nicolson technique in time, and unconditional convergence was derived using the energy method.A fast Alikhanov scheme was proposed by Ran and Lei [21] to approximate the Caputo derivative of Problem (1) with a variable coefficient, and the unconditional stability as well as the optimal convergence under the maximum norm of the scheme was demonstrated.
It should be noted that the above numerical methods are based on the assumption that the exact solution of Problem (1) is smooth.This assumption, as far as we know, is unrealistic, because at the initial time, the solution of time-fractional equations usually behaves as a weak singularity.Therefore, this paper establishes the regularity of u as well as the intermediate variable for Equation (1).By introducing two intermediate variables, w = D α/2 t (u − t φ) and p = ∆(u − t φ), we are able to rewrite a fourth-order fractional diffusion-wave Equation (1) and approximate its solution using nonuniform Alikhanov schemes in temporal direction and the finite difference method in the spatial direction.Furthermore, we provide an optimal convergence analysis of the proposed method based on the H 2 -norm which satisfies α-robust.
We organize this paper as follows.The regularity of u and the introduced variable w = D α/2 t (u − t φ) is analyzed in Section 2. In Section 3, an equivalent coupled system for Equation (1) is approximated using nonuniform Alikhanov schemes in time and the finite difference method in space.Moreover, Section 4 provides stability and convergence results in H 2 -norm for the proposed scheme.In Section 5, two numerical experiments are presented to verify the accuracy of theoretical results.The final section of the paper presents our conclusions.
In the rest of the paper, we note that generic constant C is independent of the mesh, and it takes different values depending on where it is used.Furthermore, this constant is α-robust, i.e., the constant's value remains finite as α → 2.

Regularity of the Solution
In this section, the regularity of the exact solution for original Problem (1) is analyzed.We denote ũ(t, x) := u(t, x) − t φ(x).Now, we reformulate the original Problem (1) as ũ(0, x) = φ(x) and ũt (0, x) = 0 for x ∈ Ω, (3b) Applying the separation of variables, we construct the solution to Problem (3).We let {(λ i , ϕ i ); i = 1, 2, . . .} be the eigenvalues and eigenfunctions of the following Sturm-Liouville problem: where ϕ i is the unit orthogonal basis of L 2 (Ω) for all i and 0 < λ From [22] (p.3327), we observe that ϕ i is the eigenfunctions for the eigenvalue problem, where eigenvalue λ i is defined by λ i := λ i (κ 2 λ i + c) for i = 1, 2, . . ., ∞.Now, applying the separation of variables to the equivalent Equation (3) yields where J i (t) is defined by Here, the classic Mittag-Leffler function E ρ,η can be represented as follows: Now, the following fractional power L γ is defined for each γ ∈ R by We set w = D α/2 t (u − t φ).Imitating [23] (Section 3), a priori bounds on ũ and w of Problem (1) are stated in the following lemma.
Then, the solution of (3) and the introduced variable w satisfy Proof.Applying the definition of the Caputo derivative and (4) yields where S i 1 and S i 2 are defined by Differentiating Series ( 4) and ( 6) with respect to t, x, and y, and then imitating the analysis of [23] (3.10) yields (5).

The Fully Discrete Method
In this section, we construct the fully discrete method for Problem (1).Now, we set β := α/2 in the remainder of the paper.From [15] (Lemma 2.1), the following significant property is stated to construct the numerical method.
Using Lemma 2, Equation ( 3) can be rewritten as w(0, x) = 0, p(0, x) = 0, and ũ(0, x) = φ(x) for x ∈ Ω, For n = 0, 1, . . ., N, we set t n = T(n/N) r , where N is a positive integer, and r represents the user's choice of grading parameter.We denote For any function v ∈ C[0, T] ∩ C 3 (0, T], the fractional derivative D β t v(t, x) of ( 7) can be calculated at t = t n−σ using the following nonuniform Alikhanov method: in which I 2,j v(η) is the quadratic polynomial interpolating at points t j−1 , t j and t j+1 .In (8), coefficients d σ n,n−l are defined as follows: where Now, the temporal truncation error of the nonuniform Alikhanov scheme is stated in the following two lemmas.
For spatial meshes, the uniform step sizes in each directions are defined as h x := (x r − x l )/M x and h y := (y r − y l )/M y , where M x and M y are positive constants.We let and ∂Ω h = Ω h ∩ ∂Ω, where x i := x l + ih x and y j := y l + jh y .We denote u n i,j or u n h as the nodal approximation to u(t n , x i , y j ) for all admissible i, j and n.
For each grid function v on Ω h , ∆v(t n , x i , y j ) can be calculated by where At each point (t n−σ , x i , y j ), applying the nonuniform Alikhanov scheme (8) and Scheme (9) to approximate the Caputo derivative and the Laplacian operator of System (7) yields our fully discrete Alikhanov scheme: for n = 1, . . ., N, where ũ0 h = φ(x), w 0 h = 0, and p 0 h = 0.

The Error Estimate of the Fully Discrete Alikhanov Scheme
In this section, the H 2 -norm stability result of the fully discrete Alikhanov scheme (10) is given.Moreover, we obtain an α-robust convergent result.
For any grid function v and η that vanish on ∂Ω, we define the discrete inner product, (v, η) := h x h y ∑ (x i ,y j )∈Ω h v(x i , y j )η(x i , y j ), the discrete L 2 -norm, ∥v∥ := (v, v), and the discrete for any v h and η h that vanish on ∂Ω.Now, an important property for the Alikhanov scheme ( 8) is stated.
Lemma 6 ([8] (Theorem 3.1)).We let the bounded sequences, { ξ n ≥ 0 : n ≥ 1} and { v n ≥ 0 : then, we have The next step is to consider the following general case, which includes (10) as a special case.For any 1 ≤ n ≤ N, we assume the pair of grid functions (ω n h , ν n h , µ n h ) that vanish on ∂Ω satisfy where ω n h , µ n h , and ν n h vanish on ∂Ω, and ζ n h , f n h , and ς n h are given grid functions on Ω h .Next, we present a useful bound of (15), which is applied to derive the stability result of the fully discrete Alikhanov scheme (10).
Proof.We fix n ∈ {1, 2, . . ., N}. Taking the inner product of Equations (15a)-(15c) by ω n,σ h , κ 2 ∆ h µ n,σ h , and κ 2 ∆ h ω n,σ h respectively, we have Adding the above three equations, we have where ( 11) is used.By taking the inner product of Equations ( 15b) and (15c) with ∆ h µ n,σ h and ω n,σ h , we have Adding the above two equations and applying (11), we have Substituting the above result into (17) yields Using a Cauchy-Schwarz inequality and Lemma 5, we have Applying Hölder inequality produces Furthermore, applying Lemma 6 yields We now present the stability result for the fully discrete Alikhanov scheme (10) in the following theorem.Theorem 1. Solution (w n h , ũn h , p n h ) of (10) satisfies Proof.It is clear that ( 10) is a special case of (15) with From ( 20), we obtain (19) immediately by invoking ( 12) with γ = β.
Next, we state the error analysis of our fully discrete scheme.First, we denote e n w := w n h − w(t n , x i , y j ), e n u := ũn h − ũ(t n , x i , y j ), and e n p := p n h − p(t n , x i , y j ).
Now, we begin to present the error equation of coupled system (10).For 1len ≤ N, we define From ( 7) and (10), we obtain the following error equations: where χ n−σ is defined by The optimal error estimate for Scheme (10) is given by combining the regularity results given in Section 2.
Theorem 2. We let l N = 1/lnN.At each time level t n , we let (w n , u n , p n ) and (w n h , u n h , p n h ) be the solutions of (7) and (10), respectively.Then, we have for 1 ≤ n ≤ N.
Proof.It is clear that System of error equations ( 21) is a special case of (15) with where e 0 u = 0, e 0 w = 0, and ∥v∥ ≤ C Ω ∥∆ h v∥ are used.Applying the triangle inequality and regularity Result (5), we obtain where Lemma 3 and Lemma 4 are used.Applying the Taylor expansion, we arrive at Hence, combining the above Taylor expansion with regularity Result (5a) yields Similarity to the estimation of ∥∆ h ρ s,σ u ∥, we have Substituting ( 24)-( 26) into (23) yields However, Utilizing the above bound and using (12) Invoking (12) with γ = β for the term h 2 x + h 2 y , we have Substituting ( 28) and ( 29) into (27) offers Finally, the desired Result ( 22) is obtained immediately.
Remark 1.In [23], the Caputo fractional derivative of Problem (1) is decomposed into D α−1 t v and u t by investigating the variable v = u t , then they are approximated by the nonuniform L1 scheme and the nonuniform Crank-Nicolson scheme respectively.Furthermore, the proposed numerical method achieves the min{rα/2, 3 − α} order in time, which is suboptimal.However, our convergent result given in Theorem 2 attains the 2 order in time.In addition, Theorem 2 shows that our error analysis is α-robust as α → 2 − .Nevertheless, the convergent result given in [15] contains the term Γ(1 − α/2), which blows up as α → 2 − .Corollary 1.We suppose r ≥ 4/α; then, we obtain that numerical solution ( ũn h , w n h ) satisfies for 1 ≤ n ≤ N.

Numerical Experiments
In this section, two numerical experiments with weakly singular solutions are presented to verify the theoretical result of our fully discrete Alikhanov scheme (10).
The theoretical solution of this numerical example is u(t, x, y) = (t α + t 3 ) sin x sin y, which displays weak singularity at t = 0.
In our following calculation, we estimate the global H 2 -norm error max 1≤n≤N ∥∆u n − ∆ h u n h ∥ of the Alikhanov scheme by taking r = 4/α.Table 1 shows the H 2 -norm error and the rate of the Alikhanov scheme for α = 1.2, 1.5, 1.8, where N = M x = M y is taken.The results displayed indicate that temporal convergence rates are N −2 , in agreement with Theorem 2. Table 2 shows the errors and their associated spatial orders for different α, where N = 200 and M x = M y are taken.Form Table 2, we observe O(h 2 ) convergence, again as predicted by Theorem 2.
Due to the fact that the solution of this numerical example is unknown, we use the two mesh principles given in [27] to verify the theoretical result of Theorem 2. Figure 1 displays the error max 1≤n≤N ∥∆u n − ∆ h u n h ∥ in time for different α, where N = M x = M y is chosen.It is shown that the rate of convergence attains the 2 order in time, as predicted by Theorem 2. Figure 2 presents the numerical solution at t n = 1 for α = 1.1, 1.9.The displayed result of this figure shows that the solution of this problem behaves diffusive as α → 1 + , and the characteristic of the initial solution φ(x, y) is well maintained as α → 2 − , which indicates that the propagation speed of waves is finite.Moreover, the attenuation of the numerical solution gradually slows down as α → 2 − , verifying the wave feature again.

Conclusions
This paper investigates the regularity result of Equation ( 1), which indicates that the solution behaves as a weak singularity at t = 0.In order to construct a fully discrete scheme, the nonuniform Alikhanov scheme in time is used along with the finite difference method in space.Then, the stability analysis and an α-robust optimal convergent result under H 2 -norm for the proposed scheme are obtained.Finally, two numerical examples are provided to verify the theoretical result.In the future, we will extend the proposed scheme to the Equation (1) with a nonlinear term and present the unconditional optimal error estimate.In addition, we will extend the proposed scheme to fractional equations with the Caputo-Hadamard fractional derivative.
Author Contributions: Conceptualization, methodology, investigation, writing-original draft preparation, Z.A.; validation, writing-review and editing, supervision, C.H.All authors have read and agreed to the published version of the manuscript.
Funding: The research of Chaobao Huang is supported by the National Natural Science Foundation of China under grant 12101360, the Support Plan for Outstanding Youth Innovation Team in Shandong Higher Education Institutions under grant 2022KJ184, and the Natural Science Foundation of Shandong Province under grant ZR2020QA031.

Table 1 .
Convergent results of max 1≤j≤N ∥∆u n − ∆ h u n h ∥ in temporal direction with r = 4/α.

Table 2 .
Convergent results of max 1≤j≤N ∥∆u n − ∆ h u n h ∥ in spatial direction.Let us solve the original Problem (1) with κ