Numerical solving unsteady space-fractional problems with the square root of an elliptic operator

An unsteady problem is considered for a space-fractional equation in a bounded domain. A first-order evolutionary equation involves the square root of an elliptic operator of second order. Finite element approximation in space is employed. To construct approximation in time, regularized two-level schemes are used. The numerical implementation is based on solving the equation with the square root of the elliptic operator using an auxiliary Cauchy problem for a pseudo-parabolic equation. The scheme of the second-order accuracy in time is based on a regularization of the three-level explicit Adams scheme. More general problems for the equation with convective terms are considered, too. The results of numerical experiments are presented for a model two-dimensional problem.


Introduction
Non-local applied mathematical models based on the use of fractional derivatives in time and space are actively discussed in the literature [1,2]. Many models, which are used in applied physics, biology, hydrology and finance, involve both sub-diffusion (fractional in time) and super-diffusion (fractional in space) operators. Super-diffusion problems are treated as evolutionary problems with a fractional power of an elliptic operator. For example, suppose that in a bounded domain Ω on the set of functions u(x) = 0, x ∈ ∂Ω, there is defined the operator D: Du = − u, x ∈ Ω. We seek the solution of the Cauchy problem for the equation with the square root of an elliptic operator: for a given f (x, t), u 0 (x), x ∈ Ω using the notation f (t) = f (·, t).
To solve problems with the square root of an elliptic operator, we can apply finite volume and finite element methods oriented to using arbitrary domains and irregular computational grids [3,4]. The computational realization is associated with the implementation of the matrix function-vector multiplication. For such problems, different approaches [5] are available. Problems of using Krylov subspace methods with the Lanczos approximation when solving systems of linear equations associated with the fractional elliptic equations are discussed, e.g., in [6]. A comparative analysis of the contour integral method, the extended Krylov subspace method, and the preassigned poles and interpolation nodes method for solving space-fractional reaction-diffusion equations is presented in [7]. The simplest variant is associated with the explicit construction of the solution using the known eigenvalues and eigenfunctions of the elliptic operator with diagonalization of the corresponding matrix [8,9,10]. Unfortunately, all these approaches demonstrates too high computational complexity for multidimensional problems.
We have proposed [11] a computational algorithm for solving an equation with fractional powers of elliptic operators on the basis of a transition to a pseudo-parabolic equation. For the auxiliary Cauchy problem, the standard twolevel schemes are applied. The computational algorithm is simple for practical use, robust, and applicable to solving a wide class of problems. A small number of pseudo-time steps is required to reach a steady solution. This computational algorithm for solving equations with fractional powers of operators is promising when considering transient problems.
To solve numerically evolutionary equations of first order, as a rule, two-level difference schemes are used for approximation in time. Investigation of stability for such schemes in the corresponding finite-dimensional (after discretization in space) spaces is based on the general theory of operator-difference schemes [12,13]. In particular, the backward Euler scheme and Crank-Nicolson scheme are unconditionally stable for a non-negative operator. As for one-dimensional problems for the space-fractional diffusion equation, an analysis of stability and convergence for this equation was conducted in [14] using finite element approximation in space. A similar study for the Crank-Nicolson scheme was conducted earlier in [15] using finite difference approximations in space. We separately highlight the works [16,17,18], where the numerical methods for solving onedimensional in spaces transient problems for a convection and space-fractional diffusion equation are considered.
In the present paper, we construct unconditionally stable two-and three-level schemes for the approximate solution of unsteady problems with the square root of an elliptic operator. A transition to a new temporary level involves solving the standard elliptic problem for the desired solution. The main computational costs are associated with the evaluation of the right-hand side containing the square root of an elliptic operator.
The paper is organized as follows. The formulation of an unsteady prob-lem containing the square root of an elliptic operator is given in Section 2.
Finite-element approximation in space is discussed in Section 3. In Section 4, we construct a regularized difference scheme and investigate its stability. Regularized schemes of the second-order accuracy are developed in section 5. More general problems of convection-diffusion are discussed in section 6. The results of numerical experiments are described in Section 7. The computational algorithm for solving the equation with a fractional power of an operator based on the Cauchy problem for a pseudo-parabolic equation is proposed. At the end of the work the main results of our study are summarized.

