A Mixed Element Algorithm Based on the Modiﬁed L 1 Crank–Nicolson Scheme for a Nonlinear Fourth-Order Fractional Diffusion-Wave Model

: In this article, a new mixed ﬁnite element (MFE) algorithm is presented and developed to ﬁnd the numerical solution of a two-dimensional nonlinear fourth-order Riemann–Liouville fractional diffusion-wave equation. By introducing two auxiliary variables and using a particular technique, a new coupled system with three equations is constructed. Compared to the previous space–time high-order model, the derived system is a lower coupled equation with lower time derivatives and second-order space derivatives, which can be approximated by using many time discrete schemes. Here, the second-order Crank–Nicolson scheme with the modiﬁed L 1-formula is used to approximate the time direction, while the space direction is approximated by the new MFE method. Analyses of the stability and optimal L 2 error estimates are performed and the feasibility is validated by the calculated data.


Introduction
Fourth-order fractional partial differential equations (PDEs) including fourth-order fractional subdiffusion models [1][2][3] and fourth-order fractional diffusion-wave models [2,4,5] can be founded in many fields of science and engineering. Thus far, there have been many efficient numerical algorithms for solving linear or nonlinear fourth-order fractional subdiffusion and diffusion-wave models. Liu et al. [6], Liu et al. [7], and Liu et al. [8] considered different mixed element methods to solve fourth-order nonlinear fractional subdiffusion models with the first-order time derivative and developed numerical theories including stability and convergence. Liu et al. [3] introduced a mixed element algorithm with a new approximation of the fractional derivative. Ji et al. [9], Ran et al. [10], Nandal and Pandey [11], Sun et al. [12], and Huang et al. [13] considered some difference schemes for linear or nonlinear fourth-order fractional diffusion or diffusion-wave models. Abbaszadeh and Dehghan [14] studied the direct meshless local Petrov-Galerkin method for solving fourth-order reaction-diffusion problems with a time-fractional derivative. Yang et al. [15] and Zhang et al. [16] found the numerical solutions for a fourth-order fractional model by using the orthogonal spline collocation method. Tariq and Akram [17] considered a quintic spline technique to solve a fourth-order time-fractional subdiffusion model. Guo et al. [18] and Du et al. [19] studied the LDG methods for solving some timefractional subdiffusion models with fourth-order spatial derivative terms, respectively. In [1], Nikan et al. developed a local radial basis function generated by the finite difference scheme for a time-fractional fourth-order reaction-diffusion model. In [5], Jafari et al. solved a fourth-order fractional diffusion-wave equation by the decomposition method. Hu and Zhang [2] implemented numerical calculations via finite difference methods for fourthorder time fractional subdiffusion and diffusion-wave models. Li and Wong [20] developed an efficient numerical algorithm for a fourth-order time-fractional diffusion-wave model.
Here, we propose a new mixed element algorithm to solve the following nonlinear fourth-order time-fractional diffusion-wave model: where Ω ⊂ R d (d ≤ 2) and J = (0, T] with 0 < T < ∞ are the spatial domain and time interval, respectively. u 1 (x) is an initial value function, g(x, t) is a given source term, f (u) is a polynomial function or bounded function on u satisfying f ∈ C 2 (R), and the Riemann-Liouville fractional derivative is defined by where the nonlinear fourth-order fractional diffusion-wave model (1) can be generated by the classical fourth-order hyperbolic wave equation. When β → 1 or 2 and f (u) = u 3 − u, the model (1) can be reduced to an important Cahn-Hilliard equation model [21].
Recently, Zeng and Li [22] developed a new Crank-Nicolson scheme based on a modified L1-formula, whose coefficients are different from the famous L1-formula (see [23,24] for the fractional parameter α ∈ (0, 1)). One should note that this modified L1-formula can only approximate the Caputo or Riemann-Liouville fractional derivative with parameter α ∈ (0, 1), and it cannot approximate the case β ∈ (1, 2). Here, we will develop the modified L1-formula for the case of β ∈ (1, 2) by using some techniques.
In this article, by introducing two auxiliary functions and using some techniques, we propose a new mixed element algorithm. Here, our major contributions are as follows: (1) by the introduction of two auxiliary functions, we reduce the nonlinear fourth-order time-fractional diffusion-wave model to a low-order coupled system; (2) we turn order β ∈ (1, 2) into order α ∈ (0, 1) for the Riemann-Liouville fractional derivative; (3) we approximate the derived coupled system with a fractional derivative with order α ∈ (0, 1) by the modified L1 Crank-Nicolson scheme with the developed new mixed element method; (4) we derive the stability of the new mixed element scheme and optimal error estimates in the L 2 -norm for three functions.
The structure of this article is as follows: in Section 2, we provide some numerical approximation formulas, propose a new mixed element scheme, and prove the stability of the derived scheme; in Section 3, we derive optimal error estimates for three variables; in Section 4, some numerical data are computed and discussed; Finally, in Section 5, we give some concluding remarks.

Numerical Approximation and Stability
Based on the relation between the Riemann-Liouville fractional derivative and Caputo fractional derivative, we take α = β − 1 and v = ∂u ∂t to obtain Let σ = u − f (u); (1) can be rewritten as the following coupled system: For formulating the fully discrete scheme, we insert the nodes t n = n∆t(n = 0, 1, 2, · · · , N) in time interval [0, T], where t n satisfy 0 = t 0 < t 1 < t 2 < · · · < t N = T and the time step length size ∆t = T/N, for some positive integer N. For a smooth function φ defined on the time interval [0, T], we denote φ n = φ(t n ). Now, we need to introduce some lemmas on integer and fractional derivatives.
, the following relation holds: Lemma 3 ([22]). At t k+ 1 2 , the Caputo fractional derivative has the following form: for k ≥ 0, we have , the Riemann-Liouville fractional derivative has the following approximation:

