Some Second-Order σ Schemes Combined with an H1-Galerkin MFE Method for a Nonlinear Distributed-Order Sub-Diffusion Equation

In this article, some high-order time discrete schemes with an H 1 -Galerkin mixed finite element (MFE) method are studied to numerically solve a nonlinear distributed-order sub-diffusion model. Among the considered techniques, the interpolation approximation combined with second-order σ schemes in time is used to approximate the distributed order derivative. The stability and convergence of the scheme are discussed. Some numerical examples are provided to indicate the feasibility and efficiency of our schemes.


Introduction
Fractional differential equations (FDEs) have been greatly developed, and in fact they can be found in many scientific and engineering fields, and include some practical models, such as fractional diffusion or advection models [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17], fractional water wave models [18,19], fractional telegraph models [20], fractional Cable models [21,22], and fractional control models [23]. It is very important to find the solution of FDEs. However, the analytical solutions for a large number of FDEs are hard to look for, so we have to develop efficient numerical algorithms to find numerical solutions for these model equations. Recently, distributed-order DEs, which are yielded by FDEs, are considered by more and more people.
In [24], Mainardi et al. developed a continuous distribution of fractional time derivatives. In 2009, Diethelm and Ford in [25] presented some numerical analysis for distributed-order differential equations and gave a convergence theory. In [26], Ford et al. looked for the numerical solution of a distributed-order diffusion model using an implicit finite difference approximation. In [27], Ye et al. solved a distributed-order time-fractional diffusion-wave equation by a compact difference algorithm. In [28], Meerschaert et al. studied distributed-order fractional diffusions. In [29], Luchko discussed some boundary value problems for the distributed-order time-fractional diffusion problem. Gao et al. in [30] studied the time multi-term and distributed-order fractional sub-diffusion equations. In [31], Lorenzo and Hartley studied variable order and distributed-order fractional operators. Li et al. [32] developed a novel finite volume method for a space distributed-order diffusion model. Recently, based on the approximation methods [30,33], Li et al. [34] considered the H 1 -Galerkin MFE schemes to a linear distributed-order fractional diffusion problem. In [35], Abbaszadeh and Dehghan solved a 2D distributed-order diffusion-wave model with a time-fractional derivative by using an improved meshless method.
In this paper, we will consider the following nonlinear distributed-order sub-diffusion equation where Ω = (0, 1), J = (0, T], and the nonlinear reaction term m(p) satisfies |m(p)| ≤ C|p| with |m (p)| ≤ C, where C is a positive constant. Define where ω(α) ≥ 0, 1 0 ω(α)dα = c 0 > 0. Here, we develop the H 1 -Galerkin mixed finite element (GMFE) method [18,[36][37][38][39][40][41][42] to solve the initial-boundary problem (1)-(3) of the nonlinear time distributed-order model. Compared with the standard finite element method, the H 1 -GMFE method can simultaneously obtain the numerical solution for both unknown scalar function p and its derivative u = p x . Here, motivated by the works [21,22,30], some second-order σ schemes combined with the higher accurate interpolation approximation were developed to discretize the distributed-order derivative. Meanwhile, the temporal integer derivative ∂u ∂t at any point t n−1+σ (∀σ ∈ [ 1 2 , 1]) is approximated by some second-order σ formulas which are proposed in [21]. The second-order σ formulas are called second-order θ schemes in [21]. Some terms including the first derivative and the nonlinear term can be discretized by the corresponding second-order approximation formulas, which can be found in Lemmas 8-10 in this article. Here, the σ can be selected according to the method in [30]; see also Lemmas 2-3 in our current article. The stability and optimal a priori error estimates for both p and u are derived, and some numerical tests on three parameters are implemented to verify our theoretical results.
The article is organized as follows. In Section 2, we give some lemmas. In Section 3, we analyze the stability for the H 1 -Galerkin mixed finite element scheme. In Section 4, we demonstrate the error analysis in L 2 -norm. In Section 5, some numerical experiments are presented to confirm the theoretical analysis. Throughout the paper, we denote C > 0 as a positive constant.
is monotonically decreasing and convergent to σ * .

Remark 1.
It is easy to know from Lemmas 2-3 that σ is how to be selected in the current algorithm. One can also see from Lemma 2 that the value of σ depends on two step size parameters ∆α and ∆t, and in fact this relation can also be verified by the computed data of Tables 1-3 in Section 5.
Based on the chosen σ by using Lemmas 2-3, we will give a numerical differentiation formula at the point t = t n−1+σ = (n − 1 + σ)∆t and reveal its numerical accuracy. [30]) Suppose q ∈ C 3 ([t 0 , t n ]), then we have