Problem formulation
In a bounded polygonal domain Ω ⊂ R m , m = 1, 2, 3 with the Lipschitz continuous boundary ∂Ω, we search the solution for the problem with a fractional power of an elliptic operator. Define the elliptic operator as The operator D is defined on the set of functions u(x) that satisfy on the boundary ∂Ω the following conditions: where µ(x) ≥ µ 1 > 0, x ∈ ∂Ω.
In the Hilbert space H = L 2 (Ω), we define the scalar product and norm in the standard way: For the spectral problem we have λ 1 ≤ λ 2 ≤ ..., and the eigenfunctions ϕ k , ϕ k = 1, k = 1, 2, ... form a basis in L 2 (Ω). Therefore, Let the operator A be defined in the following domain: Under these conditions the operator D is self-adjoint and positive defined: where I is the identity operator in H. For δ, we have δ = λ 1 . In applications, the value of λ 1 is unknown (the spectral problem must be solved). Therefore, we assume that δ ≤ λ 1 in (3). Let us assume for the square root of the operator D We seek the solution of the Cauchy problem for the evolutionary first-order equation with the square root of the operator D. The solution u(x, t) satisfies the equation and the initial condition The key issue in the study of the computational algorithm for solving the Cauchy problem (4), (5) is to establish the stability of the approximate solution with respect to small perturbations of the initial data and the right-hand side in various norms.

Discretization in space
To solve numerically the problem (4), (5), we employ finite-element approximations in space [19,20]. For (1) and (2), we define the bilinear form Define a subspace of finite elements V h ⊂ H 1 (Ω). Let x i , i = 1, 2, ..., M h be triangulation points for the domain Ω. Define pyramid function We define the discrete elliptic operator D as The square root of the operator D is defined similarly to D 1/2 . For the spectral problem The domain of definition for the operator D is The operator D acts on a finite dimensional space V h defined on the domain D(D) and, similarly to (3), where δ ≤ λ 1 ≤ λ 1 . For the square root of the operator D, we suppose For the problem (4), (5), we put into the correspondence the operator equation where ψ(t) = P f (t), w 0 = P u 0 with P denoting L 2 -projection onto V h . Now we obtain an elementary a priori estimate for the solution of (2), (3) assuming that the solution of the problem, coefficients of the elliptic operator, the right-hand side and initial conditions are sufficiently smooth.
Let us multiply equation (2) by w and integrate it over the domain Ω: In view of the self-adjointness and positive definiteness of the operator D 1/2 , the right-hand side can be evaluated by the inequality By virtue of this, we have The latter inequality leads us to the desired a priori estimate: We will focus on the estimate (9) for the stability of the solution with respect to the initial data and the right-hand side in constructing discrete analogs of the problem (7), (8).

Two-level scheme
Let τ be a step of a uniform grid in time such that w n = w(t n ), t n = nτ , n = 0, 1, ..., N, N τ = T . To solve numerically the problem (7), (8), we use the simplest implicit two-level scheme. We start with the simplest explicit scheme When considering the standard evolutionary problems, advantages and disadvantages of explicit schemes are well-known (see, e.g., [12,13]). For problems with the square root of the operator, the main drawback (conditional stability) for explicit schemes remains. In addition, the advantage in terms of implementation simplicity is no more. In this case, the approximate solution at a new time level is determined via (10) as Thus, we must calculate D 1/2 w n . Consider now the possibility of using implicit schemes. Let us approximate equation (7) by the backward Euler scheme: The main advantage of the implicit scheme (13) in comparison with (10) is its absolute stability. Let us derive for this scheme the corresponding estimate for stability. Multiplying equation (13) scalarly by τ w n+1 , we obtain The terms on the right side of (14) are estimated using the inequalities: The substitution into (14) leads to the following level-wise estimate: This implies the desired estimate for stability: which is a discrete analog of the estimate (9). To obtain the solution at the new time level, it is necessary to solve the problem In our case, we must calculate the values of A more complicated situation arises in the implementation of the Crank-Nicolson scheme: In this case, we have i.e., we need to evaluate both Φ(z) = (1 + 0.5τ z 1/2 ) −1 and Φ(z) = z 1/2 . The numerical implementation of the above-mentioned approximations in time for the standard parabolic problems in (7) is based on calculating the values of Φ(D)b for Φ(z) = (1 + στ z) −1 , σ = 0.5, 1 and Φ(z) = z. For problems with the square root of elliptic operators, we apply the approach proposed early in our paper [11]. It is based on the computation of Φ(D)b for Φ(z) = z −1/2 .
Now we will derive the stability conditions for the regularized scheme (11), (16) and after that we will select the appropriate regularizing operator R itself. Rewrite equation (16) in the form Multiplying it scalarly by τ (w n+1 + w n ), we get For G = G * > 0, and G = I + O(τ ), we obtain the inequality Thus, for the regularized difference scheme (11), (16) the following estimate for stability with respect to the initial data and the right-hand side holds: To select an appropriate regularizing operator R, we should take into account two conditions, i.e., first, to satisfy the inequality G > 0, and secondly, to simplify calculations. In the scheme (16), we put In view of D − 2D 1/2 + I ≥ 0, in selecting (19), we have for σ ≥ 0.25 the inequality G > 0 holds. The result of our analysis is the following statement.

