Numerical solution of time-dependent problems with fractional power elliptic operator

TL;DR

This paper develops numerical methods for solving unsteady space-fractional elliptic equations using finite element spatial discretization and Pade-type approximations for fractional powers, with demonstrated effectiveness through 2D tests.

Contribution

It introduces a novel approach combining finite element methods with Pade-type approximations for time-dependent fractional elliptic problems.

Findings

Effective finite element discretization in space.

Successful implementation of Pade-type approximations for fractional powers.

Numerical experiments confirm method accuracy for 2D problems.

Abstract

An unsteady problem is considered for a space-fractional equation in a bounded domain. A first-order evolutionary equation involves a fractional power of an elliptic operator of second order. Finite element approximation in space is employed. To construct approximation in time, standard two-level schemes are used. The approximate solution at a new time-level is obtained as a solution of a discrete problem with the fractional power of the elliptic operator. A Pade-type approximation is constructed on the basis of special quadrature formulas for an integral representation of the fractional power elliptic operator using explicit schemes. A similar approach is applied in the numerical implementation of implicit schemes. The results of numerical experiments are presented for a test two-dimensional problem.

Figures

Tables

Table 1. Table 1: The roots of equation ( 30 )

$k$	$g=1$	$g=10$	$g=100$
1	1.25578371	2.17949660	2.38090166
3	7.15579917	7.95688342	8.56783165

Table 2. Table 2: The spectrum bounds of the operator A 𝐴 A

$g$	$\delta ={\lambda }_1$	grid	${\underset{¯}{\delta }}_h$	${\overline{\delta }}_h$
		1	1.57959231369	4225.51507674
1	1.57699272630	2	1.57763558651	17104.1780271
		3	1.57715815735	74989.7519112
		1	4.76409956820	4252.23867499
10	4.75020542941	2	4.75363764524	17143.1279728
		3	4.75108440807	74989.7519112
		1	5.68846224707	7310.80621520
100	5.66869271459	2	5.67358161306	22017.7463507
		3	5.66994292109	74989.7519112

Table 3. Table 3: The upper bounds of the operator A ​ R M ​ ( A ) 𝐴 subscript 𝑅 𝑀 𝐴 AR_{M}(A)

$M$	$\overline{\gamma }(M,0.25)$	$\overline{\gamma }(M,0.5)$	$\overline{\gamma }(M,0.75)$
5	4.4602175	21.794966	142.00220
10	6.3106349	43.589932	401.45610
20	8.9256294	87.179864	1135.3565
40	12.623116	174.35973	3211.1792

Table 4. Table 4: Error of the solution for the explicit scheme on grid 2 ( μ = 10 , α = 0.5 formulae-sequence 𝜇 10 𝛼 0.5 \mu=10,\alpha=0.5 )

$M$	$N$	25	50	100	200
5	${\epsilon }_2$	0.00436648	0.00207328	0.00094905	0.00041937
	${\epsilon }_{\infty }$	0.01896717	0.00927482	0.00447550	0.00216564
10	${\epsilon }_2$	0.00507635	0.00277981	0.00164515	0.00108352
	${\epsilon }_{\infty }$	0.02186657	0.01217982	0.00738292	0.00499609
20	${\epsilon }_2$	0.00507902	0.00278251	0.00164787	0.00108627
	${\epsilon }_{\infty }$	0.02187724	0.01219044	0.00739354	0.00500671
40	${\epsilon }_2$	0.00507902	0.00278251	0.00164787	0.00108627
	${\epsilon }_{\infty }$	0.02187723	0.01219043	0.00739352	0.00500669

Table 5. Table 5: Error of the solution for various spatial grids ( μ = 10 , α = 0.5 , M = 20 formulae-sequence 𝜇 10 formulae-sequence 𝛼 0.5 𝑀 20 \mu=10,\alpha=0.5,M=20 )

grid	$N$	25	50	100	200
1	${\epsilon }_2$	0.00641387	0.00419465	0.00310833	0.00257568
	${\epsilon }_{\infty }$	0.02505950	0.01634587	0.01202969	0.00988175
