A Petrov-Galerkin Finite Element Method for Solving the Time-fractional Di ﬀ usion Equation with Interface

. Abstract Time-fractional partial di ﬀ erential equation is widely applied in a variety of disciplines, its numerical solution has at-tracted much attention from researchers in recent years. Time-fractional di ﬀ erential equations with interfaces is a more challenging problem because the governing equation has discontinuous coe ﬃ cients at interfaces and sometimes singular source term exists. In this paper, we propose a Petrov-Galerkin ﬁnite element method for solving the two-dimensional time-fractional di ﬀ usion equation with interfaces. In this method, a ﬁnite di ﬀ erence scheme is employed in time and a Petrov-Galerkin ﬁnite element method is employed in space. Extensive numerical experiments show that for a fractional di ﬀ usion equation of order α with interfaces, our method gets to (2 − α )-order accurate in the L 2 and L ∞ norm.


Introduction
Fractional differential equations are generalized differential operations from integer orders to fractional derivative operations.The problem of fractional calculus has wide applications in fractal theory, signal processing, system control, quantum mechanics, environmental science, and finance (Jiang and Qi, 2010;Mainardi, 1996).The time-fractional differential equation is obtained by replacing the integer order in the classical model to fractional derivative, the fractional derivative at a certain time depends on all the values of the function before this time point, so the fractional partial differential equation is applicable for problems with memory process, genetic property and heterogeneous material (Rossikhin and Shitikova, 1997;Ichise et al., 1971).Laplace and Fourier transform are employed to get the analytical solutions of some fractional partial differential equations, for example, equations with constant coefficients (Baeumer et al., 2005;Duan, 2005;Meerschaert et al., 2002).But for most fractional partial differential equations, analytical solutions are not available, in recent decades, the numerical solution of fractional differential equations has attracted much attention from researchers, numerical methods that are employed include finite difference method (Langlands and Henry, 2005;Liu et al., 2004;Meerschaert and Tadjeran, 2006;Su et al., 2009;Murio, 2008), finite element method (Ervin et al., 2007;Deng, 2008;Burrage et al., 2012;Jiang and Ma, 2011), discontinuous Galerkin method (Mustapha and McLean, 2011;Ji and Tang, 2012), spectral method (Zayernouri and Karniadakis, 2014;Zheng and Wei, 2010), meshless method (Liu et al., 2011;Dehghan et al., 2015), and some other methods (Ray, 2007;El-Wakil et al., 2009;Rihan, 2010).
In (Liu et al., 2002(Liu et al., , 2003;;Huang and Liu, 2005), Liu et al. derived the numerical methods for solving the fractional advection-dispersion equation.In (Liu et al., 2007), the stability and convergence of the difference methods are studied.Meerschaert et al. (Meerschaert et al., 2006) used a variation of the Grünwald finite difference approximation to solve two-dimensional fractional partial differential equations with variable coefficients.Roop (Roop, 2006) pointed out that the main applying the finite element method to the numerical solution of fractional partial differential equations is that the fractional operater is non-local, the resulting matrix is almost dense, so it needs higher computational costs.Zeng et al. (Zeng et al., 2015) provided a second-order accurate scheme for the time-fractional subdiffusion equation.The authors gave strict analysis of the stability and convergence of this method.Compared with the finite difference methods and finite element methods, the meshless method has the advantage that they use higher order trial functions, so that they usually have higher accuracy.(Liu et al., 2011) used an implicit meshless radial basis function(RBF) to solve the time fractional diffusion equations.Dehghan et al. (Dehghan et al., 2015) used the meshless method to solve the time fractional sine-Gordon and Klein-Gordon equations.In (Li and Xu, 2009), a spectral method is proposed to solve the time fractional diffusion equation.The authors formulated the time fractional fractional differential problem into an elliptic problem by finding proper spaces and norms.Zheng et al. (Zheng et al., 2015) presented a space-time spectral method for the time fractional Fokker-Planck equation with high order accuracy and efficiency.
In order to simulate problems with irregular convex domains, Fan et al. (Fan et al., 2017) proposed an unstructured mesh finite element method to solve the time-space fractional wave equation on an irregular convex domain.For many practical problems, interfaces are involved in the research domain.Interface problems attract much interest in recent years because of their wide applications in fluid dynamics, electromagnetic, biological systems (Chadam and Yin, 1993;Kandilarov, 2005), etc.The pioneering work was done by Peskin in the simulation of blood flow in heart.He proposed the "immersed boundary" method (Peskin, 1977;Peskin and Printz, 1993).In (Sussman et al., 1994), the "immersed boundary" method was combined with the level set method, resulting in a first order numerical method that is simple to implement, even in multiple spatial dimensions.The "immersed interface" method presented in (LeVeque and Li, 1994) can get second-order accuracy.(Chen and Zou, 1998) used finite element method to solve elliptic problems with the interface conditions [u] = 0 and [βu n ] 0, second order accuracy in the energy norm and nearly second order accuracy in the L 2 norm can be achieved.The boundary condition capturing method (Liu et al., 2000) uses the Ghost fluid method (GFM) (Fedkiw et al., 1999) to capture the boundary conditions.This method was sped up by a multi-grid method (Justin and Liu, 2004), and in (Macklin and Lowengrub, 2008), the method is improved to second order accuracy for smooth interfaces.In (Zhou et al., 2006), the matched interface and boundary (MIB) method was proposed to solve elliptic equations with smooth interfaces.In (Yu et al., 2007), the MIB method was generalized to treat sharp-edged interfaces.(Colella and Johansen, 1998;Oevermann and Klein, 2006;Oevermann et al., 2009) developed high order methods for elliptic equations in complex domains from the finite volume perspective.(Ying and Henriquez, 2007) is a kernel-free boundary integral (KFBI) method for solving elliptic BVPs.In (Hou and Liu, 2005), a finite element method with non-body-fitting grids is proposed to solve elliptic equations with matrix coefficients and sharp-edged interfaces.In (Wang et al., 2017(Wang et al., , 2018)), the method is improved to solve three-dimensional problems.
Although the numerical solution of fractional differential equations and interface problems have gained much attention respectively, there are few research about the fractional differential equations with interfaces.The main challenge of solving this kind of problem is that the solution at current time level, the coefficient matrix and the jump conditions all depend on the solutions at all previous time levels because of the existence of interfaces and fractional derivative in time.In this paper, we propose a Petrov-Galerkin finite element method to solve the time-fractional diffusion equation with interfaces.The main idea is to use a finite difference scheme in time and a Petrov-Galerkin finite element method with non-body-fitting grids in space.Extensive numerical experiments demonstrate the effectiveness of our method.
|∇ϕ| can be obtained on Ω, which is a unit normal vector of Γ pointing from Ω − to Ω + .
Let T > 0, we seek solutions of the variable coefficient time-fractional diffusion equation away from the interface Γ given by in which α is the order of the time-fractional derivative, x = (x, y) denotes the spatial variables and ▽ is the gradient operator.The coefficient β(x, t) is assumed to be a 2 × 2 matrix that is uniformly elliptic on each disjoint subdomain, Ω − and Ω + , and its components are continuously differentiable on each disjoint subdomain, but they may be discontinuous across the interface Γ.The right-hand side f (x, t) is assumed to lie in L 2 (Ω).
Given functions a and b along the interface Γ, we prescribe the jump conditions The " ± " superscripts refer to limits taken from within the subdomains Ω ± .For ease of discussion in this section, and for accuracy test in the numerical examples, we assume that a and b are smooth on the closure of Ω.
The boundary condition is given as for a given function g on the boundary ∂Ω.
Finally, we prescribe the initial condition as The main challenge of solving this time-fractional diffusion equation with discontinuous coefficients is how to deal with the fractional derivatives is the Caputo fractional derivatives of order α.
Specifically, when α = 0, equation ( 1) turns out to be a model of elliptic equation with discontinuous coefficients on the interface, please refer to (Hou et al., 2010) for more details.
When α = 1, equation ( 1) turns out to be a classical model of diffusion equation with discontinuous coefficients on the interface, (Wang and Shi, 2015) for more details.
In this paper, the value of α under consideration is 0 < α < 1.In (Lin and Xu, 2007), the Caputo fractional derivative is defined as From this definition, it is obvious that the Caputo fractional derivatives uses the information of the standard derivatives at all previous time steps.Let ∆t = T K be the time step, then For simplicity, we use u k to represent u(t k ).Putting equations ( 1)-( 5) together, equation ( 1) can be rewritten to the following full discrete equation where (Lin and Xu, 2007), it is a (2 − α)-order accuracy scheme.

Numerical Method
In this section, we are going to introduce our method for solving the two-dimensional time-fractional diffusion equation with interfaces using the Petrov-Galerkin finite element method in space and the finite difference method in time.By the definition of the fractional derivative and from equation ( 7), we can see that the solution at time t k depends on the solutions at all previous time steps.Since the problem under consideration is the time-fractional diffusion equation with interfaces, the coefficients β(x, t) and the jump conditions a(x, t) and b(x, t) in this paper all depend on time.In this section, we are going to construct and solve the local system.
For ease of discussion, we restrict ourselves to a rectangular domain Ω = (x min , x max ) × (y min , y max ) in the plane.Given positive integers I and J, set ∆x = (x max − x min )/I and ∆y = (y max − y min )/J.Define a uniform Cartesian grid (x i , y j ) = (x min + i∆x, y min + j∆y) for i = 0, ..., I and j = 0, ..., J.Each (x i , y j ) is called a grid point.For the case i = 0, I or j = 0, J, a grid point is called a boundary point, otherwise it is called an interior point.The grid size is defined as h = max{∆x, ∆y} > 0.
Two sets of grid functions are needed and they are denoted by Cut each rectangular region [x i , x i+1 ]×[y j , y j+1 ] into two pieces of right triangular regions: one is bounded by x = x i , y = y j and y = Collecting all these triangular regions, we obtain a uniform triangulation T h : If ϕ(x i , y j ) ≤ 0, we count the grid point (x i , y j ) as in Ω − ; otherwise we count it as in Ω + .A cell K is called a regular cell if all its vertices belong to the same subdomain, either Ω + or Ω − , otherwise it is called an interface cell.There are four kinds of interface cells, see Fig. 2. In an interface cell, we write K = K + ∪ K − .K + and K − are separated by a straight line segment, denoted by Γ h K .The two end points of the line segment Γ h K are located on interface Γ and their locations can be calculated from the linear interpolations of the discrete level-set functions ϕ h = ϕ(x i , y j ).
) is a standard continuous piecewise linear function, which is a linear function in every triangular cell and T h (ψ h ) matches ψ h on grid points.Clearly such a function set, denoted by H 1,h 0 , is a finite dimensional subspace of H 1 0 (Ω).The second extension operator U h is constructed as follows: for any u h,k ∈ H 1,h ± with u h,k = g h,k at boundary points, U h (u h,k ) is a piecewise linear function and matches u h,k on grid points.It is a linear function in each regular cell, just like the first extension operator U h (u h,k ) = T h (u h,k ) in a regular cell.In each interface cell, it consists of two pieces of linear functions, one is on K + and the other is on K − .The location of its discontinuity in the interface cell is the straight line segment Γ h K .Note that two end points of the line segment are located on interface Γ, hence the interface condition [u] = a could be and is enforced exactly at these two end points.In each interface cell, the interface condition [β ▽ u • n] = b is enforced with value b at the middle point of Γ h K .We shall construct an approximate solution to the interface problem taking into account the jump conditions.Note that the Neumann jump condition along the interface Γ is already absorbed into the weak formulation, hence, we only need to take care of the Dirichlet jump condition along the interface Γ.We shall seek an approximate solution which is piecewise linear on both subdomains Ω − and Ω + , but discontinuous along the interface Γ. Clearly, in cases (a)-(c) of Fig. 2, when a vertex (x, y) of the interface cell K is on the interface, we need to get two solutions defined at the same point (x, y) and at the same time t k .We will compute one solution u +,k (x, y), and the other solution u −,k (x, y) can be derive from the Dirichlet jump condition, such that u −,k (x, y) = u +,k (x, y) − a k (x, y).To this end, we define u k (x, y) as follows: ) can be constructed uniquely at t k , provided T h , ϕ, a k and b k are given.Proof.Based on the classification of cells, we take the following five different cases into consideration: Case 0. For the regular cell K, K ⊂ Ω + or K ⊂ Ω − , depending on which subdomain all the vertices of K belong to.For this kind of cells, no extra treatment is needed to construct U h (u h,k ).Clearly, U h (u h,k ) is a standard linear finite element function in a regular cell, where the coefficients c ±,k 1 , c ±,k 2 and c ±,k 3 can be calculated by the following system, The superscript + or − here depends on the subdomain Ω + or Ω − the cell K belongs to.
Case 1.For the case when the interface comes across one vertex of the cell K, K ⊂ Ω + or K ⊂ Ω − , U h (u h,k ) can be construct as follows, Before calculating the coefficients c ±,k 1 , c ±,k 2 and c ±,k 3 , we need to check whether an vertex is on the interface or not.If it is on the interface, then it belongs to Ω + , otherwise it depends on which domain it belongs to.From Fig. 2(a) we can see that point p 1 is on the interface.If p 2 and p 3 are in Ω + , then K ⊂ Ω + .The coefficients c +,k 1 , c +,k 2 and c +,k 3 can be calculated by the following system, Otherwise K ⊂ Ω − , we need to take the Dirichlet jump condition into consideration.Then the coefficients c −,k 1 , c −,k 2 and c −,k 3 will be calculated by the following system, Case 2. The interface covers one side of the cell K, K ⊂ Ω + or K ⊂ Ω − .For this kind of cells, U h (u h,k ) can be construct as follows, 2(b) we can see that points p 1 and p 3 are both on the interface.If p 2 belongs to Ω + , then K ⊂ Ω + .The coefficients c +,k 1 , c +,k 2 and c +,k 3 can be calculated by the following system, Otherwise K ⊂ Ω − , we need to take the Dirichlet jump condition into consideration.Then the coefficients c −,k 1 , c −,k 2 and c −,k 3 will be calculated by the following system, In the last two cases, the interface cell K is separated into two different pieces K + and K − , such that K = K + ∪ K − .In these special cases, we need to take the Dirichlet jump condition and the Neumann jump condition into consideration.Special piecewise linear polynomials satisfying interface jump conditions are employed on this cell.The key idea of our method for this kind of interface cells is to construct a piecewise function by two linear polynomials defined on K + and K − .Let p 6 be the middle point of the interface segment Γ h K on this cell.In Fig. 2(c), p 6 is the middle point of p 2 and p 5 , the Dirichlet jump condition is enforced at points p 5 and p 6 , and the Neumann jump condition is enforced at point p 6 .In Fig. 2(d), p 6 be the middle point of p 4 and p 5 , the Dirichlet jump condition is enforced at points p 4 and p 5 , and the Neumann jump condition is enforced at point p 6 .Before continue our discussion, some of the following notations are needed, Case 3. The interface comes across one vertex and cuts one side of the cell K. U h (u h,k ) can be constructed as follows, 2 and c ±,k 3 can be calculated by the following system, Then the coefficients can be calculated by the following system, Case 4. The interface cuts two sides of the cell K, see Fig. 2(d).U h (u h,k ) can be constructed as follows, If p 1 belongs to Ω + , then K + = K 145 ⊂ Ω + and K − = K 4235 ⊂ Ω − .The coefficients c ±,k 1 , c ±,k 2 and c ±,k 3 can be calculated by the following system, Then the coefficients can be calculated by the following system, From Theorem 3.1 we can see that the gradient of u on a cell K can be represented by a linear combination of Clearly, in Cases 0, 1 and 2, the coefficients of a k 4 , a k 5 , a k 6 and b k 6 are 0; in Case 3, the coefficient of a k 4 is 0; in Case 4, and the coefficient of a k 6 is 0. Based on the above discussion, we derive the following Theorem: Theorem 3.2.All coefficients c k in equation ( 8) are finite and independent of u k , a k and b k .
For simplicity, we define some notations.Let After the construction of U h (u h,k ) and with the notations defined in equation ( 9), we further derive the integral equation defined on cell K. Again, we deduce the integral equation for each case mentioned in Theorem 3.1.
In Cases 0 and 1, the integral equation can be derived directly without any special treatment, The sign of the superscript depends on which subdomain the cell K belongs to, Ω + or Ω − .
The difference between Cases 1 and 2 is that when calculating the integration, Case 1 does not need to calculate the line integration, but Case 2 does.This means in Case 2 the Neumann jump condition needs to be taken into consideration.
When K ⊂ Ω + , we have when K ⊂ Ω − , we have In Cases 3 and 4, the cell K is separated into two pieces, K = K + ∪ K − , hence the integration on this cell is conducted in two separated pieces K + and K − , Based on the above discussion, collecting the integration on every different cells (see equations ( 10)-( 13)) together, we propose the following method: Method 3.1.Find a discrete function u h,k ∈ H 1,h ± with u h,k = g h,k on boundary points so that for all ψ h ∈ H 1,h 0 , we have Remark 3.1.In our implementation, the integrals are computed with Gaussian quadrature rules.Since the interface can separate a cell into many different polygon, when calculating the integral, we further cut it into a few triangles.For each triangular cell, the midpoint of each edge is denoted by p i j .In numerical computation, the average of three f (p i j ) in each cell is utilized.

Numerical Experiments
In all numerical experiments below, the level-set function ϕ(x, y), the coefficients β ± (x, y, t) and the solutions All the examples are defined on the domain All errors in solutions are measured in the L ∞ norm and the L 2 norm in the whole domain Ω.

Conclusion
In this paper, we proposed a numerical method for solving the two-dimensional time-fractional diffusion equation with interface.Finite difference method and finite element method are employed to obtain the numerical solution in time and space.This method is capable of solving variable matrix coefficient problems with sharp-edged interfaces.We used several examples to show that our method has (2 − α) order accuracy in both the L 2 and L ∞ norm.
(a) The interface comes across one vertex of the triangle.(b) The interface covers one side of the triangle.(c) The interface comes across one vertex and one side of the triangle.(d) The interface cuts two sides of the triangle.

Figure 4 .
Figure 4. Numerical results for example 2 Example 3.This example is a "star" interface.The level-set function ϕ, the coefficients β ± and the solution u ± are given as follows: Figure 5. Numerical results for example 3

Table 1 .
Numerical results for example 1 using four kinds of grids

Table 3 .
Numerical results for example 3 using four kinds of grids t × n x × n y ∥ u h − u ∥ ∞ Order ∥ u h − u ∥