1 Introduction

Fractional calculus has a history of at least three hundred years. However,it is in the past several decades that the adoption of fractional differential equations (FDEs) for mathematical modeling has attracted great attention and become increasingly popular. Many problems in physics,chemistry,biology,and other fields can be expressed very suitably and successfully by FDEs^[1]. Unfortunately,it is difficult to obtain the analytic solutions of many FDEs,and approximate methods should be considered.

At present,some numerical methods have been developed to obtain the approximate solutions for fractional equations. Finite difference (FD) methods were studied in Refs. [2]-[3],a spectral method approximation was proposed by Lin and Xu^[4],finite element method (FEM) has also been considered to get the approximate solutions of variational problem for the FDEs^{[5, 6, 7, 8, 9]},and so on. Among these numerical methods,the FEM is more popular because of its flexibility to irregular regions and low requirement for the smoothness of solutions. However,the general FEM always includes many degrees of freedom,especially when we solve high dimensional-space problems. It is well-known that fractional differential operators are nonlocal operators. And for this reason,a common problem,which is brought about by applying the general FEM to solve FDEs,is that high consumption of memory and long time calculation inevitably occur in the computational process. To overcome this difficulty,some authors proposed “fixed memory principle”^{[1, 10]} or “nested mesh concept”^[11]. However,these measures are at the cost of certain accuracy with the reduction of memory size. It is worth considering how to choose suitable basis functions for the FEM so that the degree of freedom is decreased, the computing efficiency is improved,and at the same time,the sufficient accurate numerical solutions are guaranteed.

Holmes et al.^[12] introduced a POD method in 1996. This method was called the principalcomponent analysis in statistics^[13] and the Karhunen-Loeve expansion in pattern recognition^[14], which provided us with the means of constructing orthogonal basis to represent the given data by the least squares optimal sense. Some numerical methods combining with the POD theory for partial differential equations can effectively establish reduced-order models,thereby reducing the size of memory and accelerating computational speed. For example,Kunisch and Volkwein^{[15, 16]} presented the method,which combined the Galerkin projection with the POD theory,to solve the parabolic problem and equations in fluid dynamics. Sun et al.^[17],Luo et al.^{[18, 19]},and Luo^[20] have presented reduced-order FD methods,reduced-order finite element (FE) schemes,and reduced-order finite volume element formulations based on the POD as well as error estimates for some equations,all of which have obtained good effects. Until this day, we have not seen any other published papers of addressing to the FEM based on the POD for FDEs except Ref. [21]. In this paper,we introduce the reduced-order FEM based on the POD technique to solve the following time-fractional Tricomi-type equation:

Here,D_α^t u is the αth-order Caputo fractional derivative,and 1 < α < 2,which is defined by

Here,Γ(·) denotes the Gamma function,and γ is a real non-negative constant. Ω ⊂ R² is a bounded and polygonal domain with the boundary ∂Ω. The operator denotes the Laplacian operator with respect to space variables. For Eq. (1),the classical FEM,which is based on the FE formulation for the space and the FD for the time,was considered in Ref. [9]. Its FE scheme shows that,to obtain the solution at the time t_n,the former all solutions on time levels t₁,t₂,· · · ,t_n−1 are needed. Therefore,the memory requirement for computer has to be enlarged. What’s more,the use of traditional basis functions for the FEM leads to solve algebraic equations with a large number of unknown quantities,which requires long time to complete calculation. In this paper,d POD basis functions are constructed and adopted,where d is only the number of the first few maximal eigenvalues of matrix K_Q×Q (see Section 2). On every time level,only d unknown quantities are to be solved. Therefore,the degrees of freedom are reduced,and the computing time is cut down drastically. It should be noted that we construct the POD basis by snapshots which is taken from the solutions of former Q time levels in the general FEM,while for practical problems in the scientific research,we can get the snapshots by drawing samples from experiments or by the past data information.

The plan of this paper is as follows. In Section 2,we first recall the fully-discrete general FE formulation for the time-fractional Tricomi-type equation,then,introduce the POD method and establish the reduced-order FE scheme based on the POD. In Section 3,the stability and error estimate for the reduced-order FE solution are proved particularly. The implementation process of the method is presented in Section 4. Two numerical examples are shown in Section 5 to demonstrate that the numerical results are consistent with the theoretical analysis.

