Multilevel methods for nonuniformly elliptic operators
Long Chen, Ricardo H. Nochetto, Enrique Otarola, Abner J. Salgado
MMULTILEVEL METHODS FOR NONUNIFORMLY ELLIPTICOPERATORS
LONG CHEN, RICARDO H. NOCHETTO, ENRIQUE OT ´AROLA,AND ABNER J. SALGADO
Abstract.
We develop and analyze multilevel methods for nonuniformly el-liptic operators whose ellipticity holds in a weighted Sobolev space with an A –Muckenhoupt weight. Using the so-called Xu–Zikatanov (XZ) identity, wederive a nearly uniform convergence result, under the assumption that theunderlying mesh is quasi-uniform. We also consider the so-called α -harmonicextension to localize fractional powers of elliptic operators. Motivated by thescheme proposed in [R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDEapproach to fractional diffusion in general domains: a priori error analysis.arXiv:1302.0698, 2013] we present a multilevel method with line smoothersand obtain a nearly uniform convergence result on anisotropic meshes. Nu-merical experiments reveal a competitive performance of our method. Introduction
Recently, a great deal of attention has been paid to the study of fractional andnonlocal operators, both from the point of view of pure mathematical researchas well as motivated by several interesting applications where they constitute afundamental part of the modeling and simulation of complex phenomena that spanvastly different length scales.Fractional and nonlocal operators can be found in a number of applications suchas boundary control problems [30], finance [21, 58], electromagnetic fluids [44], im-age processing [34], materials science [4], optimization [30], porous media flow [26],turbulence [3], peridynamics [28, 29, 48], nonlocal continuum field theories [31] andothers. From this it is evident that the particular type of operator appearing in ap-plications can widely vary and that a unified analysis of their discretizations mightbe well beyond our reach. A more modest, but nevertheless quite ambitious, goal isto develop an analysis and approximation of a model operator that is representativeof a particular class. This is the purpose of our recent research program, in whichwe deal with an important nonlocal operator — fractional powers of the DirichletLaplace operator ( − ∆) s , with s ∈ (0 , Date : September 25, 2018.2010
Mathematics Subject Classification.
Key words and phrases.
Finite elements, weighted Sobolev spaces, Muckenhoupt weights,anisotropic estimates, multilevel methods.LC has been supported by NSF Grant DMS-1115961 and DOE prime award a r X i v : . [ m a t h . NA ] M a r L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO long range or anomalous diffusion is considered as in the flow in porous media [6],or the theory of stochastic processes [7].In previous work [46] we proposed a discretization technique for this operatorand provided an a priori error analysis for it. In this paper, we shall be interested infast multilevel methods for the approximate solution of the discrete problems thatarise from the discretization of the fractional Laplacian. In other words, we shallbe concerned with efficient solution techniques for discretizations of the followingproblem. Let Ω be an open and bounded subset of R n ( n ≥ ∂ Ω.Given s ∈ (0 ,
1) and a smooth enough function f , find u such that(1.1) (cid:40) ( − ∆) s u = f, in Ω ,u = 0 , on ∂ Ω . The fractional Laplacian is a nonlocal operator (see [43, 19, 18]), which is oneof the main difficulties to study and solve problem (1.1). To localize it, Caffarelliand Silvestre showed in [19] that any power of the fractional Laplacian in R n canbe determined as a Dirichlet-to-Neumann operator via an extension problem onthe upper half-space R n +1+ . For a bounded domain Ω, this result has been adaptedin [20, 52], thus obtaining an extension problem which is now posed on the semi-infinite cylinder C = Ω × (0 , ∞ ). This extension is the following mixed boundaryvalue problem:(1.2) div ( y α ∇ U ) = 0 , in C , U = 0 , on ∂ L C ,∂ U ∂ν α = d s f, on Ω × { } , where ∂ L C = ∂ Ω × [0 , ∞ ) denotes the lateral boundary of C , and(1.3) ∂ U ∂ν α = − lim y → + y α U y , is the the so-called conormal exterior derivative of U with ν being the unit outernormal to C at Ω × { } . The parameter α is defined as(1.4) α = 1 − s ∈ ( − , . Finally, d s is a positive normalization constant which depends only on s ; see [19]for details. We will call y the extended variable and the dimension n + 1 in R n +1+ the extended dimension of problem (1.2).The following simple strategy to find the solution of (1.1) has been proposedand analyzed in [46]: given a sufficiently smooth function f we solve (1.2), thusobtaining a function U = U ( x (cid:48) , y ). Setting u : x (cid:48) ∈ Ω (cid:55)→ u ( x (cid:48) ) = U ( x (cid:48) , ∈ R , weobtain the solution of (1.1).For an overview of the existing numerical techniques used to solve problemsinvolving fractional diffusion such as the matrix transference technique and thecontour integral method, we refer to [9, 46]. In addition to [46], two other worksthat deal with the discretization of fractional powers of elliptic operators havesubsequently appeared: the approach given by Bonito and Pasciak in [9] is basedon the integral formulation for self-adjoint operators discussed, for instance, in [8,Chapter 10.4]; the work by del Teso and V´azquez [27] studies the approximation ofthe α -harmonic extension problem via the finite difference method. ULTILEVEL METHODS 3
The main advantage of the algorithm proposed in [46], is that we are solving thelocal problem (1.2) instead of dealing with the nonlocal operator ( − ∆) s of problem(1.1). However, this comes at the expense of incorporating one more dimensionto the problem, thus raising the question of how computationally efficient thisapproach is. A quest for the answer is the motivation for the study of multilevelmethods, since it is known that they are the most efficient techniques for the solutionof discretizations of partial differential equations, see [15, 16, 39, 55]. Multigridmethods for equations of the type (1.2), however, are not very well understood.The purpose of this work is twofold and hinges on the multilevel frameworkdeveloped in [12, 55] and the Xu-Zikatanov identity [57]. First, we show nearlyuniform convergence of a multilevel method for a class of general nonuniformlyelliptic equations on quasi-uniform meshes. Second, we derive an almost uniformconvergence of a multilevel method for the local problem that arises from our PDEapproach to the fractional Laplacian (1.2) on anisotropic meshes [46, 47]. Theformer result assumes that the weight ω in the differential operator belongs tothe so-called Muckenhoupt class A ; see Definition 2.1 for details. A somewhatrelated work is [38] where the authors show a uniform norm equivalence for amultilevel space decomposition under the assumption that the weight belongs tothe smaller class A . Their results and techniques, however, do not apply to oursetting since, simply put, an A -weight is “almost bounded”, which is too restrictive;see Remark 2.2 for details. We make no regularity assumption on the weight ω and show that our estimates solely depend on the A -constant C ,ω . However,our results depend on the number J of levels, and thus logarithmically on themeshsize, which seems unavoidable without further regularity assumptions. Forthe fractional Laplacian, [46] shows that a quasi-uniform mesh cannot yield quasi-optimal error estimates and, consequently, the mesh in the extended dimensionmust be graded towards the bottom of the cylinder thus becoming anisotropic. Weapply line smoothers over vertical lines in the extended domain and prove that thecorresponding multigrid V -cycle converges almost uniformly.We propose an algorithm with complexity O ( M n +1 log M ) for computing anearly optimal approximation of the fractional Laplacian problem (1.1) in R n , where M denotes the number of degrees of freedom in each direction. Notice that usingthe intrinsic integral formulation of the fractional Laplacian [18, 19], a discretizationwould result in a dense matrix with O ( M n ). Special techniques such as fast multi-pole methods [37], the H -matrix methods [41] or wavelet methods [42, 51] might beapplied to reduce the complexity of storage and manipulation of the dense matrixas well as the complexity of solvers.The outline of this paper is as follows. In Section 2, we introduce the notationand functional framework we shall work with. Section 3 contains the salient resultsabout the finite element approximation of nonuniformly elliptic equations includingthe fractional Laplacian on anisotropic meshes. Here we also collect the relevantproperties of a quasi-interpolant which are crucial to obtain the convergence anal-ysis of our multilevel methods. In Section 4, we recall the theory of subspacecorrections [55] and the Xu-Zikatanov identity [57]. We present multigrid algo-rithms for nonuniformly elliptic equations discretized on quasi-uniform meshes inSection 5 and prove their nearly uniform convergence. We adapt the algorithms andanalysis of Section 5 to the fractional Laplacian discretized on anisotropic meshesin Section 6. This requires a line smoother along the extended direction. Finally, L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO to illustrate the performance of our methods and the sharpness of our results, wepresent a series of numerical experiments in Section 7.2.
Notation and preliminaries
Notation.
Throughout this work, Ω is an open, bounded and connected sub-set of R n , with n ≥
1. The boundary of Ω is denoted by ∂ Ω. Unless specifiedotherwise, we will assume that ∂ Ω is Lipschitz. We define the semi-infinite cylinder C = Ω × (0 , ∞ ) , and its lateral boundary ∂ L C = ∂ Ω × [0 , ∞ ) . Given Y >
0, we define the truncated cylinder(2.1) C Y = Ω × (0 , Y ) . The lateral boundary ∂ L C Y is defined accordingly.Throughout our discussion, when dealing with elements defined in R n +1 , we shallneed to distinguish the extended dimension. A vector x ∈ R n +1 , will be denotedby x = ( x , . . . , x n , x n +1 ) = ( x (cid:48) , x n +1 ) = ( x (cid:48) , y ) , with x i ∈ R for i = 1 , . . . , n + 1, x (cid:48) ∈ R n and y ∈ R . The upper half-space in R n +1 will be denoted by R n +1+ = { x = ( x (cid:48) , y ) : x (cid:48) ∈ R n , y ∈ R , y > } . The relation a (cid:46) b indicates that a ≤ Cb , with a constant C that does notdepend on a, b, and the important multilevel discretization parameters J and h J (see Section 4 for their definitions), but it might depend on s and Ω. The value of C might change at each occurrence.If X and Y are topological vector spaces, we write X (cid:44) → Y to denote that X iscontinuously embedded in Y . We denote by X (cid:48) the dual of X . If X is normed, wedenote by (cid:107) · (cid:107) X its norm.2.2. Weighted Sobolev spaces.
In the Caffarelli-Silvestre extension (1.2), theparameter α = 1 − s ∈ ( − , y α is degenerate ( α > α < H but rather theweighted Sobolev space H ( y α , C ), where the weight | y | α belongs to the so-calledMuckenhoupt class A ( R n +1 ); see [33, 45, 53]. For completeness, we recall thedefinition of Muckenhoupt classes. Definition 2.1 (Muckenhoupt class A p ) . Let N ≥ ω ∈ L loc ( R N ) be suchthat ω ( x ) > x ∈ R N . We say that ω ∈ A p ( R N ), 1 < p < ∞ , if there existsa positive constant C p,ω such that(2.2) sup B (cid:18) | B | (cid:90) B ω (cid:19) (cid:18) | B | (cid:90) B ω / (1 − p ) (cid:19) p − = C p,ω < ∞ , where the supremum is taken over all balls B in R N and | B | denotes the Lebesguemeasure of B . In addition, we define A ∞ ( R N ) = (cid:91) p> A p ( R N ) , and A ( R N ) = (cid:92) p> A p ( R N ) . ULTILEVEL METHODS 5 If ω belongs to the Muckenhoupt class A p ( R N ), we say that ω is an A p -weight,and we call the constant C p,ω in (2.2) the A p -constant of ω . Remark A -class) . A useful characterization of the A -Muckenhoupt class is given in [49]: ω ∈ A ( R N ) if and only if(2.3) sup B (cid:107) ω − (cid:107) L ∞ ( B ) | B | (cid:90) B ω = C ,ω < ∞ . Since α ∈ ( − , | y | α ∈ A ( R n +1 ) but | y | α / ∈ A ( R n +1 ).From the A p -condition and H¨older’s inequality follows that an A p -weight satisfiesthe so-called strong doubling property . The proof of this fact is standard; see [53,Proposition 1.2.7] for more details. Proposition 2.1 (strong doubling property) . Let ω ∈ A p ( R n ) with < p < ∞ and let E ⊂ R N be a measurable subset of a ball B ⊂ R N . Then (2.4) ω ( B ) ≤ C p,ω (cid:18) | B || E | (cid:19) p ω ( E ) . For a weight in the Muckenhoupt class A p we define weighted L p spaces asfollows. Definition 2.3 (weighted Lebesgue spaces) . Let ω ∈ A p , and let D ⊂ R N be anopen and bounded domain. For 1 < p < ∞ , we define the weighted Lebesgue space L p ( ω, D ) as the set of measurable functions u on D for which (cid:107) u (cid:107) L p ( ω,D ) = (cid:18)(cid:90) D | u | p ω dx (cid:19) /p < ∞ . Based on the fact that L p ( ω, D ) (cid:44) → L loc ( D ) (cf. [47, Proposition 2.3]), it makessense to talk about weak derivatives of functions in L p ( ω, D ). We define weightedSobolev spaces as follows. Definition 2.4 (weighted Sobolev spaces) . Let D ⊂ R N be an open and boundeddomain, ω ∈ A p with 1 < p < ∞ and m ∈ N . The weighted Sobolev space W mp ( ω, D ) is the space of functions u ∈ L p ( ω, D ) such that for any multi-index κ with | κ | ≤ m , the weak derivatives D κ u ∈ L p ( ω, D ). We endow W mp ( ω, D ) withthe following seminorm and norm | u | W mp ( ω,D ) = (cid:88) | κ | = m (cid:107) D κ u (cid:107) pL p ( ω,D ) /p , (cid:107) u (cid:107) W mp ( ω,D ) = (cid:88) j ≤ m | u | pW jp ( ω,D ) /p , respectively. We also define ◦ W mp ( ω, D ) as the closure of C ∞ ( D ) in W mp ( ω, D ).Owing to the fact that ω ∈ A p , most of the properties of classical Sobolev spaceshave a weighted counterpart; see [33, 35, 53]. In particular we have the followingresult (c.f. [53, Proposition 2.1.2, Corollary 2.1.6] and [35, Theorem 1]). Proposition 2.2 (properties of weighted Sobolev spaces) . Let D ⊂ R N be an openand bounded domain, < p < ∞ , ω ∈ A p ( R N ) and m ∈ N . The spaces W mp ( ω, D ) and ◦ W mp ( ω, D ) are complete and W mp ( ω, D ) ∩ C ∞ ( D ) is dense in W mp ( ω, D ) . L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
The Caffarelli-Silvestre extension problem.
Here we explore problem(1.2); we refer the reader to [19, 18, 46] for details. Since problem (1.2) is posed onthe unbounded domain C , it cannot be directly approximated with finite element-like techniques. However, as [46, Proposition 3.1] shows, the solution U decaysexponentially in y so that, by truncating the cylinder C to C Y and setting a vanishingDirichlet boundary condition on the upper boundary y = Y , we only incur in anexponentially small error in terms of Y [46, Theorem 3.5].Define ◦ H L ( y α , C Y ) = (cid:8) v ∈ H ( y α , C Y ) : v = 0 on ∂ L C Y ∪ Ω × { Y } (cid:9) , where α = 1 − s . Since y α ∈ A ( R n ), Proposition 2.2 shows that ◦ H L ( y α , C Y ) is aHilbert space. We also define(2.5) H s (Ω) = H s (Ω) , s ∈ (0 , ) ,H / (Ω) , s = ,H s (Ω) , s ∈ ( , , which is the natural space for the solution u of problem (1.1); let H − s (Ω) be thedual of H s (Ω). As [46, Proposition 2.5] shows, the trace operator ◦ H L ( y α , C Y ) (cid:51) w (cid:55)→ tr Ω w ∈ H s (Ω) , is well defined. Problem (1.2) then reads: find v ∈ ◦ H L ( y α , C Y ) such that(2.6) (cid:90) C Y y α ∇ v · ∇ φ = d s (cid:104) f, tr Ω φ (cid:105) H − s (Ω) × H s (Ω) , ∀ φ ∈ ◦ H L ( y α , C Y ) , where, (cid:104)· , ·(cid:105) H − s (Ω) × H s (Ω) denotes the duality pairing between H s (Ω) and H − s (Ω).Finally, we recall the exponential convergence result [46, Theorem 3.5]: (cid:107)∇ ( U − v ) (cid:107) L ( C ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H − s (Ω) , where λ denotes the first eigenvalue of the Dirichlet Laplace operator and Y is thetruncation parameter.3. Finite element discretization of nonuniformly elliptic equations
Let D be an open and bounded subset of R N ( N ≥
1) with boundary ∂D and let f ∈ L ( ω − , D ). In this section, we focus on the study of a finite element methodfor the following nonuniformly elliptic boundary value problem: find u ∈ H ( ω, D )that solves(3.1) (cid:40) − div( A ( x ) ∇ u ) = f, in D,u = 0 , on ∂D, where A : D → R N × N is symmetric and satisfies the following nonuniform ellipticitycondition(3.2) ω ( x ) | ξ | (cid:46) ξ (cid:124) A ( x ) ξ (cid:46) ω ( x ) | ξ | , ∀ ξ ∈ R N , a.e. x ∈ D. The function ω belongs to the Muckenhoupt class A , which is defined by (2.1).Examples of nonuniformly elliptic equations are the harmonic extension problem re-lated with the fractional Laplace operator [18, 19], elliptic PDEs in an axisymmetricthree dimensional domain with axisymmetric data [5, 36], and equations modelingthe motion of particles in a central potential field in quantum mechanics [2]. ULTILEVEL METHODS 7
Nonuniformly elliptic elliptic equations of the type (3.1)-(3.2) have been studiedin [33]. Given f ∈ L ( ω − , D ), there exists a unique solution u ∈ H ( ω, D ) [33,Theorem 2.2]. Notice that by taking the weight ω to be y α we see that the un-derlying differential operator in (2.6) is a particular instance of − div( A ( x ) ∇ u ) in(3.1).We define the bilinear form(3.3) a ( u, v ) = (cid:90) D A∇ u · ∇ v d x, which is clearly continuous and coercive in H ( ω, D ). Then, a weak formulation ofproblem (3.1) reads: find u ∈ H ( ω, D ) such that(3.4) a ( u, v ) = (cid:90) D f v d x, ∀ v ∈ H ( ω, D ) . Finite element approximation on quasi-uniform meshes.
Let us de-scribe the construction of the underlying finite element spaces. To avoid technicaldifficulties, we assume D to be a polyhedral domain. Let T = { T } be a mesh of D into elements T (simplices or cubes) such that¯ D = (cid:91) T ∈ T T, | D | = (cid:88) T ∈ T | T | . The partition T is assumed to be conforming or compatible, i.e., the intersectionof any two cells T and T (cid:48) in T is either empty or a common lower dimensionalelement. We denote by T the collection of all conforming meshes. We say that T is shape regular if there exists a constant σ > T ∈ T , (3.5) max { σ T : T ∈ T } ≤ σ, where σ T := h T /ρ T is the shape coefficient of T . For simplicial elements, h T =diam( T ) and ρ T is the diameter of the largest sphere inscribed in T [17, 25]. Forthe definition of h T and ρ T in the case of n -rectangles, we refer to [25].We assume that the collection of meshes T is conforming and satisfies the reg-ularity assumption (3.5), which says that the element shape does not degeneratewith refinement. A refinement method generating meshes satisfying the shape reg-ular condition (3.5) will be called isotropic refinement . A particular instance of anisotropic refinement is the so called quasi-uniform refinement. We recall that T isquasi-uniform if it is shape regular and for all T ∈ T we havemax { h T : T ∈ T } (cid:46) min { h T : T ∈ T } , where the hidden constant is independent of T . In this case, all the elements onthe same refinement level are of comparable size. We define h T = max T ∈ T h T .Given a mesh T ∈ T , we define the finite element space of continuous piecewisepolynomials of degree one(3.6) V ( T ) = (cid:8) W ∈ C ( ¯ D ) : W | T ∈ P ( T ) ∀ T ∈ T , W | ∂ Ω = 0 (cid:9) , where for a simplicial element T , P ( T ) corresponds to the space of polynomialsof total degree at most one, i.e., P ( T ), and for n -rectangles, P ( T ) stands for thespace of polynomials of degree at most one in each variable, i.e., Q ( T ). L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
The finite element approximation of u , solution of problem (3.1), is defined asthe unique discrete function U T ∈ V ( T ) satisfying(3.7) a ( U T , W ) = (cid:90) D f W, ∀ W ∈ V ( T ) . Quasi-interpolation operator.
Let us recall the main properties of thequasi-interpolation operator Π T introduced and analyzed in [47]. This operatoris based on local averages over stars, and then it is well defined for functions in L p ( ω, D ). We summarize its construction and its approximation properties as fol-lows; see [47] for details.Given a mesh T ∈ T and T ∈ T , we denote by N ( T ) the set of nodes of T . Weset N ( T ) := ∪ T ∈ T N ( T ) and ◦ N ( T ) := N ( T ) ∩ D. Then, any discrete function W ∈ V ( T ) is characterized by its nodal values on the set ◦ N ( T ). Moreover, thefunctions φ v ∈ V ( T ), v ∈ ◦ N ( T ), such that φ v ( w ) = δ vw for all w ∈ N ( T ) are thecanonical basis of V ( T ), and W = (cid:88) v ∈ ◦ N ( T ) W ( v ) φ v . Given a vertex v ∈ N ( T ), we define the star or patch around v as S v = ∪ T (cid:51) v T, and for T ∈ T we define its patch as S T = ∪ v ∈ T S v . For each vertex v ∈ N ( T ), wedefine h v = min { h T : v ∈ T } .Let ψ ∈ C ∞ ( R N ) be such that (cid:82) ψ = 1 and supp ψ ⊂ B , where B denotes theball in R N of radius r centered at zero with r ≤ /σ , with σ defined by (3.5). For v ∈ ◦ N ( T ), we define the rescaled smooth function ψ v ( x ) = 1 h N v ψ (cid:18) v − xh v (cid:19) . Given a smooth function v , we denote by P v ( x, z ) the Taylor polynomial ofdegree one of the function v in the variable z about the point x , i.e., P v ( x, z ) = v ( x ) + ∇ v ( x ) · ( z − x ) . Then, given v ∈ ◦ N ( T ) and a function v ∈ W p ( ω, D ), we define the correspondingaveraged Taylor polynomial of first degree of v about the vertex v as(3.8) Q v v ( z ) = (cid:90) P v ( x, z ) ψ v ( x ) d x. Since supp ψ v ⊂ S v , the integral appearing in (3.8) can be written over S v . More-over, integration by parts shows that Q v is well defined for functions in L ( D );see [8, Proposition 4.1.12]. Consequently, [46, Proposition 2.3] implies that Q v isalso well defined for functions in L p ( ω, D ) with ω ∈ A p ( R N ).Given ω ∈ A p ( R N ) and v ∈ L p ( ω, D ), we define the quasi-interpolant Π T v asthe unique function Π T v ∈ V ( T ) that satisfies Π T v ( v ) = Q v ( v ) if v ∈ N ( T ), andΠ T v ( v ) = 0 if v ∈ N ( T ) ∩ ∂ Ω, i.e.,Π T v ( v ) = (cid:88) v ∈ ◦ N ( T ) Q v ( v ) φ v . For this operator, [47, Section 5] proves stability and interpolation error esti-mates in the weighted L p -norm and W p -seminorm. We recall these results forcompleteness. ULTILEVEL METHODS 9
Proposition 3.1 (weighted stability and local error estimate I) . Let T ∈ T , ω ∈ A p ( R N ) and v ∈ L p ( ω, S T ) . Then, we have the following local stability bound (3.9) (cid:107) Π T v (cid:107) L p ( ω,T ) (cid:46) (cid:107) v (cid:107) L p ( ω,S T ) . If, in addition, v ∈ W p ( ω, S T ) , then we have the local interpolation error estimate (3.10) (cid:107) v − Π T v (cid:107) L p ( ω,T ) (cid:46) h v (cid:107)∇ v (cid:107) L p ( ω,S T ) . The hidden constants in both inequalities depend only on C p,ω , ψ and σ . Proposition 3.2 (weighted stability and local error estimate II) . Let T ∈ T , ω ∈ A p ( R N ) and v ∈ W p ( ω, S T ) . Then, we have the following local stability bound (3.11) (cid:107)∇ Π T v (cid:107) L p ( ω,T ) (cid:46) (cid:107)∇ v (cid:107) L p ( ω,S T ) . If, in addition, v ∈ W p ( ω, S T ) , then (3.12) (cid:107)∇ ( v − Π T v ) (cid:107) L p ( ω,T ) (cid:46) h v (cid:107) D v (cid:107) L p ( ω,S T ) . The hidden constants in both inequalities depend only on C p,ω , ψ and σ . Finite element approximation on anisotropic meshes.
Let us now focusour attention on the finite element discretization of problem (1.2). To do so, wemust first study the regularity of its solution U . An error estimate for v , solution of(2.6), depends on the regularity of U as well [46, § U is much worse in the extended direction as the following estimates from [46,Theorem 2.7] reveal (cid:107) ∆ x (cid:48) U (cid:107) L ( y α , C ) + (cid:107) ∂ y ∇ x (cid:48) U (cid:107) L ( y α , C ) (cid:46) (cid:107) f (cid:107) H − s (Ω) , (3.13) (cid:107) U yy (cid:107) L ( y β , C ) (cid:46) (cid:107) f (cid:107) L (Ω) , (3.14)where β > α + 1. This suggests that graded meshes in the extended variable y play a fundamental role.Estimates (3.13)-(3.14) motivate the construction of a mesh over C Y with cellsof the form T = K × I , where K ⊂ R n is an element that is isoparametricallyequivalent either to the unit cube [0 , n or the unit simplex in R n and I ⊂ R is aninterval. To be precise, let T Ω = { T } be a conforming and shape regular mesh ofΩ. In order to obtain a global regularity assumption for T Y , we assume that thereis a constant σ Y such that if T = K × I and T = K × I ∈ T Y have nonemptyintersection, then(3.15) h I h I ≤ σ Y , where h I = | I | . Exploiting the Cartesian structure of the mesh it is possible tohandle anisotropy in the extended variable and obtain estimates of the form (cid:107) v − Π T Y v (cid:107) L ( y α ,T ) (cid:46) h v (cid:48) (cid:107)∇ x (cid:48) v (cid:107) L ( y α ,S T ) + h v (cid:48)(cid:48) (cid:107) ∂ y v (cid:107) L ( y α ,S T ) , (cid:107) ∂ x j ( v − Π T Y v ) (cid:107) L ( y α ,T ) (cid:46) h v (cid:48) (cid:107)∇ x (cid:48) ∂ x j v (cid:107) L ( y α ,S T ) + h v (cid:48)(cid:48) (cid:107) ∂ y ∂ x j v (cid:107) L ( y α ,S T ) , with j = 1 , . . . , n + 1, where h v (cid:48) = min { h K : v (cid:48) ∈ K } , h v (cid:48)(cid:48) = min { h I : v (cid:48)(cid:48) ∈ I } and v is the solution of problem (2.6); see [46, § U yy ≈ y − α − as y ≈
0, we realize that U / ∈ H ( y α , C ) and the secondestimate is not meaningful for j = n + 1. In view of the regularity estimate (3.14)it is necessary to measure the regularity of U yy with a different weight and thuscompensate with a graded mesh in the extended dimension. This makes anisotropicestimates essential. In order to simplify the analysis and implementation of multilevel techniques, weconsider a sequence of nested discretizations. We construct such meshes as follows.First, we introduce a sequence of nested uniform partitions of the unit interval {T k } , with mesh points (cid:98) y l,k , for l = 0 , . . . , M k and k = 0 , . . . , J . Then, we obtain afamily of meshes of the interval [0 , Y ] given by the mesh points(3.16) y l,k = Y (cid:98) y γl,k , l = 0 , . . . , M k , where γ > / (1 − α ). Then, for k = 0 , . . . , J , we consider a quasi-uniform triangu-lation T Ω ,k of the domain Ω and construct the mesh T Y ,k as the tensor product of T Ω ,k and the partition given in (3.16); hence T Y ,k = M k T Ω ,k . Assuming that T Ω ,k ≈ M nk we have T Y ,k ≈ M n +1 k . Finally, since T Ω ,k is shape regular andquasi-uniform, h T Ω ,k ≈ ( T Ω ,k ) − /n . All these considerations allow us to obtainthe following result [46, Theorem 5.4 and Remark 5.5]. Theorem 3.1 (error estimate) . Denote by V T Y ,k ∈ V ( T Y ,k ) the Galerkin approxi-mation of problem (2.6) with first degree tensor product elements. Then, (cid:107)∇ ( U − V T Y ,k ) (cid:107) L ( y α , C ) (cid:46) | log( T Y ,k ) | s ( T Y ,k ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) , where Y ≈ log( T Y ,k ) . We notice that the anisotropic meshes of the cylinder C Y considered above aresemi-structured by construction. They are generated as the tensor product of anunstructured grid T Ω together with the structured mesh T k . Figure 1 shows anexample of this type of meshes in three dimensions. Figure 1.
A three dimensional graded mesh of the cylinder(0 , × (0 , Y ) with 392 degrees of freedom. The mesh is con-structed as a tensor product of a quasi-uniform mesh of (0 , with cardinality 49 and the image of the quasi-uniform partitionof the interval (0 ,
1) with cardinality 8 under the mapping (3.16).Notice that the approximation estimates (3.9)-(3.12) are local and thus validunder the weak shape regularity condition (3.15). Owing to the tensor productstructure of the mesh, we have the following anisotropic error estimate.
ULTILEVEL METHODS 11
Lemma 3.2 (weighted L anisotropic error estimate) . Let v ∈ ◦ H L ( y α , C Y ) be thesolution of problem (2.6) . Then, the quasi-interpolation operator Π T Y satisfies thefollowing error estimate (cid:107) v − Π T Y v (cid:107) L ( y α , C Y ) (cid:46) T − / ( n +1) Y (cid:0) (cid:107)∇ x (cid:48) v (cid:107) L ( y α , C Y ) + (cid:107) ∂ y v (cid:107) L ( y α , C Y ) (cid:1) . Proof.
This is a direct consequence of the results from [47, §
5] together with theCartesian structure of the mesh T Y . (cid:3) A simple application of the mean value theorem yields(3.17) y l +1 ,k − y l,k = Y M γk (cid:0) ( l + 1) γ − l γ (cid:1) ≤ γ Y M k (cid:18) l + 1 M k (cid:19) γ − ≤ γ Y M k , for every l = 0 , . . . , M k −
1, where γ > / (1 − α ) = 3 / (2 s ) according to (3.16). Inother words, since the meshsize of the quasi-uniform mesh T Ω ,k is O ( M − k ), the sizeof the partitions in the extended variable y can be uniformly controlled by h T Ω ,k for k = 0 , . . . , J . However, γ blows up as s ↓ Multilevel space decomposition and multigrid methods
In this section, we present a V -cycle multigrid algorithm based on the method ofsubspace corrections [12, 55], and we present the key identity of Xu and Zikatanov[57] in order to analyze the convergence of the proposed multigrid algorithm.4.1. Multilevel decomposition.
We follow [10, 11] to present a multilevel de-composition of the space V ( T ). Assume that we have an initial conforming mesh T made of simplices or cubes, and a nested sequence of discretizations { T k } Jk =0 where, for k > T k is obtained by uniform refinement of T k − . We then obtain anested sequence, in the sense of trees, of quasi-uniform meshes T ≤ T ≤ · · · ≤ T J = T . Denoting by h k := h T k the meshsize of the mesh T k , we have that h k (cid:104) ρ k forsome ρ ∈ (0 , J (cid:104) | log h J | . Let V k := V ( T k ) denote the correspondingfinite element space over T k defined by (3.6). We thus get a sequence of nestedspaces V ⊂ V ⊂ · · · ⊂ V J = V , and a macro space decomposition V = J (cid:88) k =0 V k . Note the redundant overlapping of the multilevel decomposition above; in partic-ular, the sum is not direct. We now introduce a space micro-decomposition. Westart by defining N k := N ( T k ) = dim V k , i.e., the number of interior vertices of themesh T k . In order to deal with point and line Gauss-Seidel smoothers, we introducethe following sets of indices: For j = 1 , . . . , M k we denote by I k,j a subset of theindex set { , , . . . , N k } , and assume I k,j satisfies M k (cid:91) j =1 I k,j = { , , . . . , N k } . The sets I k,j may overlap, i.e., given 0 < j , j ≤ M k such that j (cid:54) = j , we mayhave I k,j ∩ I k,j (cid:54) = ∅ . This overlap, however, is finite and independent of J and N k .Upon denoting the standard nodal basis of V k by φ k,i , i = 1 , . . . , N k , we define V k,j = span { φ k,i : i ∈ I k,j } and we have the space decomposition(4.1) V = J (cid:88) k =0 M k (cid:88) j =1 V k,j . Multigrid algorithm.
We now describe the multigrid algorithm for the non-uniformly elliptic problem (3.1). We start by introducing several auxiliary opera-tors. For k = 0 , . . . , J , we define the operator A k : V k → V k by( A k v k , w k ) L ( ω,D ) = a ( v k , w k ) , ∀ v k , w k ∈ V k , where the bilinear form a is defined in (3.3). Notice that this operator is symmetricand positive definite with respect to the weighted L -inner product. The projectionoperator P k : V J → V k in the a -inner product is defined by a ( P k v, w k ) = a ( v, w k ) , ∀ w k ∈ V k , and the weighted L -projection Q k : V J → V k is defined by( Q k v, w k ) L ( ω,D ) = ( v, w k ) L ( ω,D ) , w k ∈ V k . We define, analogously, the operators A k,j : V k,j → V k,j , P k,j : V k → V k,j and Q k,j : V k → V k,j . The operator A k,j can be regarded as the restriction of A k to thesubspace V k,j , and its matrix representation, which is the sub-matrix of A k obtainedby deleting the indices i / ∈ I k,j , is symmetric and positive definite. On the otherhand, the operators P k,j and Q k,j denote the projections with respect to the a - andthe weighted L -inner products into V k,j , respectively. We also remark that thematrix representation of the operator Q k,j is the so-called restriction operator, andthe prolongation operator Q Tk,j corresponds to the natural embedding V k,j (cid:44) → V k .The following property, which is of fundamental importance, will be used frequentlyin the paper(4.2) A k,j P k,j = Q k,j A k . With this notation we define a symmetric V -cycle multigrid method as in Algo-rithm 1. When m = 1, it is equivalent to the application of successive subspacecorrections (SSC) to the decomposition (4.1) with exact sub-solvers A − k,j so that the V -cycle multigrid method has a smoother at each level of block Gauss-Seidel type[11, 55]. In particular, if we consider a nodal decomposition I k,j = { j } we obtaina point-wise Gauss-Seidel smoother. On the other hand, if the indices in I k,j aresuch that the corresponding vertices lie on a straight line, we obtain the so-calledline Gauss-Seidel smoother, which will be essential to efficiently solve problem (1.2)with anisotropic elements.4.3. Analysis of the multigrid method.
In order to prove the nearly uniformconvergence of the symmetric V -cycle multigrid method without any assumptions,we rely on the following fundamental identity developed by Xu and Zikatanov [57];see also [23, 24] for alternative proofs. ULTILEVEL METHODS 13 e = MG ( r, k, m ) input : r ∈ V k — residual; k ∈ { , . . . J } — level; m ∈ N — number of smoothing steps. output : e ∈ V k — an approximate solution of A k e = r . if k = 0 then e = A − r ; end // pre-smoothing: m steps u = 0; for l ← to m do v ← u l − ; for j ← to M k do v ← v + A − k,j Q k,j ( r − A k v ); end u l ← v ; end // coarse grid correction u m +1 = u m + MG ( Q k − ( r − A k u m ) , k − , m ); // post-smoothing: m steps for l ← m + 2 to m + 1 do v ← u l − ; for j ← M k to do v ← v + A − k,j Q k,j ( r − A k v ); end u l ← v ; end // output e = u m +1 ; Algorithm 1:
Symmetric V -cycle multigrid method Theorem 4.1 (XZ Identity) . Let V be a Hilbert space with inner product ( · , · ) A and norm (cid:107) · (cid:107) A . Let V j ⊂ V be a closed subspace of V for j = 0 , . . . , J , satisfying V = J (cid:88) j =0 V j . Denote by P j : V → V j the orthogonal projection in the inner product ( · , · ) A onto V j . Then, the following identity holds (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) J (cid:89) j =0 ( I − P j ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) A = 1 −
11 + c , where (4.3) c = sup (cid:107) ν (cid:107) A =1 inf (cid:80) Ji =0 ν i = ν J (cid:88) i =0 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) P i J (cid:88) j = i +1 ν j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) A . The XZ identity given by Theorem 4.1, the properties of the interpolation opera-tor Π T defined in § V -cycle multigridmethod described in Algorithm 1, without resorting to any regularity assumptionson the solution. To see how this is possible we recall the basic ingredients in theanalysis of multilevel methods; see [10, 11, 13, 55] for details. We introduce, for k = 1 , , . . . , J , the operator K k = ( I − A − k, M k Q k, M k A k ) · · · ( I − A − k, Q k, A k )= ( I − P k, M k ) · · · ( I − P k, ) = M k (cid:89) j =1 ( I − P k,j ) , where we used (4.2) to obtain the second equality. With this notation, Algorithm 1can then be recast as a two-layer iterative scheme for the solution of A J u = f ofthe form u (cid:96) +1 = u (cid:96) + B J ( f − A J u (cid:96) ) , where the iterator B J satisfies I − B J A J = ( K mJ ) ∗ · · · ( K m ) ∗ ( I − P ) K m · · · K mJ , with M ∗ denoting the adjoint operator of M with respect to the a –inner product.Notice that I − B J A J is the so-called error transfer operator i.e., u − u (cid:96) +1 = ( I − B J A J ) (cid:0) u − u (cid:96) (cid:1) . Consequently, to show convergence of our scheme we must show that this operatoris a contraction with a contraction factor, ideally, independent of J . Owing to thefact that (cid:107) K mk (cid:107) A ≤ (cid:107) K k (cid:107) mA ≤ (cid:107) K k (cid:107) A , it suffices to consider the case m = 1, where, given an operator S , we denote by (cid:107) S (cid:107) A the operator norm induced by the bilinear form a . Therefore(4.4) (cid:107) I − B J A J (cid:107) A ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) J (cid:89) k =0 M k (cid:89) j =1 ( I − P k,j ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) A , because P is an exact solve and thus a projection, whence ( I − P ) = I − P .Notice that the right hand side of (4.4) is precisely the quantity that the XZ identityprovides a value for. In conclusion, based on Theorem 4.1, to prove the convergenceof the symmetric V -cycle multigrid method described in Algorithm 1, we mustobtain an estimate for the constant c given by (4.3), which will be the content ofthe next two sections.5. Analysis of Multigrid Methods on Quasi-uniform Grids
In this section we consider the V -cycle multigrid method described in Algorithm 1applied to solve the weighted discrete problem (3.7) on quasi-uniform meshes. Weconsider standard pointwise Gauss-Seidel smoothers and prove the convergence ofAlgorithm 1 with a nearly optimal rate up to a factor J (cid:104) | log h J | . Our maincontribution is the extension of the standard multigrid analysis [15, 16, 39, 55] toinclude weights belonging to the Muckenhoput class A ( R N ). An optimal resultfor weights in the A ( R N )-class is derived in [38]. Nevertheless, since our main ULTILEVEL METHODS 15 motivation is the fractional Laplacian, and the weight y α ∈ A ( R N ) \ A ( R N ), weneed to consider the larger class A ( R N ).5.1. Stability of the nodal decomposition in the weighted L -norm. Thefollowing result states that the nodal decomposition is stable in the weighted L -norm or, equivalently, the mass matrix for this inner product is spectrally equivalentto its diagonal. Lemma 5.1 (stability of the nodal decomposition) . Let T ∈ T be a quasi-uniformmesh, and let v ∈ V ( T ) . Then, we have the following norm equivalence (5.1) N ( T ) (cid:88) i =1 (cid:107) v i (cid:107) L ( ω,D ) (cid:46) (cid:107) v (cid:107) L ( ω,D ) (cid:46) N ( T ) (cid:88) i =1 (cid:107) v i (cid:107) L ( ω,D ) , where v = (cid:80) N ( T ) i =1 v i denotes the nodal decomposition for v , and the hidden con-stants in each inequality above only depend on the dimension and the A -constantof the weight ω .Proof. Let (cid:98) T ⊂ R N be a reference element and { (cid:98) φ , . . . , (cid:98) φ N (cid:98) T } be its local shapefunctions, where N (cid:98) T is the number of vertices of (cid:98) T . A standard argument shows (cid:98) c (cid:18)(cid:90) ˆ T (cid:98) ω (cid:19) N (cid:98) T (cid:88) i =1 (cid:98) V i ≤ (cid:107) (cid:98) v (cid:107) L ( (cid:98) ω, (cid:98) T ) ≤ (cid:98) c (cid:18)(cid:90) ˆ T (cid:98) ω (cid:19) N (cid:98) T (cid:88) i =1 (cid:98) V i , where 0 < (cid:98) c ≤ (cid:98) c , (cid:98) v = (cid:80) N (cid:98) T i =1 (cid:98) V i (cid:98) φ i and (cid:98) ω is a weight; see [32, Lemma 9.7]. Now,given T ∈ T , we denote by F T : (cid:98) T → T the mapping such that ˆ v = v ◦ F T . Sincethe A class is invariant under isotropic dilations [47, Proposition 2.1], a scalingargument shows (cid:18)(cid:90) T ω (cid:19) N T (cid:88) i =1 V i (cid:46) (cid:107) v (cid:107) L ( ω,T ) (cid:46) (cid:18)(cid:90) T ω (cid:19) N T (cid:88) i =1 V i . It remains thus to show that (cid:82) T ω ≈ (cid:82) T ωφ i . The fact that 0 ≤ φ i ≤ (cid:90) T ωφ i ≤ (cid:90) T ω. The converse inequality follows from the strong doubling property of ω given inProposition 2.1. In fact, setting E = { x ∈ T : φ i ≥ } ⊂ T , we have (cid:90) T ωφ i ≥ (cid:90) E ωφ i ≥ (cid:90) E ω ≥ C ,ω (cid:18) | E || T | (cid:19) (cid:90) T ω. Finally, notice that the supports of the nodal basis functions { φ i } N ( T ) i =1 have a finiteoverlap which is independent of the refinement level, i.e., for every i = 1 , . . . , N ( T ),the number n ( i ) = { j : supp φ i ∩ supp φ j (cid:54) = ∅} is uniformly bounded. We arriveat (5.1) summing over all the elements T ∈ T . (cid:3) With the aid of the stability of the nodal decomposition, we now show a weightedinverse inequality.
Lemma 5.2 (weighted inverse inequality) . Let T ∈ T be a quasi-uniform mesh,and let T ∈ T and v ∈ V ( T ) . Then, we have the following inverse inequality (5.2) (cid:107)∇ v (cid:107) L ( ω,T ) (cid:46) h − T (cid:107) v (cid:107) L ( ω,T ) . Proof.
Since T is quasi-uniform with meshsize h T , we have |∇ φ i | (cid:46) h − T , and (cid:90) T ω |∇ v | (cid:46) h − T N T (cid:88) i =1 V i (cid:90) T ω, where, we have used the nodal decomposition of v = (cid:80) N T i =1 V i φ i . As in the proof ofLemma 5.1, the strong doubling property of ω yields (cid:90) T ω (cid:46) C ,ω (cid:90) T ωφ i so that we obtain (cid:90) T ω |∇ v | (cid:46) C ,ω h − T N (cid:88) i =1 V i (cid:90) T ωφ i (cid:46) C ,ω h − T (cid:90) T ωv , where, in the last step, we have used (5.1). This concludes the proof. (cid:3) Convergence analysis.
We now present a convergence analysis of Algo-rithm 1 applied to solve the weighted discrete problem (3.7) over quasi-uniformmeshes and with standard pointwise Gauss-Seidel smoothers i.e., M k = N k and I k,j = { j } for j = 1 , . . . N k . The main ingredients in such analysis are the stabilityof the nodal decomposition obtained in Lemma 5.1, the weighted inverse inequalityof Lemma 5.2, and the properties of the quasi-interpolant introduced in Section 3.We follow [54, 56]. Theorem 5.3 (convergence of symmetric V -cycle multigrid) . Algorithm 1 withpoint-wise Gauss-Seidel smoother is convergent with a contraction rate δ ≤ −
11 +
CJ , where C is independent of the meshsize, and it depends on the weight ω only throughthe constant C ,ω defined in (2.2) .Proof. By the XZ identity stated in Theorem 4.1, we only need to estimate(5.3) c = sup (cid:107) v (cid:107) H ω,D ) =1 inf (cid:80) Jk =0 (cid:80) N ki =1 v k,i = v J (cid:88) k =0 N k (cid:88) i =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∇ P k,i (cid:88) ( l,j ) (cid:31) ( k,i ) v l,j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( ω,D ) , where (cid:31) stands for the so called lexicographic ordering , i.e.,( l, j ) (cid:31) ( k, i ) ⇔ (cid:26) l > k,l = k & j > i. We recall that k = 0 , · · · J , j = 1 , . . . , N k and the operator P k,i : V k → V k,i isthe projection with respect to the bilinear form a . For k = 0 , · · · J we denote byΠ T k the quasi-interpolation operator defined in § T k . Next, weintroduce the telescopic multilevel decomposition(5.4) v = J (cid:88) k =0 v k , v k = (Π T k − Π T k − ) v, Π T − v := 0 , ULTILEVEL METHODS 17 along with the nodal decomposition v k = N k (cid:88) i =1 v k,i , for each level k . Consequently, the right hand side of (5.3) can be rewritten byusing the telescopic multilevel decomposition (5.4) as follows: V k,i := (cid:88) ( l,j ) (cid:31) ( k,i ) v l,j = J (cid:88) l = k +1 N k (cid:88) j =1 v l,j + N k (cid:88) j = i +1 v k,j = J (cid:88) l = k +1 v l + N k (cid:88) j = i +1 v k,j = v − Π T k v + N k (cid:88) j = i +1 v k,j . Therefore, we have (cid:107)∇ P k,i V k,i (cid:107) L ( ω,D ) (cid:46) (cid:107)∇ P k,i ( v − Π T k v ) (cid:107) L ( ω,D ) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∇ P k,i N k (cid:88) j = i +1 v k,j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( ω,D ) (cid:46) (cid:107)∇ ( v − Π T k v ) (cid:107) L ( ω,D k,i ) + N k (cid:88) j = i +1 D k,i ∩ D k,j (cid:54) = ∅ (cid:107)∇ v k,j (cid:107) L ( ω,D ) , where D k,i = supp φ k,i . Adding over i = 1 , . . . , N k , and using the finite overlappingproperty of the sets D k,i , yields N k (cid:88) i =1 N k (cid:88) j = i +1 D k,i ∩ D k,j (cid:54) = ∅ (cid:107)∇ v k,j (cid:107) L ( ω,D ) (cid:46) N k (cid:88) i =1 (cid:107)∇ v k,i (cid:107) L ( ω,D ) , whence, the weighted inverse inequality (5.2) gives N k (cid:88) i =1 (cid:107)∇ P k,i V k,i (cid:107) L ( ω,D ) (cid:46) (cid:107)∇ ( v − Π T k v ) (cid:107) L ( ω,D ) + N k (cid:88) i =1 h − k (cid:107) v k,i (cid:107) L ( ω,D ) . We resort to the stability of the operator Π T k , Proposition 3.2, and the stability ofthe micro decomposition, Lemma 5.1, to arrive at N k (cid:88) i =1 (cid:107)∇ P k,i V k,i (cid:107) L ( ω,D ) (cid:46) (cid:107)∇ v (cid:107) L ( ω,D ) + h − k (cid:107) v k (cid:107) L ( ω,D ) . Since v k = (Π T k − Π T k − ) v , we utilize the approximation properties of Π T k , givenin Proposition 3.1, to deduce (cid:107) v k (cid:107) L ( ω,D ) ≤ (cid:107) v − Π T k v (cid:107) L ( ω,D ) + (cid:107) v − Π T k − v (cid:107) L ( ω,D ) (cid:46) h k (cid:107)∇ v (cid:107) L ( ω,D ) . This implies (cid:80) N k i =1 (cid:107)∇ P k,i V k,i (cid:107) L ( ω,D ) (cid:46) (cid:107)∇ v (cid:107) L ( ω,D ) , and adding over k from 0 to J yields c (cid:46) J , which completes the proof. (cid:3) A multigrid method for the fractional Laplace operator onanisotropic meshes
As we explained in § y . This allows us to recover an almost-optimal error estimate for the finite element approximation of problem (1.2) [46,Theorem 5.4]. In fact, finite elements on quasi-uniform meshes have poor approxi-mation properties for small values of the parameter s . The isotropic error estimatesof [46, Theorem 5.1] are not optimal, which makes anisotropic estimates essential.For this reason, in this section we develop a multilevel theory for problem (1.2)having in mind anisotropic partitions in the extended variable y and the multilevelsetting described in Section 4 for the nonuniformly elliptic equation (3.1). We shallobtain nearly uniform convergence of a V -cycle multilevel method for the problem(1.2) without any regularity assumptions. We consider line Gauss-Seidel smoothers.The analysis is an adaptation of the results presented in [54] for anisotropic ellipticequations, and it is again based on the XZ identity [57].6.1. A multigrid algorithm with line smoothers.
As W. Hackbusch rightfullyexplains [40]: “ the multigrid method cannot be understood as a fixed algorithm.Usually, the components of the multigrid iteration should be adapted to the givenproblem, [...] being the smoothing iteration the most delicate part of the multigridprocess ”.The success of multigrid methods for uniformly elliptic operators is due to thefact that the smoothers are effective in reducing the nonsmooth (high frequency)components of the error and the coarse grid corrections are effective in reducing thesmooth (low frequency) components. However, the effectiveness of both strategiesdepends crucially on several factors such as the anisotropy of the mesh. A keyingredient in the design and analysis of a multigrid method on anisotropic meshesis the use of the so called line smoothers; see [1, 14, 40, 50].Intuitively, when solving the α -harmonic extension (1.2) on graded meshes, theapproximation from the coarse grid is dominated by the larger meshsize in the x -direction and thus the coarse grid correction cannot capture the smaller scale inthe y -direction. One possible solution is the use of semi-coarsening, i.e., coarseningonly the y -direction until the meshsizes in both directions are comparable. Anothersolution is the use of line smoothing, i.e., solving sub-problems restricted to onevertical line. We shall use the latter approach which is relatively easy to implementfor tensor-product meshes.Let us describe the decomposition of V J = V ( T Y J ) that we shall use. To do so,we follow the notation of § M k to be the number of interior nodes of T Ω ,k and define, for j = 1 , . . . , M k , the set I k,j as the collection of indices for thevertices that lie on the line { v (cid:48) j } × (0 , Y ) at the level k . The decomposition is thengiven by (4.1). This decomposition is also stable, which allows us to obtain theappropriate anisotropic inverse inequalities; see Lemma 6.1 below.Owing to the nature of the decomposition, the smoother requires the evaluationof A − k,j which corresponds to the action of the operator over a vertical line. Thiscan be efficiently realized since the corresponding matrix is tri-diagonal. Lemma 6.1 (nodal stability and anisotropic inverse inequalities) . Let T Y be agraded tensor product grid, which is quasi-uniform in Ω and graded in the extendedvariable so that (3.16) holds. If v ∈ V ( T Y ) can be decomposed as v = (cid:80) M J j =1 v j , ULTILEVEL METHODS 19 then (6.1) M J (cid:88) j =1 (cid:107) v j (cid:107) L ( y α , C Y ) (cid:46) (cid:107) v (cid:107) L ( y α , C Y ) (cid:46) M J (cid:88) j =1 (cid:107) v j (cid:107) L ( y α , C Y ) . Moreover, we have the following inverse inequalities (6.2) (cid:107)∇ x (cid:48) v (cid:107) L ( y α ,T ) (cid:46) h − K (cid:107) v (cid:107) L ( y α ,T ) , (cid:107) ∂ y v (cid:107) L ( y α ,T ) (cid:46) h − I (cid:107) v (cid:107) L ( y α ,T ) , where T = K × I is a generic element of T Y .Proof. The nodal stability (6.1) follows along the same lines of Lemma 5.1 uponrealizing that the functions v j = v j ( x (cid:48) , y ) are defined on the vertical lines ( v (cid:48) j , y )with y ∈ (0 , Y ) and the index j corresponds to a nodal decomposition in Ω. More-over, noticing that |∇ x (cid:48) φ i | (cid:46) h − K and | ∂ y φ i | (cid:46) h − I , we derive (6.2) inspired inLemma 5.2. (cid:3) We examine Algorithm 1 applied to the decomposition (4.1) with exact sub-solvers on V k,j , i.e., with line smoothers; see [13, § III.12] and [54]. A key observationin favor of subspaces { V k,j } M k j =1 follows. Lemma 6.2 (nodal stability of y -derivatives) . Under the same assumptions ofLemma 6.1 we have (6.3) M J (cid:88) j =1 (cid:107) ∂ y v j (cid:107) L ( y α , C Y ) (cid:46) (cid:107) ∂ y v (cid:107) L ( y α , C Y ) (cid:46) M J (cid:88) j =1 (cid:107) ∂ y v j (cid:107) L ( y α , C Y ) . Proof.
We just proceed as in Lemma (5.1) with v replaced by ∂ y v = (cid:80) M J j =1 ∂ y v j . (cid:3) Exploiting Theorem 4.1, the properties of the quasi-interpolation operator Π T k defined in § V -cycle multigrid method. We follow [54, 56]. Theorem 6.3 (convergence of multigrid methods with line smoothers) . The sym-metric V -cycle multigrid method with line smoothing converges with a contractionrate δ ≤ −
11 +
CJ , where C is independent of the number of degrees of freedom. The constant C de-pends on the weight y α only through the constant C ,y α , and on s like C ≈ γ , where γ is the parameter that defines the graded mesh (3.16) .Proof. We use the XZ identity (4.1) and modify the arguments in the proof ofTheorem 5.3. We introduce the telescopic multilevel decomposition(6.4) v = J (cid:88) k =0 v k , v k = (Π T Y ,k − Π T Y ,k − ) v, Π T Y , − v := 0 , along with the line decomposition v k = M k (cid:88) j =1 v k,j . Following the same arguments developed in the proof of Theorem 5.3, and denoting V k,i = (cid:80) ( l,j ) (cid:31) ( k,i ) v l,j , we arrive at the inequality(6.5) M k (cid:88) i =1 (cid:107)∇ P k,i V k,i (cid:107) L ( y α , C Y ) (cid:46) (cid:107)∇ ( v − Π T Y ,k v ) (cid:107) L ( y α , C Y ) + M k (cid:88) j =1 (cid:107)∇ v k,j (cid:107) L ( y α , C Y ) , where we have used the finite overlapping property of the sets I k,j ; see § T Y ,k stated in (3.11) (see also [46, Theorems 4.7 and 4.8] and [47, Lemma 5.1])yields(6.6) (cid:107)∇ ( v − Π T Y ,k v ) (cid:107) L ( y α , C Y ) (cid:46) (cid:107)∇ v (cid:107) L ( y α , C Y ) . To estimate the second term in (6.5) we begin by noticing that(6.7) M k (cid:88) j =1 (cid:107)∇ v k,j (cid:107) L ( y α , C Y ) = M k (cid:88) j =1 (cid:107)∇ x (cid:48) v k,j (cid:107) L ( y α , C Y ) + M k (cid:88) j =1 (cid:107) ∂ y v k,j (cid:107) L ( y α , C Y ) . The first term is estimated via the first weighted inverse inequality (6.2) and thestability of the nodal decomposition (6.1), that is(6.8) M k (cid:88) j =1 (cid:107)∇ x (cid:48) v k,j (cid:107) L ( y α , C Y ) (cid:46) M k (cid:88) j =1 h (cid:48) k − (cid:107) v k,j (cid:107) L ( y α , C Y ) (cid:46) h (cid:48) k − (cid:107) v k (cid:107) L ( y α , C Y ) , where h (cid:48) k denotes the meshsize in the x (cid:48) direction at level k . The approximationproperty of Π T Y ,k stated in Lemma 3.2 (see also [47, Theorem 5.7]) and the defini-tion of v k yield (cid:107) v k (cid:107) L ( y α , C Y ) ≤ (cid:107) v − Π T Y ,k v (cid:107) L ( y α , C Y ) + (cid:107) v − Π T Y ,k − v (cid:107) L ( y α , C Y ) (cid:46) h (cid:48) k (cid:107)∇ x (cid:48) v (cid:107) L ( y α , C Y ) + h (cid:48)(cid:48) k (cid:107) ∂ y v (cid:107) L ( y α , C Y ) where h (cid:48)(cid:48) k denotes the maximal meshsize in the y direction at level k . Using (3.17)we see that h (cid:48)(cid:48) k (cid:46) γh (cid:48) k , and replacing the estimate above in (6.8), we obtain(6.9) M k (cid:88) j =1 (cid:107)∇ x (cid:48) v k,j (cid:107) L ( y α , C Y ) (cid:46) (cid:107)∇ v (cid:107) L ( y α , C Y ) , which bounds the first term in (6.7). To estimate the second term, we resort toLemma 6.2, namely(6.10) M k (cid:88) j =1 (cid:107) ∂ y v k,j (cid:107) L ( y α , C Y ) (cid:46) (cid:107) ∂ y v k (cid:107) L ( y α , C Y ) . Finally, inequalities (6.9) and (6.10) allow us to conclude M k (cid:88) j =1 (cid:107)∇ v k,j (cid:107) L ( y α , C Y ) (cid:46) (cid:107)∇ v (cid:107) L ( y α , C Y ) , which together with (6.6) yields the desired result after summing over k . (cid:3) Remark s ) . We point out the use of (3.17), which in turnimplies h (cid:48)(cid:48) k (cid:46) γh (cid:48) k , to derive (6.9). This translates into C ≈ γ in Theorem 6.3 and,since γ > / (1 − α ) = 3 / (2 s ), in deterioration of the contraction factor as s ↓ § ULTILEVEL METHODS 21 Numerical Illustrations
In this section, we present numerical experiments to support our theoreticalfindings. We consider two examples:(7.1) n = 1, Ω = (0 , u = sin(3 πx ),(7.2) n = 2, Ω = (0 , , u = sin(2 πx ) sin(2 πx ),and Y = 1. The length Y of the cylinder in the extended direction is fixed, asdiscussed in [46], so that it captures the exponential decay of the solution. All of ouralgorithms are implemented based on the MATLAB c (cid:13) software package i FEM [22].7.1.
Multigrid with line smoothers on graded meshes.
We partition Ω into auniform grid of size h T Ω , and we construct a graded mesh in the extended directionusing the mapping (3.16) with parameter γ = s + 0 . M = h . Some samplemeshes are shown in Figure 2. The mesh points are ordered column-wise so thatthe indices associated to vertical lines are easily accessible. Starting from h T = we obtain a sequence of meshes by halving the meshsize of Ω and applying themapping (3.16) in the extended direction with double number of mesh points. Figure 2.
Graded meshes for the cylinder C Y = (0 , × (0 , h T = . Themeshes in the extended direction are graded according to (3.16)with M = h and γ = s + 0 .
1. The left mesh is for s = 0 . right for s = 0 . V ( T k ) → V ( T k +1 ) for k = 0 , . . . , J − x (cid:48) -direction is obtained by standard averaging, while in theextended direction the weights must be modified to take into account the gradingof the mesh. The restriction matrix is taken as the transpose of the prolongationmatrix.As discussed in Section 6 we must use vertical line smoothers to attain efficiencyof the multigrid method. The tri-diagonal sub-matrix corresponding to one verticalline is inverted exactly by using the built-in direct solver in MATLAB c (cid:13) . Red-blackordering of the indices in the x (cid:48) -direction is used to further improve the efficiency ofthe line smoothers. We perform three pre- and post-smoothing steps, i.e., m = 3. We start with the zero initial guess and use as exit criterion that the (cid:96) -norm ofthe relative residual is smaller than 10 − .Tables 1 and 2 show the number of iterations for the implemented multigridmethod for the one and two dimensional problems, respectively. As we see, themethod converges almost uniformly with respect to the number of degrees of free-dom. Notice that the number of iterations for s = 0 .
15 is significantly larger thanthat for the remaining tested cases. This can be explained by the fact that, asTheorem 6.3 states, the contraction factor depends on γ ≈ s and thus, we observea preasymptotic regime where the number of iterations grows. This is exactly thecase for the one dimensional problem and we would expect a similar behavior inthe two dimensional case. However, since the extended problem is now in three di-mensions, the size of the problems grows rather quickly and thus our computationalresources were not sufficient to deal with the cases h T Ω = and h T Ω = . In § h T Ω DOFs s = 0 . s = 0 . s = 0 . s = 0 .
289 7 6 5 5
Table 1.
Number of iterations for a multigrid method for theone dimensional fractional Laplacian using a line smoother in theextended direction. The mesh in Ω is uniform of size h T Ω . Themesh in the extended direction is graded according to (3.16). h T Ω DOFs s = 0 . s = 0 . s = 0 . s = 0 . Table 2.
Number of iterations for a multigrid method for thetwo dimensional fractional Laplacian using a line smoother in theextended direction. The mesh in Ω is uniform of size h T Ω . Themesh in the extended direction is graded according to (3.16).We also tested a point Gauss-Seidel smoother for the one dimensional case Ω =(0 , h T Ω = 1 /
16, the corresponding V -cycle is not ableto achieve the desired accuracy in 200 iterations.7.2. Multigrid methods on quasi-uniform meshes.
Even though the approx-imation of the Caffarelli-Silvestre extension of the fractional Laplace operator on
ULTILEVEL METHODS 23 quasi-uniform meshes in the extended direction is suboptimal, let us use this prob-lem to illustrate the convergence properties of the multilevel method, developedin Section 5, for general A weights. The setting is the same as in the previoussubsection but we use a point-wise Gauss-Seidel smoother. Tables 3 and 4 show thenumber of iterations with respect to the number of degrees of freedom and s . Wesee that the convergence is almost uniform with respect to the number of unknownsas well as the parameter s ∈ (0 , h T Ω DOFs s = 0 . s = 0 . s = 0 . s = 0 .
289 12 13 13 14
Table 3.
Number of iterations for a multigrid method with point-wise Gauss-Seidel smoothers on uniform meshes for the one dimen-sional fractional Laplacian. h T Ω DOFs s = 0 . s = 0 . s = 0 . s = 0 . Table 4.
Number of iterations for a multigrid method with point-wise Gauss-Seidel smoothers on uniform meshes for the two dimen-sional fractional Laplacian.7.3.
Modified mesh grading.
Examining the proof of Theorem 6.3, we realizethat the critical step (6.9) consists in the application of inequality (3.17), namely h (cid:48)(cid:48) k (cid:46) γh (cid:48) k , which deteriorates as s becomes small because γ > / (1 − α ) = 3 / (2 s ).Numerically, this effect can be seen in Tables 1 and 2 where, for instance, thenumber of iterations needed for s = 0 .
15 is significantly larger than that for allthe other tested values; see the right mesh for s = 0 .
15 in Figure 2. As a result,the contraction rate of Theorem 6.3 becomes 1 − / (1 + CγJ ). Here we explorecomputationally how to overcome this issue. We construct a mesh such that themaximum meshsize in the extended direction is uniformly bounded, with respect to s , by the uniform meshsize in the x (cid:48) -direction without changing the ratio of degreesof freedom in Ω and the extended direction by more than a constant.Let us begin with some heuristic motivation. In order to control the aspectratio h (cid:48)(cid:48) k /h (cid:48) k uniformly on s ∈ (0 , y direction, increasing the number of degrees of freedom of T Y just by a constant. We denote by ˜ T Y the resulting mesh and we notice that V ( T Y ) ⊂ V ( ˜ T Y ). Thus, Galerkin orthogonality implies (cid:107)∇ ( v − V ˜ T Y ) (cid:107) L ( y α , C Y ) = inf (cid:110) (cid:107)∇ ( v − W ) (cid:107) L ( y α , C Y ) : W ∈ V ( ˜ T Y ) (cid:111) ≤ (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) (cid:46) ( T Y ) − n +1 ≈ ( T Y ) − n +1 . We build on this idea through a modification of the mapping function below.Let F : (0 , → (0 , Y ) be an increasing and differentiable function such that F (0) = 0 and F (1) = Y . By mapping a uniform grid of (0 ,
1) via the function F , we can construct a graded mesh with mesh points given by y l = F ( l/M ) for l = 1 , ..., M ; for instance, F ( ξ ) = Y ξ γ yields (3.16). The mean value theoremimplies y l +1 − y l = F (cid:48) ( c l ) M ≤ M max (cid:26) | F (cid:48) ( ξ ) | : ξ ∈ (cid:20) lM , l + 1 M (cid:21)(cid:27) , which shows that the map of (3.16) is not uniformly bounded with respect to s .For this reason, we instead consider the following construction: Let ( ξ (cid:63) , y (cid:63) ) ∈ (0 , , which we will call the transition point , and define the mapping F ( ξ ) = y (cid:63) Y (cid:18) ξξ (cid:63) (cid:19) γ , < ξ ≤ ξ (cid:63) , Y (cid:18) − y (cid:63) − ξ (cid:63) ( ξ − ξ (cid:63) ) + y (cid:63) (cid:19) , ξ (cid:63) < ξ < . Over the interval (0 , ξ (cid:63) ) the mapping F defines the same type of graded meshbut, over ( ξ (cid:63) ,
1) it defines a uniform mesh. Let us now choose the transition pointto obtain a bound on the derivative of F . We have(7.1) F (cid:48) ( ξ ) = γ Y y (cid:63) ξ (cid:63) (cid:18) ξξ (cid:63) (cid:19) γ − , < ξ ≤ ξ (cid:63) , Y − y (cid:63) − ξ (cid:63) , ξ (cid:63) < ξ < , so that F := max ξ ∈ [0 , | F (cid:48) ( ξ ) | = Y max (cid:26) γ y (cid:63) ξ (cid:63) , − y (cid:63) − ξ (cid:63) (cid:27) . Given ξ (cid:63) we choose y (cid:63) to have γ y (cid:63) ξ (cid:63) = − y (cid:63) − ξ (cid:63) , i.e., y (cid:63) = 11 + γ − ξ (cid:63) ξ (cid:63) . this immediately yields F ∈ C ([0 , F = γ Y y (cid:63) ξ (cid:63) = Y γξ (cid:63) + (1 − ξ (cid:63) ) γ ≤ Y − ξ (cid:63) . We can now choose ξ (cid:63) to gain control of F . For instance, ξ (cid:63) = 0 . F ≤ ξ (cid:63) = 0 .
75 that
F ≤
4. In the experiments presented below we choose ξ (cid:63) = 0 .
75. The theory presented in § ULTILEVEL METHODS 25
Figure 3.
Graded meshes for the extended domain C Y = (0 , × (0 , h T Ω = and s = 0 . Left : The grading is according to(3.16).
Right : The grading is given by the map (7.1).
Figure 4.
Graded meshes for the extended domain C Y = (0 , × (0 , h T Ω = and s = 0 . Left : The grading is according to(3.16).
Right : The grading is given by the map (7.1).so they are also capable of capturing the singular behavior of the solution U . How-ever, near the top part, the aspect ratio is uniformly controlled by a factor 4. Themodified mesh is only applied for γ >
4. For s = 0 .
3, 0 . .
8, no modificationis needed in the original mesh.Upon constructing a mesh with this modification, we can develop a V -cyclemultigrid solver with vertical line smoothers. Comparisons of this approach withthe setting of § V -cycle multigrid with vertical line smoothers. For the original graded meshes,there is a preasymptotic regime where the number of iterations increases faster thanlog J . The modification of the mesh proposed in (7.1) allows us to obtain an almostuniform number of iterations for all problem sizes without sacrificing the accuracyof the method. This is also evidenced by the computational time required to solvea problem with a fixed number of degrees of freedom. h T Ω DOFs I(o) I(m) E(o) E(m) CPU(o) CPU(m)
289 7 7 0.1556 0.1739 0.0209 0.0554
Table 5.
Comparison of the multilevel solver with vertical linesmoother over two graded meshes for the one dimensional frac-tional Laplacian, s = 0 . Legend : The original mesh, given by(3.16) is denoted by o, whereas the modification proposed in (7.1)is denoted by m. I – number of iterations, E – error in the energynorm, CPU – cpu time ( s ). h T Ω DOFs I(o) I(m) E(o) E(m) CPU(o) CPU(m)
Table 6.
Comparison of the multilevel solver with vertical linesmoother over two graded meshes for the two dimensional frac-tional Laplacian, s = 0 . Legend : The original mesh, given by(3.16) is denoted by o, whereas the modification proposed in (7.1)is denoted by m. I – number of iterations, E – error in the energynorm, CPU – cpu time ( s ). References [1] T. Apel and J. Sch¨oberl. Multigrid methods for anisotropic edge refinement.
SIAM J. Numer.Anal. , 40(5):1993–2006 (electronic), 2002.[2] D. Arroyo, A. Bespalov, and N. Heuer. On the finite element method for elliptic problemswith degenerate and singular coefficients.
Math. Comp. , 76(258):509–537 (electronic), 2007.[3] O.G. Bakunin.
Turbulence and diffusion . Springer Series in Synergetics. Springer-Verlag,Berlin, 2008. Scaling versus equations.[4] P.W. Bates. On some nonlocal evolution equations arising in materials science. In
Nonlineardynamics and evolution equations , volume 48 of
Fields Inst. Commun. , pages 13–52. Amer.Math. Soc., Providence, RI, 2006.[5] Z. Belhachmi, Ch. Bernardi, and S. Deparis. Weighted Cl´ement operator and applicationto the finite element discretization of the axisymmetric Stokes problem.
Numer. Math. ,105(2):217–247, 2006.[6] D.A. Benson, S.W. Wheatcraft, and M.M. Meerschaert. Application of a fractional advection-dispersion equation.
Water Resources Res. , 36:91–109, 2000.[7] J. Bertoin.
L´evy processes , volume 121 of
Cambridge Tracts in Mathematics . CambridgeUniversity Press, Cambridge, 1996.[8] M.ˇS. Birman and M.Z. Solomjak.
Spektralnaya teoriya samosopryazhennykh operatorov vgilbertovom prostranstve . Leningrad. Univ., Leningrad, 1980.[9] A. Bonito and J.E. Pasciak. Numerical approximation of fractional powers of elliptic opera-tors. Preprint, 2013.
ULTILEVEL METHODS 27 [10] J.H. Bramble and J.E. Pasciak. New convergence estimates for multigrid algorithms.
Math.Comp. , 49(180):311–329, 1987.[11] J.H. Bramble, J.E. Pasciak, J.P. Wang, and J. Xu. Convergence estimates for multigridalgorithms without regularity assumptions.
Math. Comp. , 57(195):23–45, 1991.[12] J.H. Bramble, J.E. Pasciak, J.P. Wang, and J. Xu. Convergence estimates for product iterativemethods with applications to domain decomposition.
Math. Comp. , 57(195):1–21, 1991.[13] J.H. Bramble and X. Zhang. The analysis of multigrid methods. In
Handbook of numericalanalysis, Vol. VII , Handb. Numer. Anal., VII, pages 173–415. North-Holland, Amsterdam,2000.[14] J.H. Bramble and X. Zhang. Uniform convergence of the multigrid V -cycle for an anisotropicproblem. Math. Comp. , 70(234):453–470, 2001.[15] A. Brandt. Multi-level adaptive solutions to boundary-value problems.
Math. Comp. , 31:333–390, 1977.[16] A. Brandt.
Multigrid techniques: 1984 guide with applications to fluid dynamics . Ges. f¨urMathematik u. Datenverarbeitung, 1984.[17] S.C. Brenner and L.R. Scott.
The mathematical theory of finite element methods , volume 15of
Texts in Applied Mathematics . Springer, New York, third edition, 2008.[18] X. Cabr´e and Y. Sire. Nonlinear equations for fractional Laplacians ii: Existence, uniquenessand qualitative properties of solutions. arXiv:1111.0796v1, 2011.[19] L. Caffarelli and L. Silvestre. An extension problem related to the fractional Laplacian.
Comm. Partial Differential Equations , 32(7-9):1245–1260, 2007.[20] A. Capella, J. D´avila, L. Dupaigne, and Y. Sire. Regularity of radial extremal solutions forsome non-local semilinear equations.
Comm. Partial Differential Equations , 36(8):1353–1384,2011.[21] P. Carr, H. Geman, D.B. Madan, and M. Yor. The fine structure of asset returns: An empiricalinvestigation.
Journal of Business , 75:305–332, 2002.[22] L. Chen. i FEM: An integrated finite element methods package in matlab. Technical report,University of California at Irvine, 2009.[23] L. Chen. Deriving the X-Z Identity from Auxiliary Space Method. In Yunqing Huang, RalfKornhuber, Olof Widlund, and Jinchao Xu, editors,
Domain Decomposition Methods in Sci-ence and Engineering XIX , pages 309–316. Springer Berlin Heidelberg, 2010.[24] D. Cho, J. Xu, and L. Zikatanov. New estimates for the rate of convergence of the methodof subspace corrections.
Numer. Math. Theory Methods Appl. , 1(1):44–56, 2008.[25] P.G. Ciarlet.
The finite element method for elliptic problems , volume 40 of
Classics in AppliedMathematics . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2002.[26] J. Cushman and T. Glinn. Nonlocal dispersion in media with continuously evolving scales ofheterogeneity.
Trans. Porous Media , 13:123–138, 1993.[27] F. del Teso and J.L. V´azquez. Finite difference method for a general fractional porous mediumequation. arXiv:1307.2474, 2013.[28] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou. Analysis and approximation of nonlocaldiffusion problems with volume constraints.
SIAM Rev. , 54(4):667–696, 2012.[29] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou. Analysis of the volume-constrainedperidynamic Navier equation of linear elasticity.
J. Elasticity , 113(2):193–217, 2013.[30] G. Duvaut and J.-L. Lions.
Inequalities in mechanics and physics . Springer-Verlag, Berlin,1976. Translated from the French by C. W. John, Grundlehren der Mathematischen Wis-senschaften, 219.[31] A.C. Eringen.
Nonlocal continuum field theories . Springer-Verlag, New York, 2002.[32] A. Ern and J.-L. Guermond.
Theory and practice of finite elements , volume 159 of
AppliedMathematical Sciences . Springer-Verlag, New York, 2004.[33] E.B. Fabes, C.E. Kenig, and R.P. Serapioni. The local regularity of solutions of degenerateelliptic equations.
Comm. Partial Differential Equations , 7(1):77–116, 1982.[34] G. Gilboa and S. Osher. Nonlocal operators with applications to image processing.
MultiscaleModel. Simul. , 7(3):1005–1028, 2008.[35] V. Gol (cid:48) dshtein and A. Ukhlov. Weighted Sobolev spaces and embedding theorems.
Trans.Amer. Math. Soc. , 361(7):3829–3850, 2009. [36] J. Gopalakrishnan and J.E. Pasciak. The convergence of V-cycle multigrid algorithms foraxisymmetric Laplace and Maxwell equations.
Math. Comp. , 75(256):1697–1719 (electronic),2006.[37] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.
J. Comput. Phys. ,73(2):325–348, 1987.[38] M. Griebel, K. Scherer, and M.A. Schweitzer. Robust norm equivalencies for diffusion prob-lems.
Mathematics of Computation , 76(259):1141–1162, February 2007.[39] W. Hackbusch.
Multigrid methods and applications , volume 4 of
Springer Series in Compu-tational Mathematics . Springer-Verlag, Berlin, 1985.[40] W. Hackbusch. The frequency decomposition multi-grid method. I. Application to anisotropicequations.
Numer. Math. , 56(2-3):229–245, 1989.[41] W. Hackbusch. A sparse matrix arithmetic based on H -matrices. I. Introduction to H -matrices. Computing , 62(2):89–108, 1999.[42] H. Harbrecht and R. Schneider. Rapid solution of boundary integral equations by waveletGalerkin schemes. In
Multiscale, nonlinear and adaptive approximation , pages 249–294.Springer, Berlin, 2009.[43] N.S. Landkof.
Foundations of modern potential theory . Springer-Verlag, New York, 1972.Translated from the Russian by A. P. Doohovskoy, Die Grundlehren der mathematischenWissenschaften, Band 180.[44] B.M. McCay and M.N.L. Narasimhan. Theory of nonlocal electromagnetic fluids.
Arch. Mech.(Arch. Mech. Stos.) , 33(3):365–384, 1981.[45] B. Muckenhoupt. Weighted norm inequalities for the Hardy maximal function.
Trans. Amer.Math. Soc. , 165:207–226, 1972.[46] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to fractional diffusion ingeneral domains: a priori error analysis. arXiv:1302.0698, 2013.[47] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Piecewise polynomial interpolation in Muck-enhoupt weighted Sobolev spaces and applications. arXiv:1402.1916, 2014.[48] S.A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces.
J.Mech. Phys. Solids , 48(1):175–209, 2000.[49] E.M. Stein.
Harmonic analysis: real-variable methods, orthogonality, and oscillatory inte-grals , volume 43 of
Princeton Mathematical Series . Princeton University Press, Princeton,NJ, 1993. With the assistance of Timothy S. Murphy, Monographs in Harmonic Analysis, III.[50] R. Stevenson. Robustness of multi-grid applied to anisotropic equations on convex domainsand on domains with re-entrant corners.
Numer. Math. , 66(3):373–398, 1993.[51] Rob Stevenson. Adaptive wavelet methods for solving operator equations: an overview. In
Multiscale, nonlinear and adaptive approximation , pages 543–597. Springer, Berlin, 2009.[52] P.R. Stinga and J.L. Torrea. Extension problem and Harnack’s inequality for some fractionaloperators.
Comm. Partial Differential Equations , 35(11):2092–2122, 2010.[53] B.O. Turesson.
Nonlinear potential theory and weighted Sobolev spaces , volume 1736 of
Lec-ture Notes in Mathematics . Springer-Verlag, Berlin, 2000.[54] Y. Wu, L. Chen, X. Xie, and J. Xu. Convergence analysis of V-cycle multigrid methods foranisotropic elliptic equations.
IMA Journal of Numerical Analysis , 32(4):1329–1347, 2012.[55] J. Xu. Iterative methods by space decomposition and subspace correction.
SIAM Review ,pages 581–613, 1992.[56] J. Xu, L. Chen, and R.H. Nochetto. Optimal multilevel methods for H (grad), H (curl), and H (div) systems on graded and unstructured grids. In Multiscale, nonlinear and adaptiveapproximation , pages 599–659. Springer, Berlin, 2009.[57] J. Xu and L. Zikatanov. The method of alternating projections and the method of subspacecorrections in Hilbert spaces.
Math. Comp. , 15:573–597, 2002.[58] C.-S. Zhang.
Adaptive finite element methods for variational inequalities: theory and appli-cations in finance . PhD thesis, University of Maryland, College Park, 2007.
ULTILEVEL METHODS 29 (L. Chen)
Department of Mathematics, University of California at Irvine, Irvine, CA92697, USA
E-mail address : [email protected] (R.H. Nochetto) Department of Mathematics and Institute for Physical Science andTechnology, University of Maryland, College Park, MD 20742, USA
E-mail address : [email protected] (E. Ot´arola) Department of Mathematics, University of Maryland, College Park, MD20742, USA
E-mail address : [email protected] (A.J. Salgado) Department of Mathematics, University of Tennessee, Knoxville, TN37996, USA
E-mail address ::