2	${\epsilon }_2$	0.00507902	0.00278251	0.00164787	0.00108627
	${\epsilon }_{\infty }$	0.02187724	0.01219044	0.00739354	0.00500671
3	${\epsilon }_2$	0.00472777	0.00241442	0.00126921	0.00069981
	${\epsilon }_{\infty }$	0.02077071	0.01081355	0.00588308	0.00342987

Table 6. Table 6: Error of the solution for the implicit scheme on grid 2 ( μ = 10 , α = 0.5 formulae-sequence 𝜇 10 𝛼 0.5 \mu=10,\alpha=0.5 )

$M$	$N$	25	50	100	200
5	${\epsilon }_2$	0.00326281	0.00111138	0.00044581	0.00078676
	${\epsilon }_{\infty }$	0.01672817	0.00699940	0.00208775	0.00313725
10	${\epsilon }_2$	0.00394436	0.00173870	0.00063823	0.00018398
	${\epsilon }_{\infty }$	0.01975034	0.01002383	0.00511271	0.00264674
10	${\epsilon }_2$	0.00394696	0.00174126	0.00064056	0.00018400
	${\epsilon }_{\infty }$	0.01976054	0.01003442	0.00512346	0.00265756
10	${\epsilon }_2$	0.00394696	0.00174126	0.00064056	0.00018399
	${\epsilon }_{\infty }$	0.01976037	0.01003433	0.00512339	0.00265750

Equations

$$\frac{du}{dt}+{\mathcal{A}}^{\alpha }u=f(t),0<\alpha <1,0<t\le T,$$

$$u(0)=u_0,$$

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

$$k(\mathbf{x})\frac{\partial u}{\partial n}+g(\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}+g(\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})=\{u|u(\mathbf{x})\in L_2(\Omega ),\sum _{k=0}^{\infty }|(u,{\phi }_k){|}^2{\lambda }_k<\infty \}.$$

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

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

$$\frac{du}{dt}+{\mathcal{A}}^{\alpha }u=f(t),0<t\le T,$$

$$u(0)=u_0.$$

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

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

$$(Ay,v)=a(y,v),\forall y,v\in V^h.$$

$$A{\tilde{\phi }}_k={\tilde{\lambda }}_k$$

$${\tilde{\lambda }}_1\le {\tilde{\lambda }}_2\le ...\le {\tilde{\lambda }}_{M_h},\Vert {\tilde{\phi }}_k\Vert =1,k=1,2,...,M_h.$$

$$D(A)=\{y|y\in V^h,\sum _{k=0}^{M_h}|(y,{\tilde{\phi }}_k){|}^2{\tilde{\lambda }}_k<\infty \}.$$

$$A=A^{\ast }\ge {\underline{\delta }}_{h}I,{\underline{\delta }}_h>0,$$

$$A^{\alpha }y=\sum _{k=1}^{M_h}(y,{\tilde{\phi }}_k){\tilde{\lambda }}_k^{\alpha }{\tilde{\phi }}_k.$$

$$\frac{dw}{dt}+A^{\alpha }w=\psi (t),0<t\le T,$$

$$w(0)=w_0,$$

$$\left(\frac{dw}{dt},w\right)+(A^{\alpha }w,w)=(\psi ,w).$$

$$\left(\frac{dw}{dt},w\right)\le (\psi ,w).$$

$$(\psi ,w)\le \Vert \psi \Vert \Vert w\Vert .$$

$$\frac{d}{dt}\Vert w\Vert \le \Vert \psi \Vert .$$

$$\Vert w(t)\Vert \le \Vert w_0\Vert +{\int }_0^t\Vert \psi (\theta )\Vert d\theta .$$

$$\frac{w^{n+1}-w^n}{\tau }+A^{\alpha }{w}^n={\psi }^n,n=0,1,...,N-1,$$

Full text

††journal: arXiv

1 Introduction