2 Establishment of reduced-order FEM

The Sobolev space used as follows is standard (see Ref. [22]). Let (·,·) represent the L² inner product,and its corresponding L² norm is denoted by || · ||. We denote by || · ||_r the norm in H^r(Ω).

2.1 Recall general FE scheme

We first partition the time interval [0,T] by t_i (i = 0,1,· · · ,N),where 0 = t₀ < t₁ < · · · < t_N = T ,t_n = nk (n = 1,2,· · · ,N),and . Then,for a fixed n,the Caputo fractional derivative (2) is approximated as follows^[7]:

where ,C_u denotes a constant which is relevant to u,and M1 is the upper bound of . When j = 1,u(t−1) can be approximated by u(0) − ku_t(0) = φ₀ − kφ₁,where t₋₁ = −_k.

We can verify easily that b_j have the following properties:

Let be a quasi-uniform triangle subdivision to the space region Ω. The maximal length of the sides of the triangulations is denoted by h. Let S^h ⊂ H₀ ¹(Ω) be a kind of FE space,and its degree of accuracy is not lower than one,i.e.,

Here,P_r is the set of all polynomials whose degree is up to r.

Let u_h⁻¹ and u_hⁿ denote the approximations to u(x,y,t₋₁) and u(x,y,t_n),respectively. Set f ⁿ = f(x,y,t_n). Then,the fully-discrete formulation of (1) can be stated as following.

To get the approximate solution u_hⁿ ∈ S^h (1 ≤ n ≤ N),

holds for any v_h ∈ S^h,where γ₀ = (nk)^2γ,and α₀ = k^αΓ(3 − α). u0 h may be taken as P _φ0, where P is the L₂ projection operator that satisfies (P u − u,v) = 0,∀v ∈ Pr(Ω),and u_h⁻¹ can be taken as P(φ₀ − kφ₁),or else.

The stability of problem (4) is proved in Ref. [9],and the following error estimate is presented.

Lemma 1^[9] Let uⁿ = u(x,y,t_n) be the exact solution of (1) and u_hⁿ be the finite element approximate solution of the fully-discrete formulation (4). Assume that uⁿ satisfy u_tⁿ ∈ L²(0,T ; H^r+1(Ω)) ∩ L∞(0,T ; H_r+1(Ω)),u_ttn ∈ L²(0,T ; L₂(Ω)),with φ₀,φ₁ ∈ H^r+1(Ω). Then, u_hⁿ satisfies the following error estimates:

2.2 POD method

Suppose that we have got the snapshots U_i (i = 1,2,· · · ,Q),at least one of which is assumed to be non-zero. Let

Here,ν is called the space generated by snapshots {U_i}_i^Q =1. Suppose that q = dim ν and {η_j}_j^q =1 is an orthogonal basis of ν. Then,every member of the snapshots can be expressed by

where (U_i,η_j)_X = (▽U_i,▽η_j).

Definition 1^[18] The POD method is to find the orthogonal basis η_k (k = 1,2,· · · ,q) which satisfies that for each d (1 ≤ d ≤ q),the mean square error is minimum on the average, that is,

subject to where ||U_i||_X = ||▽U_i||. A solution set {η_k}_k=1^d of (6) and (7) is called a POD basis of rank d.

Introduce a matrix K = (K_ij)_Q×Q associated with the snapshots {U_i}_i=1^Q by

Since the matrix K is positively semi-definite and has the rank q,the solutions of (6) and (7) exist,and the following results can be obtained.

Lemma 2^{[15, 16]} Let λ₁ ≥ λ₂ ≥ · · · ≥ λ_q > 0 denote the positive eigenvalues of matrix K,and ω¹,ω²,· · · ,ω^q are the associated orthogonal eigenvectors. Then,a POD basis of rank d ≤ q is expressed by

where (ωⁱ)_j is the jth component of eigenvector ωⁱ. Furthermore,it holds the following error estimate:

where λ_l depends on the steps h and k since ν relies on them.