Remark 1.
In [22], the authors provided the L1-formula above, which is different from the usual L1-formula and called the modified L1-formula.
By the approximation scheme above, we arrive at where we have the following mixed weak formulation: , based on the mixed weak formulation above, we formulate the following new mixed element system: Lemma 5 (See [22]). For b n−j , the following important inequality holds: where C is a positive constant that is independent of space-time step length sizes h and ∆t.
Proof. Applying the Taylor formula, we have Noting that ∆t = T N and n − j ≤ N, we have Substitute (15) into (14) to obtain the conclusion.
Next, we will prove the stability.

A Priori Error Estimate
Now, we provide two projection operators [25] to derive a priori error estimates of our mixed finite element method. Lemma 6. Define the L 2 projection P h : L 2 (Ω) → L h as with the estimate inequality u − P h u + u t − P h u t ≤ Ch m+1 u m+1 , ∀u ∈ L 2 (Ω).
(34) Lemma 7. Define the elliptic projection Q h : with the following inequality: In what follows, we derive the proof of error estimates in L 2 -norm in detail.
there exists a positive constant C that is independent of space-time step length sizes (h, ∆t) and we have for n ≥ 0 where, for the Caputo fractional derivative, we take µ as 0; for the Riemann-Liouville fractional derivative, we take µ as 1.
Proof. For convenience, we write Applying triangle inequality, we have Using Lemmas 6 and 7, we arrive at the estimates of E n , F n , and H n . Consequently, in the discussion below, we only need to derive the estimates of E n , F n , and H n . Using projections (33) and (35), we have error equations as follows: In (39), we set ϕ h = E n+ 1 2 , χ h = F n+ 1 2 , and ψ h = H n+ 1 2 , and add the resulting equations to obtain (a n−j − a n−j+1 )F j− 1 2 − (a n − b n )F Now, we need to estimate all terms on the right-hand side of (40). Using Cauchy-Schwarz inequality, we have Applying the mean value theorem and Cauchy-Schwarz inequality, we have where we use the boundedness of f u (u n ) ∞ and the following bounded inequality: where one can apply inverse inequality [25], and use a similar method as the one in [7,26]. Making use of (9), (3), Cauchy-Schwarz inequality, as well as Young inequality, we have Making a combination for (41)-(44) and using (18), we have With given conditions E 0 = 0, F 0 = 0, H 0 = 0, we use (23) to arrive at Sum for (46) with respect to n to arrive at For n = 0, we use a similar derivation to the one of n ≥ 1 and apply triangle inequality to arrive at Substitute (48) into (47) and use the Gronwall lemma to obtain Combining (49), (34), and (36) with (38) and noting that α = β − 1, we complete the proof of the theorem.

Remark 2.
Compared with the classical mixed element method for fourth-order partial differential equations, our method can approximate simultaneously three variables with optimal error estimates in L 2 -norm. More importantly, we can obtain directly optimal error estimates in L 2 -norm for auxiliary variables in solving fourth-order PDEs, which are difficult to achieve by using classical mixed element methods [6][7][8].

Numerical Tests
Here, we will verify the theoretical results by numerical computing. In (1), we take space domainΩ = [0, 1] 2 , time intervalJ = [0, 1], nonlinear term f (u) = u 3 − u, initial conditions with u(x, y, 0) = 0, u 1 (x, y) = 0, and exact solution u = t 3 sin(2πx) sin(2πy); we can obtain the source term g(x, y, t) and two auxiliary variables v = 3t 2 sin(2πx) sin(2πy) and σ = −t 3 sin(2πx) sin(2πy)(8π 2 + t 6 sin(2πx) 2 sin(2πy) 2 − 1). In the following numer-ical calculations, the order of convergence in space is calculated by the following formula with a sufficiently small time step size ∆t where h k (k = 1, 2) represents different space mesh step lengths. For implementing the new mixed element algorithm, we approximate the spatial direction by the finite element method with the basis function P(x, y) = a + bx + cy + dxy and discretize the time direction by using the modified L1 Crank-Nicolson scheme. In Table 1, by taking the fixed time mesh parameter ∆t = 1/200, changed spatial step length sizes h = √ 2/9, √ 2/16 and √ 2/25, and different parameters β = 1.1, 1.5, 1.9, we show the L 2 -norm error estimates and second-order convergence data in space. In Tables 2 and 3, we compute the convergence results v and σ, respectively. From Tables 1-3, one can see that the numerical method is effective for solving nonlinear fourth-order fractional diffusion-wave equation models with a smooth solution.

Concluding Remarks
From the calculated results in Tables 1-3, one can see that our method for solving fourth-order fractional diffusion-wave equations in this article can obtain optimal error estimates in L 2 -norm for three variables, which is in agreement with the derived theoretical results. These results for auxiliary variables are difficult to achieve directly by using classical mixed element methods [6][7][8].
In the future, we will improve our mixed element method by combining other techniques [7,27,28] with high-order time approximate schemes and develop their optimal numerical theories. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All the data are computed by using our mixed element method.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.