Many applied mathematical models involve both sub-diffusion (fractional in time) and super-diffusion (fractional in space) operators (see, e.g., Podlubny (1998); Uchaikin (2013)). Super-diffusion problems are treated as evolutionary problems with a fractional power of an elliptic operator. For example, suppose that in a bounded domain $\Omega$ on the set of functions $u(\bm{x})=0,\ \bm{x}\in\partial\Omega$, there is defined the operator $\mathcal{A}$: $\mathcal{A}u=-\triangle u,\ \bm{x}\in\Omega$. We seek the solution of the Cauchy problem for the equation with the fractional power elliptic operator:

[TABLE]

[TABLE]

for a given $f(\bm{x},t)$, $u_{0}(\bm{x}),\ \bm{x}\in\Omega$ using the notation $f(t)=f(\cdot,t)$.

For approximation in space, we can apply finite volume or finite element methods oriented to using arbitrary domains and irregular computational grids (Knabner and Angermann (2003); Quarteroni and Valli (1994)). After this, we formulate the corresponding Cauchy problem with a fractional power of a self-adjoint positive definite discrete elliptic operator (see Bonito and Pasciak (2015); Szekeres and Izsák (2016)) — a fractional power of a symmetric positive definite matrix (Higham (2008)).

In the study of difference schemes for time-dependent problems of BVP for PDE, the general theory of stability (well-posedness) for operator-difference schemes (Samarskii (2001); Samarskii et al. (2002)) is in common use. At the present time, the exact (matching necessary and sufficient) conditions for stability are obtained for a wide class of two- and three-level difference schemes considered in finite-dimensional Hilbert spaces. We emphasize a constructive nature of the general theory of stability for operator-difference schemes, where stability criteria are formulated in the form of operator inequalities, which are easy to verify. In particular, the backward Euler scheme and Crank-Nicolson scheme are unconditionally stable for a non-negative operator.

Problems in numerical solving unsteady problems with fractional powers of operators appear in using the simplest explicit approximations in time. A practical implementation of such approach requires the matrix function-vector multiplication. For such problems, different approaches (see Higham (2008)) are available. Algorithms for solving systems of linear equations associated with fractional elliptic equations that are based on Krylov subspace methods with the Lanczos approximation are discussed, e.g., in Ilić et al. (2008). 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 Burrage et al. (2012). The simplest variant is associated with the explicit construction of the solution using the eigenvalues and eigenfunctions of the elliptic operator with diagonalization of the corresponding matrix (Bueno-Orovio et al. (2014); Ilic et al. (2006)). Unfortunately, all these approaches demonstrate very high computational complexity for multidimensional problems.

There does exist a general approach to solve approximately equations involving a fractional power of operators based on an approximation of the original operator and then taking fractional power of its discrete variant. Using the Dunford-Cauchy formula the elliptic operator is represented as a contour integral in the complex plane. Further, applying appropriate quadratures with integration nodes in the complex plane, it is necessary to select a proper method that involves only inversion of the original operator. The approximate operator is treated as a sum of resolvents (Gavrilyuk et al. (2004, 2005)) ensuring the exponential convergence of quadrature approximations. In Bonito and Pasciak (2015), there was presented a more promising variant of using quadrature formulas with nodes on the real axis, which are constructed on the basis of the corresponding integral representation for the power operator (see Krasnoselskii et al. (1976); Carracedo et al. (2001)). In this case, the inverse operator of the problem has an additive representation, where each term is an inverse of the original elliptic operator. A similar rational approximation to the fractional Laplacian operator is studied in Aceto and Novati (2017).

We have proposed (Vabishchevich (2015)) a numerical algorithm to solve an equation for fractional power elliptic operators that is based on a transition to a pseudo-parabolic equation. For an auxiliary Cauchy problem, the standard two-level 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 time steps is required to find a solution. This computational algorithm for solving equations with fractional power operators is promising for transient problems.

