High order numerical schemes for solving fractional powers of elliptic operators

TL;DR

This paper develops high-order finite difference schemes for numerically solving fractional powers of elliptic operators, improving accuracy and efficiency in modeling nonlocal diffusion processes.

Contribution

It introduces a fourth-order finite difference scheme with an optimal weight parameter for pseudo-parabolic problems related to fractional elliptic operators.

Findings

The scheme achieves fourth-order accuracy.

Theoretical analysis confirms stability and convergence.

Computational experiments demonstrate improved performance.

Abstract

In many recent applications when new materials and technologies are developed it is important to describe and simulate new nonlinear and nonlocal diffusion transport processes. A general class of such models deals with nonlocal fractional power elliptic operators. In order to solve these problems numerically it is proposed (Petr N. Vabishchevich, Journal of Computational Physics. 2015, Vol. 282, No.1, pp.289--302) to consider equivalent local nonstationary initial value pseudo-parabolic problems. Previously such problems were solved by using the standard implicit backward and symmetrical Euler methods. In this paper we use the one-parameter family of three-level finite difference schemes for solving the initial value problem for the first order nonstationary pseudo-parabolic problem. The fourth-order approximation scheme is developed by selecting the optimal value of the weight parameter. The results of the theoretical analysis are supplemented by results of extensive computational experiments.

Equations

$$\mathcal{A}u=-\mathrm{div}(k(\mathbf{x})\mathrm{grad}u)+c(\mathbf{x})u$$

$$k(\mathbf{x})\frac{\partial u}{\partial n}+\mu (\mathbf{x})u=0,\mathbf{x}\in \partial \Omega ,$$

$$(u,v)={\int }_{\Omega }u(\mathbf{x})v(\mathbf{x})d\mathbf{x},\Vert u\Vert =(u,u{)}^{1/2}.$$

$$\mathcal{A}{\phi }_k={\lambda }_k{\phi }_k,\mathbf{x}\in \Omega ,$$

$$k(\mathbf{x})\frac{\partial {\phi }_k}{\partial n}+\mu (\mathbf{x}){\phi }_k=0,\mathbf{x}\in \partial \Omega .$$

$${\lambda }_1\le {\lambda }_2\le ....$$

$$u=\sum _{k=1}^{\infty }(u,{\phi }_k){\phi }_k.$$

$$D(\mathcal{A})=\big{\{}u:\;u(\bm{x})\in L_{2}(\Omega),\;\;\sum_{k=0}^{\infty}|(u,\varphi_{k})|^{2}\lambda_{k}<\infty\big{\}}.$$

$$\mathcal{A}={\mathcal{A}}^{\ast }>\delta \mathcal{I},\delta >0,$$

$${\mathcal{A}}^{\alpha }u=\sum _{k=0}^{\infty }(u,{\phi }_k){\lambda }_k^{\alpha }{\phi }_k,$$

$${\mathcal{A}}^{\alpha }u=f.$$

$$a(u,v)={\int }_{\Omega }\left(k\mathrm{grad}u\mathrm{grad}v+cuv\right)d\mathbf{x}+{\int }_{\partial \Omega }\mu uvd\mathbf{x}.$$

$$a(u,u)>\delta \Vert u{\Vert }^2.$$