Lemma 4. (See
where L 1,n (s) is a linear polynomial on [t n−1 , t n ] satisfying L 1,n (t n−1 ) = q(t n−1 ), L 1,n (t n ) = q(t n ), Then, we can get Based on Lemma 1 and Lemma 5, we easily prove that the following error result holds.

Lemma 5.
The numerical differentiation formula Dq(t n−1+σ ) can be rewritten as When n = 1, When n ≥ 2, According to [30], we can get the following two properties of the coefficients {κ n k }.

Numerical Approximation and Analysis of Stability
By introducing u = ∂p ∂x , the nonlinear distributed-order fractional sub-diffusion Equation (1) can be rewritten as We formulate the weak formulation for system (16) at t = t n−1+σ : Find a pair: {p, u} where Using Lemma 1, we discretize the integral term in Equation (17) where R 1 = O(∆α 2 ) and is defined in Lemma 1.
Based on Lemmas 2-10 and formula (18), Equation (17) can be changed as follows: When n ≥ 2, where Set V h and W h as finite dimensional subspaces of H 1 0 and H 1 , respectively, and have the following approximation properties, where 1 ≤ p ≤ ∞, k and r are positive integers.
where W m,p (Ω), 1 ≤ p ≤ ∞ is the classical Sobolev spaces and we denote it as W m,p with norm · m,p .
When n ≥ 2, Next, we will give the analysis of stability for Schemes (23)- (24).

Theorem 1.
For systems (23)-(24), the following stability holds: From Lemma 11 and [34], we have Substitute formula (27) into Equation (26) to get By multiplying 4∆t on both sides of inequality (28) and summing from 2 to n, n = 2, 3, · · · , Z (Z ≤ N), we can arrive at By using Lemma 11 and Cauchy-Schwarz inequality as well as Young inequality, formula (29) can be transformed into By a simple derivation, the third term on the left-hand side of inequality (30) can be written as Substitute formula (31) into inequality (30) to derive Noting that we have By multiplying 2∆t on both sides of formula (35) and using Cauchy-Schwarz inequality with Young inequality, we arrive at By using 2∆t(κ 1 Thus, we have Noting that 2σ − 1 ≥ 0, we derive the first term on the right-hand side of formula (32) as follows: Substituting formula (39) into formula (32), we can get Notice that the following inequality holds: Combine formula (40) with formula (41) to arrive at In Equation (24) (a), we take v h = p n−1+σ h and use Poincaré inequality to arrive at Substitute formula (43) into formula (42) to obtain Using Gronwall lemma, we have (45) Combining formula (43) with formula (45) and using Lemmas 6-7, we have Using formula (45) and formula (46), we accomplish the proof.

A Priori Error Analysis
Following Reference [37], we introduce elliptic projections {p h ,û h } : [0, T] → V h × W h and get the following results: To simplify the notations, denote Theorem 2. Letting {p(t n ), u(t n )} be the solution of systems (19)- (20) and {p n h , u n h } be the solution of systems (23)- (24), there exist constants C independent of h, ∆t, and α such that (a) u n − u n h j ≤ C(∆α 2 + ∆t 2 + h min(k+1,r+1−j) ), j = 0, 1, Proof. When n ≥ 2, the error equations are as follows: When n = 1, we can obtain the following error equations: Setting v h = ϕ n−1+σ , w h = η n−1+σ in Equation (49) and using formula (47), we can write Equation (49) as Using a similar method to the one in [34] and Lemma 11, we have Substitute formula (52) into Equation (51) (b) to get Making use of a similar derivation as that in [34], we obtain Then, we have Substitute formula (55) into formula (53) to obtain Using a simple derivation, we rewrite the second term on the right-hand side of formula (56) as Substituting formula (57) into formula (56), we have Multiplying both sides of Equation (58) by 4∆t and summing from 2 to n, n = 2, 3, 4, · · · Z (Z ≤ N), we obtain For formula (59), we have From [21], we can get Substituting formula (61) into formula (60) and using Cauchy-Schwarz inequality as well as Young inequality, we can get When n = 1, setting w h = η σ in Equation (50) (b), we arrive at Multiply formula (63) by 2∆t and use Cauchy-Schwarz inequality and Young inequality to get Thus, we have Substituting formula (65) into formula (62) and using Lemma 11, we have By using a similar process as in [21,34], we have Notice from Lemmas 6-7 that we can get Take v h = ϕ n in Equation (49) (a) and use Poincaré inequality to arrive at Substitute formula (70) into formula (69) to get Using formula (70) and Gronwall lemma, we can get Using triangle inequality, we finish the proof of the theorem.

Numerical Example
In this section, we consider a numerical example in space-time domain [0, 1] × [0, 1 2 ] to verify the theoretical results of the H 1 -Galerkin mixed finite element method. We choose ω(α) = Γ(5 − α) in this numerical test. By choosing the nonlinear term m(p) = sin(p), and the exact solution Here, for showing the feasibility and validity of our numerical method, we consider a linear element and provide the computing result by the Matlab procedure. In Table 1, with ∆α = 1 500 , h = 1 500 and changed ∆t = 1 10 , 1 20 , 1 40 , the errors and approximating second-order convergence rates in time are obtained, which are in agreement with the theoretical results. In Table 2, we obtain the distributed-order convergence results under the case ∆t 2 ≈ h 2 ≈ ∆α 4 . In Table 3, we show the errors and approximating second-order spatial convergence rates with ∆t = 1 500 , ∆α = 1 500 and changed h = 1 10 , 1 20 , 1 40 , which tell one that the results are in agreement with the case based on the linear element. Furthermore, the comparison between the numerical solution p h and the exact solution p at time t = 0.2, 0.3, 0.4, 0.5 with ∆t = 1 100 , ∆α = 1 8 and h = 1 15 is made in Figure 1. Similarly, in Figure 2, the comparison between the exact solution u and the numerical solution u h is also shown based on the same discrete parameters as shown in Figure 1.