As for one-dimensional problems for the space-fractional diffusion equation, an analysis of stability and convergence for this equation was conducted in Jin et al. (2014) using finite element approximation in space. A similar study for the Crank-Nicolson scheme was conducted earlier in Tadjeran et al. (2006) using finite difference approximations in space. We highlight separately the works Huang et al. (2008); Sousa (2012); Meerschaert and Tadjeran (2004), where numerical methods for solving one-dimensional transient problems of convection and space-fractional diffusion equation are considered.

In Vabishchevich (2016c), an unsteady problem is considered for a space-fractional diffusion equation in a bounded domain. To construct approximation in time, regularized two-level schemes are used (see Vabishchevich (2014)). The numerical implementation is based on solving the equation with the fractional power of the elliptic operator using an auxiliary Cauchy problem for a pseudo-parabolic equation (Vabishchevich (2015)). Some more general unsteady problems are considered in Vabishchevich (2016a, b).

In the present work, standard two-level schemes are applied to solve numerically a Cauchy problem for an evolutionary equation of first order with a fractional power of the operator. The numerical implementation is based on the rational approximation of the operator at a new time-level. When implementing the explicit scheme, the fractional power of the operator is approximated on the basis of Gauss-Jacobi quadrature formulas for the corresponding integral representation. In this case, we have (see Frommer et al. (2014)) a Pade-type approximation of the power function with a fractional exponent. A similar approach is used when considering implicit schemes.

The paper is organized as follows. The formulation of an unsteady problem containing a fractional power of an elliptic operator is given in Section 2. Here finite element approximations in space are also discussed. In Section 3, we construct the explicit approximation in time and investigate its stability. The numerical implementation is based on the rational approximation of the fractional power operator. Implicit schemes are considered in Section 4. The results of numerical experiments are described in Section 5. At the end of the work the main results of our study are summarized.

2 Problem formulation

In a bounded polygonal domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ with the Lipschitz continuous boundary $\partial\Omega$, we search the solution for the problem with a fractional power of an elliptic operator. Define the elliptic operator as

[TABLE]

with coefficients $0<\underline{k}\leq k({\bm{x}})\leq\overline{k}$, $c({\bm{x}})\geq 0$. The operator $\mathcal{A}$ is defined on the set of functions $u({\bm{x}})$ that satisfy on the boundary $\partial\Omega$ the following conditions:

[TABLE]

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

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

[TABLE]

For the spectral problem

[TABLE]

[TABLE]

we have

[TABLE]

and the eigenfunctions $\varphi_{k},\ \|\varphi_{k}\|=1,\ k=1,2,...$ form a basis in $L_{2}(\Omega)$. Therefore,

[TABLE]

Let the operator $\mathcal{A}$ be defined in the following domain:

[TABLE]

The operator $\mathcal{A}$ is self-adjoint and positive definite:

[TABLE]

where $I$ is the identity operator in $H$. For $\delta$, we have $\delta=\lambda_{1}$. In applications, the value of $\lambda_{1}$ is unknown (the spectral problem must be solved). Therefore, we assume that $\delta\leq\lambda_{1}$ in (3). Let us assume for the fractional power of the operator $\mathcal{A}$:

[TABLE]

More general and mathematically complete definition of fractional powers of elliptic operators is given in Carracedo et al. (2001); Yagi (2009).

We seek the solution of a Cauchy problem for the evolutionary first-order equation with the fractional power of the operator $\mathcal{A}$. The solution $u(\bm{x},t)$ satisfies the equation

[TABLE]

and the initial condition

[TABLE]

The key issue in the study of computational algorithms 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.

To solve numerically the problem (4), (5), we employ finite element approximations in space (see, e.g., Brenner and Scott (2008); Thomée (2006)). For (1) and (2), we define the bilinear form

[TABLE]

By (3), we have

[TABLE]

Define the subspace of finite elements $V^{h}\subset H^{1}(\Omega)$ and the discrete elliptic operator $A$ as

[TABLE]

The fractional power of the operator $A$ is defined similarly to $\mathcal{A}^{\alpha}$. For the spectral problem

[TABLE]

we have

[TABLE]