$$\chi_{i}(\bm{x}_{j})=\left\{\begin{array}[]{ll}1,&\mathrm{if~{}}i=j,\\ 0,&\mathrm{if~{}}i\neq j.\end{array}\right.$$

$$v(\mathbf{x})=\sum _{i=i}^{M_h}{v}_i{\chi }_i(\mathbf{x}),$$

$$a(u,v)=(Au,v),\forall u,v\subset V_h.$$

$$A=A^{\ast }>\delta I,\delta >0.$$

$$A^{\alpha }v=\psi ,$$

$$\Vert v\Vert \le {\delta }^{-\alpha }\Vert \psi \Vert .$$

$$w(t)={\delta }^{\alpha }(t(A-\delta I)+\delta I{)}^{-\alpha }w(0).$$

$$w(1)={\delta }^{\alpha }{A}^{-\alpha }w(0)$$

$$B(t)\frac{dw}{dt}+Dw=0,0<t\le 1,$$

$$B=\frac{1}{\alpha }(tD+\delta I),D=A-\delta I.$$

$$w(0)={\delta }^{-\alpha }\psi .$$

$$D=D^{\ast }>0.$$

$$\Vert w(t)\Vert \le \Vert w(0)\Vert .$$

$$w_n=w(t_n),t_n=n\tau ,n=0,1,...,N,N\tau =1.$$

$$(t_{n+\sigma }D+\delta I)\frac{w_{n+1}-w_n}{\tau }+\alpha D{w}_{n+\sigma }=0,n=0,1,...,N-1,$$

$$w_0={\delta }^{-\alpha }\psi .$$

$$t_{n+\sigma }=\sigma t_{n+1}+(1-\sigma )t_n,w_{n+\sigma }=\sigma w_{n+1}+(1-\sigma )w_n.$$

Full text

††journal: Journal of Computational and Applied Mathematics

1 Introduction

In many recent applications the new mathematical models are proposed, which are based on fractional derivative equations in time and space coordinates [2, 3, 4]. Very different applied mathematical models of physics, biology or finance describe a subdiffusion (represented by fractional-in-time derivatives) or superdiffusion (represented by fractional-in-space derivatives) models. The latter problems are often simulated by using fractional power elliptic operators.

Different numerical techniques, such as finite difference, finite volume methods, can be used to approximate problems with fractional power elliptic operators. In this paper we will use the method of finite elements, since this method is well-suited to solve problems in non-regular domains and to use non-uniform adaptive grids [5, 6]. The implementation of such algorithms require to compute the action of a matrix (operator) function on a vector $\varPhi(A)b$, where $A$ is a given matrix (operator) and $b$ is a given vector. For example, in order to compute the solution of the discrete fractional order elliptic problem, we get $\varPhi(z)=z^{-\alpha}$, where $0<\alpha<1$. There exist various approaches how to compute $\varPhi(A)b$ [7].

The most important class of iterative methods for this purpose are Krylov subspace methods. They are used to solve systems of linear equations obtained after approximation of fractional power elliptic problems (see, e.g. [8]). A comparison of different approaches to solve fractional-in-space reaction-diffusion equations is done in [9]. In particular the integral and adaptively preconditioned Lanczos method are analyzed.

The most straightforward algorithm to solve such systems is to construct explicitly eigenvectors and eigenvalues of the given discrete elliptic operator and to diagonalize the matrix $A$ [10, 11, 12]. But we should note that the direct implementation of this approach is very expensive for general elliptic operators in multidimensional domains. It requires the computation of all eigenvectors and eigenvalues of very large matrices.

A general approach to solve fractional power elliptic problems is based on some approximation of the nonlocal operator.

One can adopt a general approach to solve numerically equations involving fractional power of operators by a popular method is to split the task to solve numerically equations involving fractional power into two steps. First the original elliptic operator is approximated and then the fractional power of its discrete variant is taken. Using Dunford-Cauchy formula the elliptic operator is represented as a contour integral in the complex plane. Then applying appropriate quadratures with integration nodes in the complex plane we get a method that involves only inversion of the original elliptic operator. The approximate operator is treated as a sum of resolvents [13, 14], ensuring the exponential convergence of quadrature approximations.

In paper [15] a more promising quadratures algorithm is proposed, when the integration nodes are selected in the real axis. The new method is based on the integral representation of the power operator [16]. In this case the inverse operator of the fractional power elliptic problem is treated as a sum of inverse operators of elliptic operators.

Such a rational approximation is obtained when the fractional power of the operator is approximated by using the Gauss-Jacobi quadrature formulas for the corresponding integral representation. In this case, we have (see [17, 18]) a Pade-type approximation of the power function with a fractional exponent. The optimal rational approximations are investigated in [19, 20].

A separate class of methods approximates the solution of fractional power elliptic problem by some auxiliary problem of high dimension. In [21] it is shown that the solution of the fractional Laplacian problem can be obtained as a solution of the elliptic problem on the semi-infinite cylinder domain. This idea is used to construct numerical algorithms for solving stationary and non-stationary problems with fractional power elliptic operators [22, 23].

In [24], for solving fractional power elliptic problems we have proposed a numerical algorithm on the basis of a transition to a pseudo-parabolic equation, so called Cauchy problem method. The computational algorithm is simple for practical use, robust, and applicable to solving a wide class of problems. We have used this algorithm also for solving the nonstationary problem with fractional power elliptic operators [25].

For the auxiliary Cauchy problem, standard two-level schemes are applied. Depending on the weight parameters the first and second order accuracy of the approximation is obtained. For many applied problems a small number of pseudo-time steps is sufficient to get a good approximation of the solution of the discrete fractional equation. The efficiency of this algorithm is improved in [26], where a special graded grid in pseudo-time is used.

Another possibility to increase the accuracy of approximations is to use high order discrete schemes for solving the auxiliary pseudo-parabolic equation. In this paper we propose and investigate a fourth order three-level scheme.

The paper is organized as follows. In Section 2 a problem for a fractional power of elliptic operator is formulated. In Section 3 the Cauchy problem method is given. The main results are described in Section 4, where unconditionally stable fourth-order three-level scheme is proposed and investigated. Section 4 provides results of computational experiments, they illustrate the theoretical results on the approximation accuracy of fractional power problems. A model two dimensional problem is solved by using different numerical schemes. At the end of the work the main results of our study are summarized.

2 Problem Formulation

In a bounded domain $\Omega\subset{\mathbb{R}}^{d}$, $d=2,3$ with the Lipschitz continuous boundary $\partial\Omega$ we solve the boundary value problem for the fractional power elliptic operator. The following elliptic operator is defined by:

[TABLE]

where $0<k_{1}\leq k({\bm{x}})\leq k_{2}$, $c({\bm{x}})\geq c_{0}>0$. On $\partial\Omega$ the functions $u({\bm{x}})$ satisfy the boundary conditions

[TABLE]

where $\mu({\bm{x}})\geq 0,\ {\bm{x}}\in\partial\Omega$.

In the Hilbert space $H=L_{2}(\Omega)$ we define the scalar product and norm:

[TABLE]

Next we introduce the eigenvalue problem [27] for (1), (2): find $\varphi_{j}\in H$ and $\lambda_{j}\in\mathbb{R}$ so that

[TABLE]

[TABLE]

The eigenvectors are numbered in such a way, that

[TABLE]

This spectral problem has full set of eigenfunctions $\{\varphi_{k},\;\|\varphi_{k}\|=1,\;k=1,2,...\}$ that span the space $L_{2}(\Omega)$:

[TABLE]

We assume that the operator $\mathcal{A}$ is defined on the domain

[TABLE]

Then $\mathcal{A}$ is self-adjoint and coercive

[TABLE]

where $\mathcal{I}$ is the identity operator in $H$. For $\delta$ we have $\delta<\lambda_{1}$. In most applications, the value of $\lambda_{1}$ is unknown and it should be obtained numerically by solving the eigenvalue problem. In our analysis we assume that a reliable positive bound from below is known $\delta<\lambda_{1}$ in (3).

The fractional power of $\mathcal{A}$ is defined by

[TABLE]

where $0<\alpha<1$. Now we define the boundary value problem for the fractional power of $\mathcal{A}$. The solution $u(\bm{x})$ satisfies the equation

[TABLE]

We approximate the problem (4) by using the finite element method [28]. For the elliptic problem (1), (2) the bilinear form is defined by

[TABLE]

Due to (3), we have that

[TABLE]

We consider a standard sub-space of finite elements $V_{h}\subset H^{1}(\Omega)$. Let us consider a triangulation of the domain $\Omega$ into triangles and let $\bm{x}_{i}$, $i=1,2,...,M_{h}$ be the vertexes of these triangles. As a nodal basis we take the functions $\chi_{i}(\bm{x})\subset V_{h}$, $i=1,2,...,M_{h}$:

[TABLE]

Then for $v\in V_{h}$ we have

[TABLE]

where $v_{i}=v(\bm{x}_{i}),\ i=1,2,...,M_{h}$. We define the discrete elliptic operator $A$

[TABLE]

Similar to (3), the following estimates are valid for $A$:

[TABLE]

The corresponding finite element approximation of equation (4) is: find $v\subset V_{h}$

[TABLE]

where $\psi=Pf$ and $P$ is the $L_{2}$ projection on $V_{h}$. In view of (5), for the solution (6) we get the following simple a priori estimate:

[TABLE]

3 Cauchy problem method

For solving numerically problem (6) we use the Cauchy problem method, proposed in [24]. This method is based on the equivalence of (6) to an auxiliary pseudo-time evolutionary problem. Assume that

[TABLE]

Therefore

[TABLE]

and then $v=w(1)$ if $w(0)=\delta^{-\alpha}\psi$. The function $w(t)$ satisfies the evolutionary equation

[TABLE]

where

[TABLE]

We supplement (8) with the initial condition

[TABLE]

By (5), we get

[TABLE]

The solution of equation (6) can be defined as the solution of the Cauchy problem (8), (9) at the final pseudo-time moment $t=1$.

For the solution of the problem (8), (9), it is possible to obtain various a priori estimates. Here we restrict only to a simple estimate that is consistent with the estimate (7):

[TABLE]

In order to prove (11), it is sufficient to multiply scalarly equation (8) by $\alpha w+tdw/dt$.

To solve numerically the problem (8), (9), the simple implicit two-level Euler scheme can be used [29]. Let $\tau$ be the step-size of a uniform grid in time:

[TABLE]

Let us approximate equation (8) by the implicit two-level scheme

[TABLE]

[TABLE]

We use the notation

[TABLE]

For sufficiently smooth $w(t)$ and $\sigma=0.5$ (the Crank-Nicolson type scheme), the difference scheme (12), (13) approximates the problem (8), (9) with the second order, and with the first order for all other values of $\sigma$.

Theorem 1

For $\sigma\geq 0.5$ the difference scheme (12), (13) is unconditionally stable with respect to the initial data. The approximate solution satisfies the estimate

[TABLE]

Proof 1

In order to prove this result we rewrite equation (12) in the following form:

[TABLE]

Multiplying scalarly it by

[TABLE]

and due to (10), we get

[TABLE]

Since

[TABLE]

then for $\sigma\geq 0.5$ we get the inequalities

[TABLE]

Applying these estimates recursively we prove the validity of (14).

4 Three level schemes

In this section we consider high order schemes. They are based on three level finite difference schemes. For solving problem (8), (9) we use the symmetrical scheme

[TABLE]

with the given initial conditions

[TABLE]

We note that $\overline{w}_{1}$ should be computed by applying some two level numerical algorithm and the accuracy of this approximation should be the same as of the main scheme (15). More details will be given below.

It is well-known that for sufficiently smooth solutions the symmetrical scheme (15) approximates problem (8)–(10) with the second order accuracy.

Next we formulate the stability conditions for the scheme (15), (16). Here we use the general stability results for operator-difference schemes [29, 30].

Let $S$ be a self-adjoint positive operator in $H$. Then we introduce the new Hilbert space $H_{S}$, generated by operator $S$, it consists of elements from $H$ equipped with the energy norm

[TABLE]

Theorem 2

For $\sigma>0.25$ the three-level scheme (15), (16) is unconditionally stable with respect to the initial data. The approximate solution satisfies the estimate

[TABLE]

where

[TABLE]

Proof 2

We rewrite equation (15) in the following form

[TABLE]

let us introduce new functions

[TABLE]

Taking into account that

[TABLE]

we write (19) in the following form

[TABLE]

Multiplying it by

[TABLE]

and taking a discrete inner product, in view that $D$ is self-adjoint, we get

[TABLE]

Since $B_{n}>0$, then we easily get the estimates (17).

Next we consider how to define the initial condition (16) for $w_{1}$. A general approach is to use some two-level solver for $t\in[0,\tau]$. For example, it is possible to apply the symmetrical scheme (12), with $\sigma=0.5$:

[TABLE]

It follows from Theorem 1, that the scheme (21) is unconditionally stable

[TABLE]

and its approximate solution $w_{1}$ converges to $w(t_{1})$ with second order.

It is interesting to see if some explicit schemes can be used to find the initial condition for $w_{1}$. One possibility is to consider the explicit forward Euler scheme

[TABLE]

here the equality $B(0)=\alpha^{-1}\delta I$ is used. In general for sufficiently smooth solutions the $O(\tau^{2})$ accuracy is expected for $w_{1}$. Let us denote the error of the solution of (23) $\tilde{z}_{n}=w_{n}-w(t_{n})$, $n=0,1$. The function $z_{n}$ satisfies the equation

[TABLE]

where $\psi_{1}$ is the standard approximation error. Since $\tilde{z}_{0}=0$, then with sufficiently small step $\tau$ the error $\tilde{z}_{1}$ can be estimated as

[TABLE]

For a sufficiently smooth solution of (8), we have that $\|\psi_{1}\|\leq C\tau$.

More interesting second order explicit schemes can be constructed by using the well-known method described in [29]. The accuracy of the basic forward Euler scheme (23) is increased by using the differential properties of the solution of equation (8)

[TABLE]

We rewrite (24) in the following form

[TABLE]

Then the stability estimate (22) is valid if $\|R\|\leq 1$. For a self-adjoint operator $R$ this estimate is equivalent to the following two-side estimates

[TABLE]

Due to (25) the right inequality of (26) can be written as

[TABLE]

Then the following restrictions on the time step are obtained

[TABLE]

The left inequality can be rewritten as

[TABLE]

For the given values of powers $\alpha$ this inequality is always valid, since

[TABLE]

Due to the obtained stability restrictions (27) the explicit scheme (24) is not recommended for solving real applications.

It is important to note that in the family of second order unconditionally stable three-level schemes (15), (16) it is possible to find such a value of the parameter $\sigma$ which leads to the high order accuracy scheme.

Using the Taylor expansions we get the relations

[TABLE]

[TABLE]

Then the residual of the scheme can be written as

[TABLE]

Taking the solution of (8) we get

[TABLE]

Differentiation of the equation (8) leads to the equality

[TABLE]

Differentiating once more and taking into account linearity of $B$ with respect to $t$ we obtain

[TABLE]

Thus the third order derivative of the solution can be written as

[TABLE]

Substituting this relation into (28) we get the equation

[TABLE]

We approximate the second order derivative in (29) by the standard central difference formula

[TABLE]

Then from (29) we get the semi-difference scheme of the fourth approximation order

[TABLE]

where the optimal weight parameter $\sigma_{0}$ is given by

[TABLE]

Theorem 3

The three-level high order difference scheme (30), (31) is unconditionally stable with respect to the initial data.

Proof 3

The scheme (30), (31) belongs to the family of three-level weighted schemes (15), (16). Thus it is stable if $\sigma_{0}>0.25$. This condition is satisfied for all $0<\alpha<4$, including the considered fractional powers $0<\alpha<1$.

The implementation of the high order three-level difference scheme (30), (31) requires to specify the second initial condition $w_{1}$. It should be computed with the same fourth order accuracy. We do not have any robust, unconditionally stable and efficient two-level high-order difference scheme. In all computations presented in the next section the initial condition $w_{1}$ is computed by using the symmetrical two-level scheme (12) and a sufficiently fine grid is constructed on the time interval $[0,\tau]$.

Let this interval be divided into $m$ sub-intervals. The approximate solutions $w_{\beta/m}$ are computed for time moments $t_{\beta/m}=\beta\tau/m$, $\beta=1,\ldots,m$ by using the following scheme

[TABLE]

For sufficiently smooth solutions of the differential problem (8) and if $m\sim N$ then the solution $w_{1}$ is computed with the required accuracy $\mathcal{O}(\tau^{4})$. We note that the computational complexity of the three-level algorithm is increased approximately twice if such approach is applied to compute the initial condition $w_{1}$.

5 Numerical Experiments

Here we present results of the numerical solution of a model problem (1), (2), (3) in two spatial dimensions, where the computational domain is a unit square

[TABLE]

The coefficients of the operator $\mathcal{A}$ and the right-hand side function $f$ in the equation (4) are defined as

[TABLE]

The piecewise linear continuous $P_{1}$ Lagrange elements are used to approximate the elliptic operator. The domain $\Omega$ is covered by the uniform grid with $50$ intervals in each direction.

The accuracy of different approximations in time will be estimated by a reference solution. It was obtained using the symmetrical two-level scheme (12) with $\sigma=0.5$ and taking a sufficiently small time step: $N=5000$. The relative errors of the approximate solution in the norm $\|\cdot\|$ of space $L_{2}(\Omega)$ and in the norm $\|\cdot\|_{1}$ $H^{1}(\Omega)$ are defined by

[TABLE]

where $\bar{w}_{N}$ is the reference solution. In Fig. 1 we show the reference solution for various values of the fractional power parameter $\alpha$.

For the two-level weighted difference scheme (12) the errors of the solution are presented in Figs. 2 – 4. As it follows from the theoretical analysis the accuracy of approximation is essentially increased for the values of parameter $\sigma$ in the neighbourhood of 0.5.

We also note that the accuracy of the approximate solution is better for larger values of $\alpha$. This result is explained by the increased smoothness of the solution for larger values of $\alpha$.

The main goal of this paper is to investigate the accuracy of the three-level difference scheme (15), (16). First we have used the two-level symmetrical scheme (21) to compute the initial condition for $w_{1}$. For sufficiently smooth solutions it defines the initial condition with $O(\tau^{2})$ accuracy. Results of computations for various weight parameters $\sigma$ are shown in Figs. 5 – 7. It is clearly seen that the accuracy of the approximation is increased for the optimal weight parameter $\sigma_{0}$.

Next we have investigated the influence of the initial condition for $w_{1}$. The accuracy of the approximation is further increased when the initial condition is computed using the algorithm (32). Results of computations for various weight parameters $\sigma$ are shown in Figs. 8 – 10.

Here we note that the observed convergence rates of the two-level and three-level schemes depend on the discrete regularity of the solution of the discrete fractional power problem and they are not reaching the maximal possible convergence rates of these schemes. As expected from the theoretical analysis (see, e.g. [26]), the convergence rate is increased for larger values of $\alpha$. The dependence on the regularity of the solution can be reduced by using geometrically refined time grids.

6 Conclusions

1. We have formulated the problem of finding the high order difference schemes for solving the nonstationary Cauchy type problem which is equivalent to the fractional power elliptic problem. The high order approximations are used to approximate the time dependence of the solution, while the elliptic operator is approximated by the standard finite element scheme.

2. The sufficient stability conditions are given for the two-level discrete schemes with weight parameters. The second order accuracy is proved for the symmetrical Crank-Nicolson type scheme.

3. The family of three-level symmetrical discrete schemes is constructed and investigated. It is proved that the second order approximation is valid for sufficiently smooth solutions. The initial condition on the first time level is computed by using the symmetrical two-level scheme.

4. It is shown that for a special weight parameter $\sigma_{0}$ we get the fourth-order three-level scheme. The value of this optimal parameter depends on the fractional power $\alpha$ of the elliptic operator. The initial condition on the first time level of the main grid is computed by using the symmetrical two-level scheme with a specially selected fine time grid.

5. The theoretical results are illustrated by results of numerical experiments. A two-dimensional problem is solved for the elliptic operator with the discontinuous sink term coefficient.

Acknowledgements

This work of second author was supported by the mega-grant of the Russian Federation Government (# 14.Y26.31.0013).