The scheme of the second-order accuracy
The regularized scheme (11), (16) has the first-order accuracy in time. Let us consider the possibility to construct second-order regularized schemes in the class of three-level schemes.
Among multilevel difference schemes we can select the explicit Adams methods. To start with the generating scheme, we take the explicit three-level scheme with specified w 0 , w 1 . It is possible to define w 1 via the explicit scheme First of all, we formulate a condition for the stability of the scheme. The scheme (20) may be written in the canonical form for three-level difference schemes [12,13]: Taking into account w n+1 − w n = 1 2 (w n+1 − w n−1 ) + 1 2 (w n+1 − 2w n + w n−1 ), 3w n − w n−1 = 1 2 (w n+1 − w n−1 ) − 1 2 (w n+1 − 2w n + w n−1 ) + 2w n , the scheme (20) obtains the representation (21), wherẽ The necessary and sufficient condition for stability of the three-level scheme (20) with self-adjoint operators (see [12,13]) has the form By (22), the explicit scheme (20) is stable under the conditioñ This gives the constraint on a time step: In compare with the two-level explicit scheme, this reduce the permissible time step by a factor of two. The construction of a three-level regularized scheme is conducted similarly to the scheme (11), (16). Remaining in the class of the schemes with the secondorder approximation, instead of (22), we use the regularized scheme Theorem 2. The regularized scheme (23) with the regularizer S = στ 2 D is unconditionally stable for σ ≥ 1/4.

Proof.
The scheme (23) may be written in the canonical form (21) with For S ≥ 0 the condition for stability of (22) holds if Thus, the regularized scheme (23) is unconditionally stable if σ ≥ 1/4 in (19) .
The computational implementation of the scheme (19), (23) is similar to the numerics of a two-level scheme (16), (19).

