A Singularly Perturbed Boundary Value Problems with Fractional ...
A Singularly Perturbed Boundary Value Problems with Fractional Powers of Elliptic Operators
arXiv:1604.04427v1 [cs.NA] 15 Apr 2016
Petr N. Vabishchevich1,2 1
2
Nuclear Safety Institute, 52, B. Tulskaya, 115191 Moscow, Russia North-Eastern Federal University, 58, Belinskogo, 677000 Yakutsk, Russia
Abstract. A boundary value problem for a fractional power 0 < ε < 1 of the second-order elliptic operator is considered. The boundary value problem is singularly perturbed when ε → 0. It is solved numerically using a time-dependent problem for a pseudo-parabolic equation. For the auxiliary Cauchy problem, the standard two-level schemes with weights are applied. The numerical results are presented for a model two-dimensional boundary value problem with a fractional power of an elliptic operator. Our work focuses on the solution of the boundary value problem with 0 < ε 1.
1
Introduction
Non-local applied mathematical models based on the use of fractional derivatives in time and space are actively discussed in the literature [2, 11]. 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 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 A: Au = −4u, x ∈ Ω. We seek the solution of the problem for the equation with the fractional power of an elliptic operator: Aε u = f, with 0 < ε < 1 for a given f (x), x ∈ Ω. To solve problems with the fractional power of an elliptic operator, we can apply finite volume or finite element methods oriented to using arbitrary domains discretized by irregular computational grids [12, 15]. The computational realization is associated with the implementation of the matrix function-vector multiplication. For such problems, different approaches [7] 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 [10]. 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 [6]. 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 [5, 8, 9]. Unfortunately, all these approaches demonstrates too high computational complexity for multidimensional problems. We have proposed [19] 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-state solution. This computational algorithm for solving equations with fractional powers of operators is promising when considering transient problems. The boundary value problem for the fractional power of an elliptic operator is singularly perturbed when ε → 0. To solve it numerically, we focus on numerical methods that are designed for classical elliptic problems of convection-diffusionreaction [14, 17]. In particular, the main features are taken into account via using locally refining grids. The standard strategy of goal-oriented error control for conforming finite element discretizations [1, 3] is applied.
2
Problem formulation
In a bounded polygonal domain Ω ⊂ Rm , 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 Au = −div(k(x)grad u)
(1)
with coefficient 0 < k1 ≤ k(x) ≤ k2 . The operator A is defined on the set of functions u(x) that satisfy on the boundary ∂Ω the following conditions: u(x) = 0,
x ∈ ∂Ω.
(2)
In the Hilbert space H = L2 (Ω), we define the scalar product and norm in the standard way: Z (u, v) = u(x)v(x)dx, kuk = (u, u)1/2 . Ω
For the spectral problem Aϕk = λk ϕk , ϕk (x) = 0,
x ∈ Ω, x ∈ ∂Ω,
we have λ1 ≤ λ2 ≤ ...,
and the eigenfunctions ϕk , kϕk k = 1, k = 1, 2, ... form a basis in L2 (Ω). Therefore, ∞ X u= (u, ϕk )ϕk . k=1
Let the operator A be defined in the following domain: D(A) = {u | u(x) ∈ L2 (Ω),
∞ X
|(u, ϕk )|2 λk < ∞}.
k=0
Under these conditions the operator A is self-adjoint and positive defined: A = A∗ ≥ δI,
δ > 0,
(3)
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 fractional power of the operator A ∞ X (u, ϕk )λεk ϕk . Aε u = k=0
We seek the solution of the problem with the fractional power of the operator A. The solution u(x) satisfies the equation Aε u = f,
(4)
with 0 < ε < 1 for a given f (x), x ∈ Ω. The key issue in the study of the computational algorithm for solving the problem (4) is to establish the stability of the approximate solution with respect to small perturbations of the right-hand side in various norms. In view of (3), the solution of the problem (4) satisfies the a priori estimate kuk ≤ δ −ε kf k,
(5)
which is valid for all 0 < ε < 1. The boundary value problem for the fractional power of the elliptic operator (4) demonstrates a reduced smoothness when ε → 0. For the solution, we have (see, e.g., [20]) the estimate kuk2ε ≤ Ckf k, with 0 ≤ ε < 1/2, is k · k2ε is the norm in H 2ε (Ω). For the limiting solution, we have u0 (x) = f (x), x ∈ Ω. Thus, a singular behavior of the solution of the problem (4) appears with ε → 0 and is governed by the right-hand side f (x).
3
Discretization in space
To solve numerically the problem (4), we employ finite-element approximations in space [4, 18]. For (1) and (2), we define the bilinear form Z a(u, v) = k grad u grad v. Ω
By (3), we have a(u, u) ≥ δkuk2 . Define a subspace of finite elements V h ⊂ H01 (Ω). Let xi , i = 1, 2, ..., Mh be triangulation points for the domain Ω. Define pyramid function χi (x) ⊂ V h , i = 1, 2, ..., Mh , where 1, if i = j, χi (xj ) = 0, if i 6= j. For v ∈ Vh , we have Mh X
v(x) =
vi χi (x),
i=i
where vi = v(xi ), i = 1, 2, ..., Mh . We have defined Lagrangian finite elements of first degree, i.e., based on the piecewise-linear approximation. We will also use Lagrangian finite elements of second degree defined in a similar way. We define the discrete elliptic operator A as (Ay, v) = a(y, v),
∀ y, v ∈ V h .
The fractional power of the operator A is defined similarly to Aε . For the spectral problem ek Aϕ ek = λ we have e1 ≤ λ e2 ≤ ... ≤ λ eM , λ h
kϕ ek k = 1,
k = 1, 2, ..., Mh .
The domain of definition for the operator A is h
D(A) = {y | y ∈ V ,
Mh X
ek < ∞}. |(y, ϕ ek )|2 λ
k=0
The operator A acts on a finite dimensional space V h defined on the domain D(A) and, similarly to (3), A = A∗ ≥ δI,
δ > 0,
e1 . For the fractional power of the operator A, we suppose where δ ≤ λ1 ≤ λ Aε y =
Mh X k=1
eε ϕ (y, ϕ ek )λ k ek .
(6)
For the problem (4), we put into the correspondence the operator equation for w(t) ∈ V h : Aε w = ψ, (7) where ψ = P f with P denoting L2 -projection onto V h . For the solution of the problem (6), (7), we obtain (see (5)) the estimate kwk ≤ δ −ε kψk,
(8)
for all 0 < ε < 1.
4
Singularly perturbed problem for a diffusion-reaction equation
The object of our study is associated with the development of a computational algorithm for approximate solving the singularly perturbed problem (4). After constructing a finite element approximation, we arrive at equation (7). Features of the solution related to a boundary layer are investigated on a model singularly perturbed problem for an equation of diffusion-reaction. The key moment is associated with selecting adaptive computational grids (triangulations). In view of ε Aε = (exp(ln A)) = I + ε ln A + O(ε2 ), we put the problem (7) into the correspondence with solving the equation εAu + u = ψ.
(9)
The equation (9) corresponds to solving the Dirichlet problem (see the condition (2)) for the diffusion-reaction equation − ε div(k(x)grad u) + u = f (x),
x ∈ Ω.
(10)
Basic computational algorithms for the singularly perturbed boundary problem (2), (10) are considered, for example, in [14, 17]. In terms of practical applications, the most interesting approach is based on an adaptation of a computational grid to peculiarities of the problem solution via a posteriori error estimates. Among main approaches, we highlight the strategy of the goal-oriented error control for conforming finite element discretizations [1,3], which is applied to approximate solving boundary value problems for elliptic equations. The strategy of goal-oriented error control is based on choosing a calculated functional. The accuracy of its evaluation is tracked during computations. In our Dirichlet problem for the second-order elliptic equation, the solution is varied drastically near the boundary. So, it seems natural to control the accuracy of calculations for the normal derivatives of the solution (fluxes) across the boundary or a portion of it. Because of this, we put Z G(u) = − εk(x)(grad u · n)d x, ∂Ω
where n is the outward normal to the boundary. An adaptation of a finite element mesh is based on an iterative local refinement of the grid in order to evaluate the goal functional with a given accuracy η on the deriving approximate solution uh , i.e., |G(u) − G(uh )| ≤ η. To conduct our calculations, we used the FEniCS framework (see, e.g., [13]) developed for general engineering and scientific calculations via finite elements. Features of the goal-oriented procedure for local refinement of the computational grid are described in [16] in detain. Here, we consider only a key idea of the adaptation strategy of finite element meshes, which is associated with selecting the goal functional. The model problem (2), (10) is considered with k(x) = 1,
f (x) = (1 − x1 )x22 ,
in the unit square (Ω = (0, 1) × (0, 1)). The threshold of accuracy for calculating the functional G(u) is defined by the value of η = 10−5 . As an initial mesh, there is used the uniform grid obtained via division by 8 intervals in each direction (step 0 — 128 cells). First, Lagrangian finite elements of first order have been used in our calculations. For this case, the improvement of the goal functional during the iterative procedure of adaptation is illustrated by the data presented in Table 1. Table 2 demonstrates values of the goal functional G(uh ) calculated on the final computational grid, the number of vertices of this final grid and the number of adaptation steps for solving the problem at various values of the small parameter ε. These numerical results demonstrate the efficiency of the proposed strategy for goal-oriented error control for conforming finite element discretizations applied to approximate solving singular perturbed problems of diffusion-reaction (2), (10). Table 1. Calculation of the goal functional during adaptation steps ε Step of adaptation s 0 1 2 3 4 5 6 7 8 9 10
10−1 G(uh ) 0.087608 0.110432 0.116155 0.119766 0.122702 0.125653 0.127950 0.128835 0.129542 0.129940 0.130149
Mh 81 97 140 222 384 694 1235 2179 3841 6540 11040
10−3 G(uh ) Mh 0.0056973 81 0.0107507 98 0.0129506 132 0.0155597 195 0.0175113 305 0.0194985 466 0.0210232 754 0.0221562 1279 0.0229284 2132 0.0234492 3753 0.0237487 6626
10−5 G(uh ) Mh 0.00006643 81 0.00015584 95 0.00023996 120 0.00035644 164 0.00050472 225 0.00068154 349 0.00090839 550 0.00115091 853 0.00137740 1242 0.00161273 1865 0.00181249 2711
Table 2. Adaptation for various values of ε ε 10−1 10−2 10−3 10−4 10−5
Goal functional G(uh ) 0.130396 0.064867 0.024191 0.008061 0.002580
Number of vertices 51868 72297 90170 67476 99003
Number of adaptation steps s 13 14 15 16 18
Next, similar results have been obtained using Lagrangian finite elements of second order. For this case, summary data are presented in Table 3. As expected, the desired accuracy η = 10−5 is reached on adaptive meshes of smaller sizes than in the case of Lagrangian finite elements of first order (see Table 2 for a comparison). Table 3. Adaptation for Lagrangian elements of second order Goal functional G(uh ) 0.130423 0.064884 0.024184 0.008076 0.002574
ε 10−1 10−2 10−3 10−4 10−5
5
Number of vertices 3574 5137 6573 12775 18501
Number of adaptation steps s 7 8 9 11 12
Numerical algorithm for the problem with a fractional power
An approximate solution of the problem (7) is sought as a solution of an auxiliary pseudo-time evolutionary problem [19]. Assume that y(t) = δ ε (t(A − δI) + δI)−ε y(0). Therefore y(1) = δ ε A−ε y(0) and then w = y(1). The function y(t) satisfies the evolutionary equation (tD + δI)
dy + εDy = 0, dt
where D = A − δI.
0 < t ≤ 1,
(11)
By (6), we get D = D∗ > 0.
(12)
We supplement (11) with the initial condition y(0) = δ −ε ψ.
(13)
The solution of equation (7) can be defined as the solution of the Cauchy problem (11)–(13) at the final pseudo-time moment t = 1. For the solution of the problem (11), (13), it is possible to obtain various a priori estimates. The elementary estimate that is consistent with the estimate (8) have the form ky(t)k ≤ ky(0)k. (14) To get (14), multiply scalarly equation (11) by εy + tdy/dt. To solve numerically the problem (11), (13), we use the simplest implicit two-level scheme. Let τ be a step of a uniform grid in time such that y n = y(tn ), tn = nτ , n = 0, 1, ..., N, N τ = 1. Let us approximate equation (11) by the implicit two-level scheme (tσ(n) D + δI)
y n+1 − y n + εDy σ(n) = 0, τ
n = 0, 1, ..., N − 1,
y 0 = δ −ε ψ.
(15) (16)
We use the notation tσ(n) = σtn+1 + (1 − σ)tn ,
y σ(n) = σy n+1 + (1 − σ)y n .
For σ = 0.5, the difference scheme (15), (16) approximates the problem (11), (12) with the second order by τ , whereas for other values of σ, we have only the first order. Theorem 1. For σ ≥ 0.5 the difference scheme (15), (16) is unconditionally stable with respect to the initial data. The approximate solution satisfies the estimate ky n+1 k ≤ ky 0 k, n = 0, 1, ..., N − 1. (17) Proof. Rewrite equation (15) in the following form: y n+1 − y n y n+1 − y n δ + D εy σ(n) + tσ(n) = 0. τ τ Multiplying scalarly it by εy σ(n) + tσ(n) in view of (12), we arrive at
y n+1 − y n , τ
y n+1 − y n σ(n) ,y τ
≤ 0.
We have y
n+1 1 y − yn 1 = σ− τ + (y n+1 + y n ). 2 τ 2
σ(n)
If σ ≥ 0.5, then ky n+1 k ≤ ky n k,
n = 0, 1, ..., N − 1.
Thus, we obtain (17). The key point in approximate solving singularly perturbed boundary value problems is associated with mesh adaptation. In the case of solving the problem (4), we use finite element approximations and proceed to the problem (7) and then formulate the Cauchy problem (11), (13) approximated by the scheme (15), (16). In our case, singularity is associated only with spatial variables. The decomposition of the solution of the problem (11), (13) by eigenfunctions of the operator A results in y(t) =
Nh X
ak (t)ϕ ek .
k=1
For coefficients ak (t), we get ek − δ)t)−ε , ak (t) = (ψ, ϕ ek )(δ + (λ
k = 1, 2, ..., Mh .
Because of this, errors in specifying the initial conditions monotonically decrease for increasing t. A similar behavior demonstrates an approximate solution of the Cauchy problem (11), (13) obtained using the fully implicit scheme with σ = 1 in (15), (16). For the Crank-Nicolson scheme (i.e., σ = 0.5 in (15), (16)), we cannot guarantee a monotone decrease of errors in time, but the error at t = 1 will not be more than at t = 0. The practical significance of such an analysis is that it provides us a simple adaptation strategy for computational grids in solving the problem (11), (13), namely, spatial mesh adaptation is conducted at the first time step of calculations.
6
Solution of a model problem
Below, there are presented some results of numerical solving the problem (7) for small values of ε. A computational algorithm must track a singular behavior of the solution, which is directly related to the singular behavior of the righthand side f (x). Let us consider the problem (2), (10) in the unit square Ω = (0, 1) × (0, 1) with 1 − x2 x1 2 x2 − exp − . k(x) = 1, f (x) = 1 − x1 − exp − µ µ The singularity of the right-hand side (the singularity of a numerical solution of the problem with a fractional power of an elliptic operator) results from existing a boundary layer at low values of µ.
An adaptation of the computational grid is performed during the calculation of the first time step using the two-level scheme (15), (16). For the basic variant, it is assumed that ε = 10−2 , µ = 10−2 , the initial uniform spatial grid contains 8 intervals in each direction and the time step is τ = 10−2 . The parameter δ = 2π 2 corresponds the minimal eigenvalue of the elliptic operator A. Mesh adaptation is carried out taking into account peculiarities of the right-hand side and the goal functional defined in the form Z G(u; t = τ ) = −
k(x)(grad u · n)d x. ∂Ω
Next, the problem (2), (10) is solved using the derived grid in space and the uniform grid in time. Thus, we apply the simplest one-stage starting adaptation of the computational grid for numerical solving the unsteady problem. Lagrangian finite elements of second order are used. For time-stepping, the Crank-Nicolson (σ = 0.5 in (15)) scheme is utilized. The sequence of calculated adaptive grids is shown in Fig. 1. Note that this sequence is weakly dependent on the choice of a time step. The goal functional dynamics for different levels of adaptation is presented in Table 4. The problem is solved with different values of ε.
Table 4. Calculation of the goal functional during adaptation steps ε Step of adaptation s 0 1 2 3 4 5 6 7 8 9 10
10−1 G(uh ; t = τ ) Mh 16.0955 289 24.3875 315 31.1692 399 37.0996 559 42.0854 833 45.7594 1270 48.6087 1849 50.3070 2753 51.2621 4067 51.9766 5965 52.2862 9201
10−2 G(uh ; t = τ ) Mh 21.8932 289 33.7810 315 43.8893 399 52.8249 559 60.4373 837 66.2019 1282 70.4881 1885 73.1491 2778 74.6778 4120 75.8362 5968 76.3402 9235
10−3 G(uh ; t = τ ) Mh 22.5773 289 34.9016 315 45.4226 399 54.7328 559 62.2786 834 68.4363 1264 73.2101 1889 75.9998 2774 77.5762 4125 78.7648 6028 79.2942 9261
Acknowledgements This work was supported by the Russian Foundation for Basic Research (projects 14-01-00785, 15-01-00026).
0 — 128 cells
1 — 140 cells
2 — 180 cells
3 — 256 cells
4 — 388 cells
5 — 599 cells
6 — 886 cells
7 — 1313 cells
Fig. 1. The grid obtained at succesive steps of adaptation
References 1. Ainsworth, M., Oden, J.T.: A Posteriori Error Estimation in Finite Element Analysis. Wiley, New York (2000) 2. Baleanu, D.: Fractional Calculus: Models and Numerical Methods. World Scientific, New York (2012) 3. Bangerth, W., Rannacher, R.: Adaptive Finite Element Methods for Differential Equations. Birkh¨ auser, Basel (2003) 4. Brenner, S.C., Scott, L.R.: The mathematical theory of finite element methods. Springer, New York (2008) 5. Bueno-Orovio, A., Kay, D., Burrage, K.: Fourier spectral methods for fractionalin-space reaction-diffusion equations. BIT Numerical Mathematics pp. 1–18 (2014) 6. Burrage, K., Hale, N., Kay, D.: An efficient implicit fem scheme for fractional-inspace reaction-diffusion equations. SIAM Journal on Scientific Computing 34(4), A2145–A2172 (2012) 7. Higham, N.J.: Functions of matrices: theory and computation. SIAM, Philadelphia (2008) 8. Ilic, M., Liu, F., Turner, I., Anh, V.: Numerical approximation of a fractional-inspace diffusion equation, I. Fractional Calculus and Applied Analysis 8(3), 323–341 (2005) 9. Ilic, M., Liu, F., Turner, I., Anh, V.: Numerical approximation of a fractional-inspace diffusion equation. II with nonhomogeneous boundary conditions. Fractional Calculus and applied analysis 9(4), 333–349 (2006) 10. Ili´c, M., Turner, I.W., Anh, V.: A numerical solution using an adaptively preconditioned lanczos method for a class of linear systems related with the fractional poisson equation. International Journal of Stochastic Analysis 2008, 26 pages (2008) 11. Kilbas, A.A., Srivastava, H.M., Trujillo, J.J.: Theory and Applications of Fractional Differential Equations. North-Holland mathematics studies, Elsevier, Amsterdam (2006) 12. Knabner, P., Angermann, L.: Numerical methods for elliptic and parabolic partial differential equations. Springer Verlag, New York (2003) 13. Logg, A., Mardal, K.A., Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer, Berlin (2012) 14. Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Fitted Numerical Methods For Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions. World Scientific, New Jersey 15. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations. Springer-Verlag, Berlin (1994) 16. Rognes, M.E., Logg, A.: Automated goal-oriented error control i: Stationary variational problems. SIAM Journal on Scientific Computing 35(3), C173–C193 (2013) 17. Roos, H.G., Stynes, M., Tobiska, L.: Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Springer, Berlin (2008) 18. Thom´ee, V.: Galerkin finite element methods for parabolic problems. Springer Verlag, Berlin (2006) 19. Vabishchevich, P.N.: Numerically solving an equation for fractional powers of elliptic operators. Journal of Computational Physics 282(1), 289–302 (2015) 20. Yagi, A.: Abstract Parabolic Evolution Equations and Their Applications. Springer, Berlin (2009)
arXiv:1604.04427v1 [cs.NA] 15 Apr 2016
Petr N. Vabishchevich1,2 1
2
Nuclear Safety Institute, 52, B. Tulskaya, 115191 Moscow, Russia North-Eastern Federal University, 58, Belinskogo, 677000 Yakutsk, Russia
Abstract. A boundary value problem for a fractional power 0 < ε < 1 of the second-order elliptic operator is considered. The boundary value problem is singularly perturbed when ε → 0. It is solved numerically using a time-dependent problem for a pseudo-parabolic equation. For the auxiliary Cauchy problem, the standard two-level schemes with weights are applied. The numerical results are presented for a model two-dimensional boundary value problem with a fractional power of an elliptic operator. Our work focuses on the solution of the boundary value problem with 0 < ε 1.
1
Introduction
Non-local applied mathematical models based on the use of fractional derivatives in time and space are actively discussed in the literature [2, 11]. 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 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 A: Au = −4u, x ∈ Ω. We seek the solution of the problem for the equation with the fractional power of an elliptic operator: Aε u = f, with 0 < ε < 1 for a given f (x), x ∈ Ω. To solve problems with the fractional power of an elliptic operator, we can apply finite volume or finite element methods oriented to using arbitrary domains discretized by irregular computational grids [12, 15]. The computational realization is associated with the implementation of the matrix function-vector multiplication. For such problems, different approaches [7] 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 [10]. 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 [6]. 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 [5, 8, 9]. Unfortunately, all these approaches demonstrates too high computational complexity for multidimensional problems. We have proposed [19] 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-state solution. This computational algorithm for solving equations with fractional powers of operators is promising when considering transient problems. The boundary value problem for the fractional power of an elliptic operator is singularly perturbed when ε → 0. To solve it numerically, we focus on numerical methods that are designed for classical elliptic problems of convection-diffusionreaction [14, 17]. In particular, the main features are taken into account via using locally refining grids. The standard strategy of goal-oriented error control for conforming finite element discretizations [1, 3] is applied.
2
Problem formulation
In a bounded polygonal domain Ω ⊂ Rm , 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 Au = −div(k(x)grad u)
(1)
with coefficient 0 < k1 ≤ k(x) ≤ k2 . The operator A is defined on the set of functions u(x) that satisfy on the boundary ∂Ω the following conditions: u(x) = 0,
x ∈ ∂Ω.
(2)
In the Hilbert space H = L2 (Ω), we define the scalar product and norm in the standard way: Z (u, v) = u(x)v(x)dx, kuk = (u, u)1/2 . Ω
For the spectral problem Aϕk = λk ϕk , ϕk (x) = 0,
x ∈ Ω, x ∈ ∂Ω,
we have λ1 ≤ λ2 ≤ ...,
and the eigenfunctions ϕk , kϕk k = 1, k = 1, 2, ... form a basis in L2 (Ω). Therefore, ∞ X u= (u, ϕk )ϕk . k=1
Let the operator A be defined in the following domain: D(A) = {u | u(x) ∈ L2 (Ω),
∞ X
|(u, ϕk )|2 λk < ∞}.
k=0
Under these conditions the operator A is self-adjoint and positive defined: A = A∗ ≥ δI,
δ > 0,
(3)
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 fractional power of the operator A ∞ X (u, ϕk )λεk ϕk . Aε u = k=0
We seek the solution of the problem with the fractional power of the operator A. The solution u(x) satisfies the equation Aε u = f,
(4)
with 0 < ε < 1 for a given f (x), x ∈ Ω. The key issue in the study of the computational algorithm for solving the problem (4) is to establish the stability of the approximate solution with respect to small perturbations of the right-hand side in various norms. In view of (3), the solution of the problem (4) satisfies the a priori estimate kuk ≤ δ −ε kf k,
(5)
which is valid for all 0 < ε < 1. The boundary value problem for the fractional power of the elliptic operator (4) demonstrates a reduced smoothness when ε → 0. For the solution, we have (see, e.g., [20]) the estimate kuk2ε ≤ Ckf k, with 0 ≤ ε < 1/2, is k · k2ε is the norm in H 2ε (Ω). For the limiting solution, we have u0 (x) = f (x), x ∈ Ω. Thus, a singular behavior of the solution of the problem (4) appears with ε → 0 and is governed by the right-hand side f (x).
3
Discretization in space
To solve numerically the problem (4), we employ finite-element approximations in space [4, 18]. For (1) and (2), we define the bilinear form Z a(u, v) = k grad u grad v. Ω
By (3), we have a(u, u) ≥ δkuk2 . Define a subspace of finite elements V h ⊂ H01 (Ω). Let xi , i = 1, 2, ..., Mh be triangulation points for the domain Ω. Define pyramid function χi (x) ⊂ V h , i = 1, 2, ..., Mh , where 1, if i = j, χi (xj ) = 0, if i 6= j. For v ∈ Vh , we have Mh X
v(x) =
vi χi (x),
i=i
where vi = v(xi ), i = 1, 2, ..., Mh . We have defined Lagrangian finite elements of first degree, i.e., based on the piecewise-linear approximation. We will also use Lagrangian finite elements of second degree defined in a similar way. We define the discrete elliptic operator A as (Ay, v) = a(y, v),
∀ y, v ∈ V h .
The fractional power of the operator A is defined similarly to Aε . For the spectral problem ek Aϕ ek = λ we have e1 ≤ λ e2 ≤ ... ≤ λ eM , λ h
kϕ ek k = 1,
k = 1, 2, ..., Mh .
The domain of definition for the operator A is h
D(A) = {y | y ∈ V ,
Mh X
ek < ∞}. |(y, ϕ ek )|2 λ
k=0
The operator A acts on a finite dimensional space V h defined on the domain D(A) and, similarly to (3), A = A∗ ≥ δI,
δ > 0,
e1 . For the fractional power of the operator A, we suppose where δ ≤ λ1 ≤ λ Aε y =
Mh X k=1
eε ϕ (y, ϕ ek )λ k ek .
(6)
For the problem (4), we put into the correspondence the operator equation for w(t) ∈ V h : Aε w = ψ, (7) where ψ = P f with P denoting L2 -projection onto V h . For the solution of the problem (6), (7), we obtain (see (5)) the estimate kwk ≤ δ −ε kψk,
(8)
for all 0 < ε < 1.
4
Singularly perturbed problem for a diffusion-reaction equation
The object of our study is associated with the development of a computational algorithm for approximate solving the singularly perturbed problem (4). After constructing a finite element approximation, we arrive at equation (7). Features of the solution related to a boundary layer are investigated on a model singularly perturbed problem for an equation of diffusion-reaction. The key moment is associated with selecting adaptive computational grids (triangulations). In view of ε Aε = (exp(ln A)) = I + ε ln A + O(ε2 ), we put the problem (7) into the correspondence with solving the equation εAu + u = ψ.
(9)
The equation (9) corresponds to solving the Dirichlet problem (see the condition (2)) for the diffusion-reaction equation − ε div(k(x)grad u) + u = f (x),
x ∈ Ω.
(10)
Basic computational algorithms for the singularly perturbed boundary problem (2), (10) are considered, for example, in [14, 17]. In terms of practical applications, the most interesting approach is based on an adaptation of a computational grid to peculiarities of the problem solution via a posteriori error estimates. Among main approaches, we highlight the strategy of the goal-oriented error control for conforming finite element discretizations [1,3], which is applied to approximate solving boundary value problems for elliptic equations. The strategy of goal-oriented error control is based on choosing a calculated functional. The accuracy of its evaluation is tracked during computations. In our Dirichlet problem for the second-order elliptic equation, the solution is varied drastically near the boundary. So, it seems natural to control the accuracy of calculations for the normal derivatives of the solution (fluxes) across the boundary or a portion of it. Because of this, we put Z G(u) = − εk(x)(grad u · n)d x, ∂Ω
where n is the outward normal to the boundary. An adaptation of a finite element mesh is based on an iterative local refinement of the grid in order to evaluate the goal functional with a given accuracy η on the deriving approximate solution uh , i.e., |G(u) − G(uh )| ≤ η. To conduct our calculations, we used the FEniCS framework (see, e.g., [13]) developed for general engineering and scientific calculations via finite elements. Features of the goal-oriented procedure for local refinement of the computational grid are described in [16] in detain. Here, we consider only a key idea of the adaptation strategy of finite element meshes, which is associated with selecting the goal functional. The model problem (2), (10) is considered with k(x) = 1,
f (x) = (1 − x1 )x22 ,
in the unit square (Ω = (0, 1) × (0, 1)). The threshold of accuracy for calculating the functional G(u) is defined by the value of η = 10−5 . As an initial mesh, there is used the uniform grid obtained via division by 8 intervals in each direction (step 0 — 128 cells). First, Lagrangian finite elements of first order have been used in our calculations. For this case, the improvement of the goal functional during the iterative procedure of adaptation is illustrated by the data presented in Table 1. Table 2 demonstrates values of the goal functional G(uh ) calculated on the final computational grid, the number of vertices of this final grid and the number of adaptation steps for solving the problem at various values of the small parameter ε. These numerical results demonstrate the efficiency of the proposed strategy for goal-oriented error control for conforming finite element discretizations applied to approximate solving singular perturbed problems of diffusion-reaction (2), (10). Table 1. Calculation of the goal functional during adaptation steps ε Step of adaptation s 0 1 2 3 4 5 6 7 8 9 10
10−1 G(uh ) 0.087608 0.110432 0.116155 0.119766 0.122702 0.125653 0.127950 0.128835 0.129542 0.129940 0.130149
Mh 81 97 140 222 384 694 1235 2179 3841 6540 11040
10−3 G(uh ) Mh 0.0056973 81 0.0107507 98 0.0129506 132 0.0155597 195 0.0175113 305 0.0194985 466 0.0210232 754 0.0221562 1279 0.0229284 2132 0.0234492 3753 0.0237487 6626
10−5 G(uh ) Mh 0.00006643 81 0.00015584 95 0.00023996 120 0.00035644 164 0.00050472 225 0.00068154 349 0.00090839 550 0.00115091 853 0.00137740 1242 0.00161273 1865 0.00181249 2711
Table 2. Adaptation for various values of ε ε 10−1 10−2 10−3 10−4 10−5
Goal functional G(uh ) 0.130396 0.064867 0.024191 0.008061 0.002580
Number of vertices 51868 72297 90170 67476 99003
Number of adaptation steps s 13 14 15 16 18
Next, similar results have been obtained using Lagrangian finite elements of second order. For this case, summary data are presented in Table 3. As expected, the desired accuracy η = 10−5 is reached on adaptive meshes of smaller sizes than in the case of Lagrangian finite elements of first order (see Table 2 for a comparison). Table 3. Adaptation for Lagrangian elements of second order Goal functional G(uh ) 0.130423 0.064884 0.024184 0.008076 0.002574
ε 10−1 10−2 10−3 10−4 10−5
5
Number of vertices 3574 5137 6573 12775 18501
Number of adaptation steps s 7 8 9 11 12
Numerical algorithm for the problem with a fractional power
An approximate solution of the problem (7) is sought as a solution of an auxiliary pseudo-time evolutionary problem [19]. Assume that y(t) = δ ε (t(A − δI) + δI)−ε y(0). Therefore y(1) = δ ε A−ε y(0) and then w = y(1). The function y(t) satisfies the evolutionary equation (tD + δI)
dy + εDy = 0, dt
where D = A − δI.
0 < t ≤ 1,
(11)
By (6), we get D = D∗ > 0.
(12)
We supplement (11) with the initial condition y(0) = δ −ε ψ.
(13)
The solution of equation (7) can be defined as the solution of the Cauchy problem (11)–(13) at the final pseudo-time moment t = 1. For the solution of the problem (11), (13), it is possible to obtain various a priori estimates. The elementary estimate that is consistent with the estimate (8) have the form ky(t)k ≤ ky(0)k. (14) To get (14), multiply scalarly equation (11) by εy + tdy/dt. To solve numerically the problem (11), (13), we use the simplest implicit two-level scheme. Let τ be a step of a uniform grid in time such that y n = y(tn ), tn = nτ , n = 0, 1, ..., N, N τ = 1. Let us approximate equation (11) by the implicit two-level scheme (tσ(n) D + δI)
y n+1 − y n + εDy σ(n) = 0, τ
n = 0, 1, ..., N − 1,
y 0 = δ −ε ψ.
(15) (16)
We use the notation tσ(n) = σtn+1 + (1 − σ)tn ,
y σ(n) = σy n+1 + (1 − σ)y n .
For σ = 0.5, the difference scheme (15), (16) approximates the problem (11), (12) with the second order by τ , whereas for other values of σ, we have only the first order. Theorem 1. For σ ≥ 0.5 the difference scheme (15), (16) is unconditionally stable with respect to the initial data. The approximate solution satisfies the estimate ky n+1 k ≤ ky 0 k, n = 0, 1, ..., N − 1. (17) Proof. Rewrite equation (15) in the following form: y n+1 − y n y n+1 − y n δ + D εy σ(n) + tσ(n) = 0. τ τ Multiplying scalarly it by εy σ(n) + tσ(n) in view of (12), we arrive at
y n+1 − y n , τ
y n+1 − y n σ(n) ,y τ
≤ 0.
We have y
n+1 1 y − yn 1 = σ− τ + (y n+1 + y n ). 2 τ 2
σ(n)
If σ ≥ 0.5, then ky n+1 k ≤ ky n k,
n = 0, 1, ..., N − 1.
Thus, we obtain (17). The key point in approximate solving singularly perturbed boundary value problems is associated with mesh adaptation. In the case of solving the problem (4), we use finite element approximations and proceed to the problem (7) and then formulate the Cauchy problem (11), (13) approximated by the scheme (15), (16). In our case, singularity is associated only with spatial variables. The decomposition of the solution of the problem (11), (13) by eigenfunctions of the operator A results in y(t) =
Nh X
ak (t)ϕ ek .
k=1
For coefficients ak (t), we get ek − δ)t)−ε , ak (t) = (ψ, ϕ ek )(δ + (λ
k = 1, 2, ..., Mh .
Because of this, errors in specifying the initial conditions monotonically decrease for increasing t. A similar behavior demonstrates an approximate solution of the Cauchy problem (11), (13) obtained using the fully implicit scheme with σ = 1 in (15), (16). For the Crank-Nicolson scheme (i.e., σ = 0.5 in (15), (16)), we cannot guarantee a monotone decrease of errors in time, but the error at t = 1 will not be more than at t = 0. The practical significance of such an analysis is that it provides us a simple adaptation strategy for computational grids in solving the problem (11), (13), namely, spatial mesh adaptation is conducted at the first time step of calculations.
6
Solution of a model problem
Below, there are presented some results of numerical solving the problem (7) for small values of ε. A computational algorithm must track a singular behavior of the solution, which is directly related to the singular behavior of the righthand side f (x). Let us consider the problem (2), (10) in the unit square Ω = (0, 1) × (0, 1) with 1 − x2 x1 2 x2 − exp − . k(x) = 1, f (x) = 1 − x1 − exp − µ µ The singularity of the right-hand side (the singularity of a numerical solution of the problem with a fractional power of an elliptic operator) results from existing a boundary layer at low values of µ.
An adaptation of the computational grid is performed during the calculation of the first time step using the two-level scheme (15), (16). For the basic variant, it is assumed that ε = 10−2 , µ = 10−2 , the initial uniform spatial grid contains 8 intervals in each direction and the time step is τ = 10−2 . The parameter δ = 2π 2 corresponds the minimal eigenvalue of the elliptic operator A. Mesh adaptation is carried out taking into account peculiarities of the right-hand side and the goal functional defined in the form Z G(u; t = τ ) = −
k(x)(grad u · n)d x. ∂Ω
Next, the problem (2), (10) is solved using the derived grid in space and the uniform grid in time. Thus, we apply the simplest one-stage starting adaptation of the computational grid for numerical solving the unsteady problem. Lagrangian finite elements of second order are used. For time-stepping, the Crank-Nicolson (σ = 0.5 in (15)) scheme is utilized. The sequence of calculated adaptive grids is shown in Fig. 1. Note that this sequence is weakly dependent on the choice of a time step. The goal functional dynamics for different levels of adaptation is presented in Table 4. The problem is solved with different values of ε.
Table 4. Calculation of the goal functional during adaptation steps ε Step of adaptation s 0 1 2 3 4 5 6 7 8 9 10
10−1 G(uh ; t = τ ) Mh 16.0955 289 24.3875 315 31.1692 399 37.0996 559 42.0854 833 45.7594 1270 48.6087 1849 50.3070 2753 51.2621 4067 51.9766 5965 52.2862 9201
10−2 G(uh ; t = τ ) Mh 21.8932 289 33.7810 315 43.8893 399 52.8249 559 60.4373 837 66.2019 1282 70.4881 1885 73.1491 2778 74.6778 4120 75.8362 5968 76.3402 9235
10−3 G(uh ; t = τ ) Mh 22.5773 289 34.9016 315 45.4226 399 54.7328 559 62.2786 834 68.4363 1264 73.2101 1889 75.9998 2774 77.5762 4125 78.7648 6028 79.2942 9261
Acknowledgements This work was supported by the Russian Foundation for Basic Research (projects 14-01-00785, 15-01-00026).
0 — 128 cells
1 — 140 cells
2 — 180 cells
3 — 256 cells
4 — 388 cells
5 — 599 cells
6 — 886 cells
7 — 1313 cells
Fig. 1. The grid obtained at succesive steps of adaptation
References 1. Ainsworth, M., Oden, J.T.: A Posteriori Error Estimation in Finite Element Analysis. Wiley, New York (2000) 2. Baleanu, D.: Fractional Calculus: Models and Numerical Methods. World Scientific, New York (2012) 3. Bangerth, W., Rannacher, R.: Adaptive Finite Element Methods for Differential Equations. Birkh¨ auser, Basel (2003) 4. Brenner, S.C., Scott, L.R.: The mathematical theory of finite element methods. Springer, New York (2008) 5. Bueno-Orovio, A., Kay, D., Burrage, K.: Fourier spectral methods for fractionalin-space reaction-diffusion equations. BIT Numerical Mathematics pp. 1–18 (2014) 6. Burrage, K., Hale, N., Kay, D.: An efficient implicit fem scheme for fractional-inspace reaction-diffusion equations. SIAM Journal on Scientific Computing 34(4), A2145–A2172 (2012) 7. Higham, N.J.: Functions of matrices: theory and computation. SIAM, Philadelphia (2008) 8. Ilic, M., Liu, F., Turner, I., Anh, V.: Numerical approximation of a fractional-inspace diffusion equation, I. Fractional Calculus and Applied Analysis 8(3), 323–341 (2005) 9. Ilic, M., Liu, F., Turner, I., Anh, V.: Numerical approximation of a fractional-inspace diffusion equation. II with nonhomogeneous boundary conditions. Fractional Calculus and applied analysis 9(4), 333–349 (2006) 10. Ili´c, M., Turner, I.W., Anh, V.: A numerical solution using an adaptively preconditioned lanczos method for a class of linear systems related with the fractional poisson equation. International Journal of Stochastic Analysis 2008, 26 pages (2008) 11. Kilbas, A.A., Srivastava, H.M., Trujillo, J.J.: Theory and Applications of Fractional Differential Equations. North-Holland mathematics studies, Elsevier, Amsterdam (2006) 12. Knabner, P., Angermann, L.: Numerical methods for elliptic and parabolic partial differential equations. Springer Verlag, New York (2003) 13. Logg, A., Mardal, K.A., Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer, Berlin (2012) 14. Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Fitted Numerical Methods For Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions. World Scientific, New Jersey 15. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations. Springer-Verlag, Berlin (1994) 16. Rognes, M.E., Logg, A.: Automated goal-oriented error control i: Stationary variational problems. SIAM Journal on Scientific Computing 35(3), C173–C193 (2013) 17. Roos, H.G., Stynes, M., Tobiska, L.: Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Springer, Berlin (2008) 18. Thom´ee, V.: Galerkin finite element methods for parabolic problems. Springer Verlag, Berlin (2006) 19. Vabishchevich, P.N.: Numerically solving an equation for fractional powers of elliptic operators. Journal of Computational Physics 282(1), 289–302 (2015) 20. Yagi, A.: Abstract Parabolic Evolution Equations and Their Applications. Springer, Berlin (2009)