2.3 Reduced-order FE scheme

We first take u^h_i (i = 1,2,· · · ,Q,Q << N),which are the solutions of (4),as the snapshots for convenience. Let S^d = span{η₁,η₂,· · · ,η_d} (1 ≤ d ≤ q). Define the Ritz-projection P ^d : S^h → S^d denoted by

From the linear operator theory,there exists a unique extension P ^h : H₀ ¹ → S^h of P ^d such that P ^h|_S^h = P ^d : S^h → S^d defined by

where u ∈ H¹ ₀(Ω). If u ∈ H ^s(Ω),there holds the following estimates (see Ref. [23]):

and

Lemma 3^[18] For each d (1 ≤ d ≤ q),it holds the following estimates for projection operator P ^d :

where u_hⁱ ∈ ν is the solution of (4).

Now the reduced-order FEM based on the POD theory for Eq. (1) can be stated as follows.

Find u_dⁿ ∈ S^d (n = 1,2,· · · ,N) such that for any v_d ∈ S^d,

Here,u_hⁿ (n = 1,2,· · · ,Q) are chosen from the first Q solutions of (4).

Since (u_dⁿ ,v_d) + α₀(▽u_dⁿ ,▽v_d) is positive definite,the problem (13) has a unique solution u_dⁿ ∈ S^d.

Remark 1 We construct the POD basis by snapshots which are taken from the solutions of former Q time levels in the general FEM. However,the snapshots can be obtained by drawing samples from experiments or by the past data information for practical problems. Suppose that we take S^h as the piecewise linear functions space,that is r = 1 in S^h. If we use the FE scheme (13) to solve the approximate solution of Eq. (1),the number of unknown quantities equals to the amount of inner vertices of triangles in (see Ref. [23]). While if we solve the equation by reduced-order FE scheme (13),the amount of total degrees of freedom is only d (d ≤ q ≤ Q << N ),where d is just the number of former bigger positive eigenvalues of matrix K_Q×Q. Therefore,it is tiny (see the numerical example in Section 5). On the other hand,we have found from the discrete formulations (4) and (13) that the solutions of all former time levels which satisfy t < t_n must be saved and added up in order to get the solution at the time t_n. Consequently,when we solve the numerical solution of Eq. (1) by applying reduced-order FE scheme (13),its computational quantities are largely decreased since the degrees of freedom in each time level are much less than those by using the classic FE formulation.

3 Stability analysis and error estimate for reduced-order solution

In this section,we develop stability analysis as well as the error estimate for the problem (13). From here,we will denote a generic constant by C,which may not be the same at different positions,but it is always independent of the steps k and h. For convenience and without lose of generality,we assume that φ₀(x,y) and φ₁(x,y) are two zero functions in the following stability analysis.

Theorem 1 The reduced-order FE scheme (13) has a unique solution u_dⁿ ∈ S^d which satisfies the following stability result:

Let u_hⁿ and u_dⁿ be the solutions of problem (4) and problem (13),respectively. If k = O(h),and Q = O(N ^1/2),it holds the error estimate

Proof Take vd = u_dⁿ in (13),we can get

Therefore,

In the following,we apply the mathematical induction method to prove the result (14). In fact,when n = 1,it is obvious that

Assuming that (14) holds for n = 2,3,· · · ,m (m < N),we want to verify that it is also true when n = m + 1. Considering (3),it holds

Therefore,

Then,(14) is proved.

In order to carry out the error estimate,we subtract (13) from (4) and get the error equation

According to (8) and (15),we get

By the Hölder inequality and (10),we obtain that when k = O(h),

Therefore,

Setting Q = O(N ^1/2) and using triangle inequality,(16),(11),and (12),we have

Once again,we can derive the following estimate by the mathematical induction method:

Since it is similar to the above proof,we omit the process.

According to the triangle inequality and the estimate results in Lemma 1 and Theorem 1, the error between the exact solution and reduced-order solution can be easily obtained.