The domain of definition for the operator $A$ is

[TABLE]

The operator $A$ acts on the finite dimensional space $V^{h}$ defined on the domain $D(A)$ and, similarly to (3),

[TABLE]

where $\underline{\delta}_{h}\leq\lambda_{1}\leq\widetilde{\lambda}_{1}$. For the fractional power of the operator $A$, we suppose

[TABLE]

The use of finite element approximations for fractional power elliptic operators is discussed in detail, for instance, in the works Acosta and Borthagaray (2017); Szekeres and Izsák (2016).

For the problem (4), (5), we put into the correspondence the operator equation for $w(t)\in V^{h}$:

[TABLE]

[TABLE]

where $\psi(t)=Pf(t)$, $w_{0}=Pu_{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 $\Omega$:

[TABLE]

In view of the self-adjointness and positive definiteness of the operator $A^{\alpha}$, we have

[TABLE]

The right-hand side can be evaluated by the inequality

[TABLE]

By virtue of this, we have

[TABLE]

The latter inequality leads us to the desired a priori estimate:

[TABLE]

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).

3 Explicit scheme

To solve numerically the problem (7), (8), we use simplest explicit and implicit two-level schemes. Let $\tau$ be a step of a uniform grid in time such that $w^{n}=w(t^{n}),\ t^{n}=n\tau$, $n=0,1,...,N,\ N\tau=T$. It seems reasonable to begin with the simplest explicit scheme

[TABLE]

[TABLE]

Advantages and disadvantages of explicit schemes for the standard parabolic problem ($\alpha=1$) are well-known, i.e., these are a simple computational implementation and a time step restriction (see, e.g., Samarskii (2001); Samarskii et al. (2002)). In our case ($\alpha\neq 1$), the main drawback (conditional stability) remains, whereas the advantage in terms of implementation simplicity does not exist. The approximate solution at a new time-level is determined via (10) as

[TABLE]

Thus, we must calculate $A^{\alpha}w^{n}$. In view of these problems, considering the scheme (10), it is more correct to speak about the scheme with the explicit approximations in time in contrast to the standard fully explicit scheme.

The numerical implementation of (12) is based on the following representation:

[TABLE]

We construct a numerical algorithm that employ the rational approximation of the operator

[TABLE]

In this case, we solve standard problems that are related to the operator $A$.

We use an approximation for $A^{-\beta}$ based on integral representation of a self-adjoint and positive definite operator $A$ (see, e.g., Krasnoselskii et al. (1976); Carracedo et al. (2001)):

[TABLE]

The approximation of $A^{-\beta}$ is based on the use of one or another quadrature formulas for the right-hand side of (13). Various possibilities in this direction are discussed in Bonito and Pasciak (2015). One possibility is special Gauss-Jacobi quadrature formulas considered in Frommer et al. (2014); Aceto and Novati (2017). Just this approximation of the operator $A^{-\beta}$ is used in the present work.

To achieve higher accuracy in approximating the the right-hand side (13), it is natural to focus on the use of Gauss quadrature formulas. Some possibilities of constructing quadratures directly for half-infinite intervals are investigated, for example, in the work Gautschi (1991). The classical Gauss quadrature formulas can be used via introducing a new variable of integration.

Suppose (see Frommer et al. (2014)) that

[TABLE]

From (13), we have

[TABLE]

To approximate the right-hand side of (14), we apply the Gauss-Jacobi quadrature formula with the weight $(1-\eta)^{\tilde{\alpha}}(1+\eta)^{\tilde{\beta}})$ (see Ralston and Rabinowitz (2001)):

[TABLE]

Here $\eta_{1},\eta_{2},...,\eta_{M}$ are the roots of the Jacobi polynomial $J_{M}(\eta;\tilde{\alpha},\tilde{\beta})$ of degree $M$. The weights $\omega_{1},\omega_{2},...,\omega_{M}$ are given by the formula

[TABLE]

where $\Gamma$ denotes the gamma function.

For the fractional power of the operator $A$, from (14), we get