Convection-diffusion problems
Convection-diffusion problems are typical for mathematical models of fluid mechanics. Heat transfer as well as impurities spreading are occurred not only due to diffusion, but result also from medium motion. Principal features of physical and chemical phenomena observing in fluids and gases [22,23] are generated by media motion resulting from various forces. Computational algorithms for the numerical solution of such problems are of great importance; they are discussed in many publications (see, e.g., [24,25]). Here we discuss the problem of constructing unconditionally stable schemes for convection-diffusion problems, where diffusive transport is described by the square root of an elliptic operator.
The convection is provided by transient field of medium velocity v(x, t). On the boundary the following condition is satisfied v(x, t) = 0, x ∈ ∂Ω.
We take the convection operator in the skew-symmetric form [21,26]: Using the restrictions (24) the convection operator is skew-symmetric: We seek the solution of the Cauchy problem for the convection-diffusion equation. The solution u(x, t) satisfies the equation and the initial condition (5). The main feature of the problem is that the evolutionary first-order equation has the square root of the operator D.
The discrete convection operator is introduced similarly. Taking into account (24), (25) we set We define the operator C as Similarly to (26) we have For the problem (5), (27), we put into the correspondence the operator equation for w(t) ⊂ V h : taking into account the previously introduced notation and the initial condition (8).
Let us represent the two-level regularized scheme (8), (29) in the form (I +R) w n+1 − w n τ +C w n+1 + w n 2 +D 1/2 w n = ψ n+1 , n = 0, 1, ..., N −1. (30) The stability of this scheme is considered in the same way as for the scheme (16). The key moment is the use of the skew-symmetry property of convective transport operator (28). The scheme (30) may be written as Multiplying it scalarly by τ (w n+1 + w n ), in view of (28), we again obtain where operator G is defined according to (17).
At the new time level we solve the problem Therefore, it is necessary to solve the discrete elliptic problem with the operator and compute D −β w n .
The computational algorithm is based on the consideration of the auxiliary Cauchy problem [11].
Computational implementation of regularized schemes is based on the previously proposed approach involving the solution of the auxiliary Cauchy problem for a pseudo-parabolic equation. It makes possible to construct an unconditionally stable scheme, to obtain an appropriate estimate of stability, and to conduct numerical experiments for the unsteady problem of fractal diffusion. Any other known approaches can also be applied. Here we focus on numerical algorithms with the calculation of D −α y, α = 1/2. But it is possible to employ methods based on an approximation of the operator D −α . In this connection, an approach of great interest is to approximate the operator D −α using operator terms (I + γ i D) −1 [27].
Assume that y(s) = δ 1/2 (s(D − δI) + δI) −1/2 y(0), then for the determination of g n , we can put The function y(s) satisfies the evolutionary equation where G = D − δI ≥ 0. Thus, the calculation of D −1/2 w n is based on the solution of the Cauchy problem (31), (32) within the unit interval for the pseudoparabolic equation.
The test problem is constructed using the exact solution of the problem in the unit circle. The computational domain is a quarter of the circle (see Fig. 1).
We consider the case, where the solution depends only on r, and r = (x 2 1 +x 2 2 ) 1/2 . By virtue of this The solution of the spectral problem is well-known (see, e.g., [28,29]). Eigenfunctions are represented as zero-order Bessel functions: whereas eigenvalues λ k = ν 2 k , where ν k , k = 1, 2, ... are roots of the equation The general solution of the homogeneous (f (t) = 0) Cauchy problem for equation (4) is To verify the accuracy of the approximate solution of the time-dependent equation with the square root of an elliptic operator, we use the exact solution The values of the roots ν 1 , ν 3 for different values of the boundary condition µ are given in Predictions were performed on a sequence of refining grids, which are shown in Figs. 5-7. Suppose that the parameter δ in the estimate (6) for µ = 1, 10, 100 is equal to 1. The scheme (35) is used at each time level on sufficiently fine grid with K = 100. An approximate solution is compared with the exact solution We estimated the accuracy of the regularized scheme (11), (16) on various grids in time and space. The weighting parameter σ was chosen to be 0.25. Table 2 demonstrates predictions for the model problem with µ = 10. The dependence of the accuracy on the boundary condition (on the parameter µ) is illustrated by Table 3. The calculations were performed on grid 2.

Conclusion
1. There is considered a nonclassical problem with the initial data, which is described by an evolutionary equation of first order with the square root of an elliptic operator. The multidimensional problem is approximated in space using standard finite-element piecewise-linear approximations. A priori estimate of stability with respect to the initial data and the righthand side is provided. 2. A regularized two-level scheme of the first-order accuracy in time is constructed. Its unconditional stability is shown. The numerical implementation involves the solution of the standard grid elliptic problem at each time level. In addition, in these problems, we need to solve an operator equation for the square root of an elliptic operator. 3. The possibility of constructing schemes of the second-order accuracy in time for evolutionary problems with the square root of an elliptic operator is studied. The improvement of accuracy is achieved via regularization of a three-level scheme. 4. If convective transport is taken into account, then the skew-symmetric operator is added to the equation. A two-level regularized scheme with the skew-symmetric approximation of convective transport is proposed. It demonstrates unconditional stability. 5. The efficiency of the proposed regularized scheme was demonstrated through the numerical solution of a two-dimensional test problem.