Theorem 2 Let uⁿ = u(t_n) and u_dⁿ be the solutions of the problem (1) and (13),respectively. Assume uⁿ satisfying u_tⁿ ∈ L²(0,T ; H^r+1(Ω)) ∩ L∞(0,T ; H^r+1(Ω)) and u_ttⁿ ∈ L²(0,T ; L₂(Ω)). If k = O(h),Q = O(N 1/2),there holds the following error estimate for n = 1,2,· · · ,Q,Q + 1,· · · ,N:

Remark 2 The condition Q = O(N ^1/2) in Theorem 1 and Theorem 2 shows that only the solutions on the former Q time levels of the general FE scheme need to be taken as snapshots. The condition k = O(h) in theorems is just an assumption for the theoretical analysis. It can be changed properly on other cases. The error estimate results in two theorems also provide a measurement for deciding the number of POD basis,that is,the number d of POD basis should satisfy

4 Implementation of algorithm

In this section,we give the implementation process for solving the problem (13),which can be divided into six steps.

Step 1 Solve the equation (4) for given f(x,y,t) when n = 1,2,· · · ,Q and obtain the snapshots U_i(x,y) = u_hⁱ (i = 1,2,· · · ,Q).

Step 2 Compute the matrix K = (K_ij)_Q×Q,where .

Step 3 Solve the rank q and eigenvalues of matrix K and obtain eigenvalues λ₁ ≥ λ₂ ≥ · · · ≥ λ_q > 0. The associated orthogonal eigenvectors ω¹,ω²,· · · ,ω^q also need to be solved out.

Step 4 Choose the number d ≤ q such that the following inequality is satisfied

Then,the POD basis is generated,i.e.,

Step 5 Set S^d = span{η₁(x,y),η₂(x,y),· · · ,η_d(x,y)},solve Eq. (13) with only d degrees of freedom,and get its solution u_dⁿ (n = 1,2,· · · ,N).

Step 6 For given error ,u_dⁿ is the solution of the problem (13),else let U_i = u_dⁱ (i = n − Q,n − Q + 1,· · · ,n − 1),and repeat Step 1 to Step 5.

5 Numerical examples

In this section,we provide two numerical examples to verify our theoretical analysis results further. Assume that is a triangular partition to the region [0,1] × [0,1],and S^h is the space of piecewise linear and continuous functions,which means that we all take r = 1 in the following examples. Let and .

Example 1 Consider the following problem in the two-dimensional space:

with sin(2πx) sin(2πy) + 8π²t³ sin(2πx) sin(2πy). The exact solution of this equation is u(x,y,t) = t² sin(2πx) sin(2πy).

The error estimate results when T = 1 are listed in Table 1 for different α,M,N,and Q,but d = 2 is only needed. Under the condition k = O(h) and Q = O(N ^1/2),we find that the convergent accuracy in the L² norm for the reduced-order FE solution solved by the formulation (13) almost arrives at second-order for different values of α. The influence of space is obvious since our problem is two-dimensional in space. Although we only adopt d = 2 POD basis functions,high accuracy is still obtained,and meanwhile,the calculation quantities and computer memory are both saved.

Next,we choose another example to show the efficiency of the reduced-order FEM. To this end,we solve the following equation by formulas (4) and (13),respectively.

Example 2

We take the source term

Then,the exact solution of equation (17) is

For α = 1.5,we obtain the approximate solution when T = 2. The FE solution u_h and reduced-order FE solution ud are displayed in Fig. 1,and their error maps with the exact solution are exhibited in Fig. 2. In these computational process,we take M = 64,N = 200,Q = 15, which satisfies Q = O(N ^1/2),k = O(h). d is set to be 6 since it should satisfy d ≤ q and . We can see from the figures that the reducedorder FE solution is still convergent without renewing the POD basis and has the similar effect with the FE solution. In computing time,the classical FE formulation (4),which includes 64 × 64 degrees of freedom on every time level,the required time of solving system of equations is 8 min,while for the reduced-order FEM (13) with 6 POD basis functions,the corresponding time is only 2 s. That means,to solve Eq. (17),the required time of solving equations by the classical FE formulation (4) is 240 times that of the reduced-order FEM based on the POD with only 6 degrees of freedom on every time level. Although the construction of POD basis functions also consumes time,it is much less than that computed by the standard FEM in the whole process since Q = O(N ^1/2).