[TABLE]

where

[TABLE]

In view of (15), the approximate solution of the problem $r^{n}=A^{\alpha-1}w^{n}$ is associated with solving $M$ standard problems $r^{n}_{m}=(c_{m}I+A)^{-1}w^{n},\ m=1,2,...,M$.

Instead of (12), we employ the scheme

[TABLE]

Let us consider the stability conditions for the scheme (11), (16). For a finite-dimensional self-adjoint operator $A$, in addition to the lower bound (6), the following upper bound holds:

[TABLE]

where $\overline{\delta}_{h}=\mathcal{O}(h^{-2})$. Thus

[TABLE]

Similar estimates we have also for $AR_{M}(A)$:

[TABLE]

with some $\underline{\gamma}_{h},\overline{\gamma}_{h}>0$.

Theorem 1

If

[TABLE]

then the scheme (11), (16) is stable in $H$ and the solution satisfies the following estimate:

[TABLE]

Proof 1

From (16), we directly obtain

[TABLE]

By (6), (17), (18), we get

[TABLE]

under the restrictions (19). In view of this, we have the level-wise estimate

[TABLE]

that results in the estimate (20) for the stability of the solution on the right-hand side and the initial conditions.

It should be noted that the estimate (20) for the difference scheme (11), (16) is a discrete analog of the a priori estimate (9) for the problem (7), (8).

Special attention (see Frommer et al. (2014); Aceto and Novati (2017)) should be given to the problem of choosing the parameter $\mu$ in (14). Taking into account the definition of the operator $A$, we are interested in the best approximation of $A^{-\beta}$ for the smallest (principal) eigenvalue $\widetilde{\lambda}_{1}$. In Frommer et al. (2014), there is established a remarkable fact that $R_{M}(z)$ corresponds to a Pade-type approximation for the function $z^{-\beta}$ with expansion point $\mu$. Thus, the optimal choice corresponds to the selection $\mu=\underline{\delta}_{h}$. In this case, in (18), we have $\underline{\gamma}=\underline{\delta}_{h}^{\alpha}$. The computational complexity of finding $\underline{\delta}_{h}=\widetilde{\lambda}_{1}$ (the principal eigenvalue of a discrete self-adjoint elliptic operator of second order) is not high. To this end, it is possible to use standard algorithms (see, e.g., Saad (2011)) and the corresponding software (see Hernandez et al. (2005)).

The function $zR_{M}(z)$ for $z\geq z_{0}>0$ is a positive and monotonically increasing function. In view of this, taking into account (15), we have

[TABLE]

where

[TABLE]

From theorem 1, it follows that for the explicit scheme (11), (16) the time-step restrictions do not depend on discretization in space, but depend on the power $\alpha$ and the number of approximation nodes $M$.

4 Implicit scheme

Unconditionally stable schemes are constructed on the basis of implicit approximations in time. Here we consider standard two-level schemes with weights (Samarskii (2001); Samarskii et al. (2002)). For a constant weight parameter $\sigma\ (0<\sigma\leq 1)$, we define

[TABLE]

Instead of (10), let us consider the implicit scheme

[TABLE]

For $\sigma=1/2$, the difference scheme (21) is the symmetric scheme (the Crank-Nicolson scheme). It approximates the differential problem with the second order by $\tau$, whereas for other values of $\sigma$, we have only the first order.

Rewrite the scheme (21) in the form

[TABLE]

In view of this, the transition to a new time-level involves the solution of the problem

[TABLE]

For this problem, we construct the rational approximation of the operator

[TABLE]

The necessary approximation is based on the following integral representation:

[TABLE]

taken from the work Carracedo et al. (2001) (see Proposition 5.3.2).

Using the new variable $\theta$, from (22), we arrive at the representation

[TABLE]

where

[TABLE]

Further, the Gauss quadrature formula is used (see Gautschi (2004)) with the weight function

[TABLE]

From (23), we get

[TABLE]

Thereby $R_{M}(A;0)=R_{M}(A)$.

For $\sigma>0$, an approximate solution is obtained from

[TABLE]

Theorem 2

The difference scheme (11), (25) for $\sigma\geq 0.5$ and

[TABLE]

is unconditionally stable in $H$ and its solution satisfies the a priori estimate

[TABLE]

Proof 2

The use of (25) means that instead of (21) we employ the scheme

[TABLE]

where, in view of (26), we have

[TABLE]

Multiplying (29) scalarly in $H$ by $w^{n+\sigma}$, we get

[TABLE]

Taking into account

[TABLE]

in view of $0.5\leq\sigma\leq 1$, for the left-hand side of (29), we obtain

[TABLE]

For the right-hand side of (29), we have

[TABLE]

In view of this, from (29), we get the inequality

[TABLE]

which provides the estimate (27).

5 Numerical experiments

The test problem under the consideration is constructed using the exact solution of the problem in the unit circle (see Vabishchevich (2016b)). The computational domain is a quarter of the circle (see Fig. 1). Consider the equation

[TABLE]

with the boundary conditions

[TABLE]

[TABLE]

We study the case, where the solution depends only on $r$, and $r=(x_{1}^{2}+x_{2}^{2})^{1/2}$. By virtue of this

[TABLE]

[TABLE]

The solution of the spectral problem

[TABLE]

[TABLE]

is well-known (see, e.g., Polyanin (2002); Carslaw and Jaeger (1986)). Eigenfunctions are represented as zero-order Bessel functions:

[TABLE]

whereas eigenvalues $\lambda_{k}=\nu_{k}^{2}$, where $\nu_{k},\ k=1,2,...$ are roots of the equation

[TABLE]

The general solution of the homogeneous ($f(t)=0$) Cauchy problem for equation (4) is

[TABLE]

To study the accuracy of the approximate solution of the time-dependent equation with the fractional power of an elliptic operator, we use the exact solution

[TABLE]

The values of the roots $\nu_{1},\ \nu_{3}$ for different values of the boundary condition $\mu$ are given in Table 1. The exact solution for $T=0.25$ at different values of $g$ and $\alpha$ is shown in Figs. 2 and 3.

Predictions were performed on a sequence of refining grids, which are shown in Fig. 4. The numerical solution is compared with the exact one at the final time moment $u(\bm{x},T)$. Error estimation is performed in $L_{2}(\Omega)$ and $L_{\infty}(\Omega)$:

[TABLE]

The finite element approximation in space is based on the use of continuous $P_{1}$ Lagrange element, namely, piecewise-linear elements. The calculations were performed using the computing platform FEniCS for solving partial differential equations (website http://fenicsproject.org, Logg et al. (2012); Alnæs et al. (2015)). To solve spectral problems with symmetrical matrices, we use the SLEPc library (Scalable Library for Eigenvalue Problem Computations, http://slepc.upv.es, Hernandez et al. (2005)). We apply the Krylov-Schur algorithm, a variation of the Arnoldi method, proposed by Stewart (2001).

Table 2 presents the lower and upper bounds of the operator spectrum (see (6), (17)) on various grids for different values of the parameter $g$ in the boundary condition. A comparison with the exact minimum eigenvalue demonstrates good accuracy in evaluation of $\underline{\delta}_{h}$, which increases with the grid refinement. It is easy to see a significant dependence of the maximum eigenvalue on the grid size.

Peculiarities of the approximation (15) are illustrated by the accuracy of the approximation of the function $z^{-\beta}$ for $z\geq z_{0}=\delta$ (see (3)). Figure 5 shows the absolute error arising from the approximation of $z^{-\beta}$ by the function $R_{M}(z)$ for $\beta=0.5$ and $g=10$. In this case, $\mu=\delta$ (see Table 1) and $R_{M}(z_{0})=z_{0}^{-\beta}$. We see higher accuracy near the left boundary $z=z_{0}$, whereas for large $z$, the approximation accuracy decreases. The effect of increasing accuracy with increasing number of nodes of the quadrature formula is clearly observed.

Decreasing the approximation accuracy at $z\approx z_{0}$, we can increase the accuracy for other values of $z$. For example, Fig. 6 demonstrates the approximation accuracy for $\mu=50$. In this case $R_{M}(\mu)=\mu^{-\beta}$. We need to have good accuracy for small $z$, and therefore, in calculations, we are guided by the choice of $\mu=\delta$. The dependence of the approximation accuracy of the function $z^{-\beta}$ on the value of $\beta$ is shown in Figs. 7 and 8. In these figures, we observe a significant drop in accuracy with decreasing $\beta$.

The numerical implementation of the explicit scheme (11), (16) involves the approximation of the operator $A^{\alpha}$ by the expression $AR_{M}(A)(\beta=1-\alpha)$. Peculiarities of this approximation at $\alpha=0.5$, $g=10$ are shown in Fig. 9. It should be noted that the operator $AR_{M}(A)$ is bounded and the constant $\overline{\gamma}(M,\alpha)$ on the right-hand side of (18) for the corresponding values of $M$ is presented via dotted lines.

The upper bounds of the operator $AR_{M}(A)$ are given in Table 3 for $g=10$. Increasing $\overline{\gamma}(M,\alpha)$ with increasing the number of quadrature formula nodes $M$ results from increasing the accuracy of approximation of the unbounded operator $A^{\alpha}$. As $\alpha$ decreases, the value of $\overline{\gamma}(M,\alpha)$ decreases drastically.

Now we present numerical results obtained using the explicit scheme (11), (16). We confine ourselves to the case $\alpha=0.5$ with the value of the boundary condition parameter $g=10$. It is interesting to identify the dependence of accuracy on grids in space and time. In our case, we should also study the influence of the number of quadrature formula nodes $M$. Table 4 demonstrates the numerical solution convergence for decreasing the time step and increasing the accuracy of approximation of the fractional power operator. Increasing the accuracy of an approximate solution due to spatial grid refinement is shown in Table 5.

The numerical implementation of implicit schemes is associated (see (21)) with the function $(\nu+z^{\alpha})^{-1}$. It should be noted that $\nu=(\sigma\tau)^{-1}$ and therefore, this parameter is large enough in numerical solving unsteady problems. The function $(\nu+z^{\alpha})^{-1}$, which corresponds to our test problem for $\alpha=0.5$, is shown in Fig. 10. As we noted earlier, if $\nu=0$, then the optimal value is $\mu=\delta$. This value is also used in our calculations for $\nu>0$.

Figure 11 shows the approximating function $R_{M}(z;\nu)$ for $\nu=200$ with $\mu=\delta$. The approximation accuracy for various values of $\nu$ is presented in Figs. 12 and 13. Operator approximations were designed using the package ORTHPOL (see Gautschi (1994)) developed for constructing Gauss quadrature formulas with an arbitrary weight function.

The accuracy of the approximate solution of the test problem obtained using the implicit scheme (11), (25) was investigated for $\alpha=0.5$ and $g=10$. For the fully implicit scheme ($\sigma=1$), Table 6 demonstrates the dependence of the solution accuracy on the grid in time for various numbers of the quadrature formula nodes $M$.

Conclusion

There is considered a nonclassical problem with the initial data, which is described by an evolutionary equation of first order with a fractional power of an elliptic operator. The multidimensional problem is approximated in space using standard finite element piecewise-linear approximations. An a priori estimate for stability with respect to the initial data and the right-hand side is provided. 2. 2.

The explicit scheme is implemented using a Pade-type approximation for the fractional power elliptic operator. Sufficient conditions for the stability of the explicit scheme are formulated. They do not depend on spatial grid steps. 3. 3.

Rational approximation is employed to implement implicit schemes. It is based on a computational generation of Gauss quadrature formulas for an integral representation of the operator of transition to a new time-level. 4. 4.

Possibilities of the proposed algorithms were demonstrated through numerical solving a test two-dimensional problem.

Acknowledgements

The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).