Adaptive Multi-level Algorithm for a Class of Nonlinear Problems

: In this article, we propose an adaptive mesh-refining based on the multi-level algorithm and derive a unified a posteriori error estimate for a class of nonlinear problems. We have shown that the multi-level algorithm on adaptive meshes retains quadratic convergence of Newton’s method across different mesh levels, which is numerically validated. Our framework facilitates to use the general theory established for a linear problem associated with given nonlinear equations. In particular, existing a posteriori error estimates for the linear problem can be utilized to find reliable error estimators for the given nonlinear problem. As applications of our theory, we consider the pseudostress-velocity formulation of Navier–Stokes equations and the standard Galerkin formulation of semilinear elliptic equations. Reliable and efficient a posteriori error estimators for both approximations are derived. Finally, several numerical examples are presented to test the performance of the algorithm and validity of the theory developed.


Introduction
Adaptive mesh-refining plays an important practical role in accurate and efficient calculation of the numerical solutions of partial differential equations.A posteriori error estimates are an essential ingredient of adaptivity.They are composed of the known values such as the given data of the problem and the computed numerical solutions.Various a posteriori error estimates for finite element methods for second-order elliptic problems have already been studied in [1,2,8,17,19], for the Stokes problem in [9,26], and for the Navier-Stokes problem in [12,14].Particularly, the nonlinear elliptic equations were considered in [19].There, a priori and a posteriori error estimates measured in the L m (Ω)-norm in the framework of Brezzi-Rappaz-Raviart (BRR) (cf.[4]) were derived.The framework was designed for approximation of branches of nonsingular solutions for a class of nonlinear problems.The theories developed by BRR were applied to study mixed approximation of scalar elliptic problems with gradient nonlinearities in [16,19].The theories applied to the velocity-pressure formulation of the Navier-Stokes equations can be found in [13].The authors of [12,14] constructed a posteriori error estimators for the two-level method for the Navier-Stokes equations.
In this paper, we derive a unified a posteriori error estimate for the approximate solutions of the multilevel algorithm developed in [22] for a class of nonlinear problems in the framework of BRR.The multi-level algorithm using a two-grid idea is an efficient numerical method designed to resolve nonlinearities.The two-grid algorithm (see [20,27,28]) is first applied, and in the subsequent process, uniform mesh refinement is exploited to deliver a desired accuracy and convergence; as the mesh is being refined, the solution on a given mesh is exploited as an accurate starting iterate for solutions on the next mesh level.The multi-level algorithm has quadratic convergence in the sense of equation (2.5) in Theorem 2.5.We apply this strategy to efficiently compute numerical solutions by adaptive mesh refinement.We emphasize that the multi-level algorithm on adaptive meshes retains quadratic convergence of Newton's method across different mesh levels, which is validated from the numerical results presented in the last section.The BRR framework facilitates to use the general theory established for a linear problem associated with given nonlinear equations.In particular, existing a posteriori error estimates for the linear problem can be utilized to find reliable error estimators for the nonlinear problem.We mention that this approach is quite general and can be applied to any numerical scheme for nonlinear problems satisfying the BRR framework.
As applications, we consider the pseudostress-velocity formulation of the stationary Navier-Stokes equations (NSE) and the standard Galerkin formulation of semilinear elliptic equations (SEE).Finite element methods for the velocity-pressure and the stress-velocity-pressure formulation of the NSE have been analyzed for decades (cf.[3,13]).The stress-velocity-pressure formulation is the original physical equations for incompressible Newtonian flows induced by the conservation of momentum and the constitutive law.The mixed finite element methods for the pseudostress-velocity formulation have been recently established in for the Stokes [6,10] and the Navier-Stokes equations [5,21].The pseudostress-velocity formulation has several advantages.The pseudostress is nonsymmetric unlike stress tensor and approximated by the RT k elements, Raviart-Thomas elements of index k ⩾ 0.Moreover, the approximation of the pressure, the velocity gradient, or even the stress can be algebraically obtained from the approximate value of the pseudostress.
For adaptivity, we derive reliable and efficient residual-based a posteriori error estimators for our multilevel algorithm for the NSE.We mention existing works on the adaptive two-level algorithm.The authors in [12,14] derived reliable a posteriori error estimators for the two-level methods for the velocity-pressure formulation of the NSE.In comparison to usual adaptive finite element methods [10,18,19,23,25], the estimators for two-level algorithm contain an additional term which measures the difference between coarse and fine grid solutions.If the additional term measures the difference between consecutive levels and is of higher order (and negligible), then the link between levels does not effect the computation, which is computationally more advantageous.Indeed, this additional term is negligible in our multi-level setting in contrast to [14].See Remark 3.19 for further details.As the second application of the unified framework, we consider the standard Galerkin method for the SEE and derive reliable a posteriori error estimators.It turns out that behavior of the additional term arising from SEE is similar to that of NSE.Note that our algorithm can be applied to other problems such as the elliptic problem with gradient nonlinearities (cf.[16,19,20,24]).
The remainder of this paper is organized as follows.In the next section, we introduce approximation of branches of nonsingular solutions based on the BRR theory and multi-level algorithm analyzed in [22].We derive a unified a posteriori error estimate for the approximate solutions of the multi-level algorithm for a class of nonlinear problems in the framework of BRR.Section 3 and Section 4 are devoted to applying the algorithms to the NSE and the SEE, respectively, and constructing a posteriori error estimators of residual type.Section 5 presents our adaptive multi-level algorithm.In Section 6, numerical results are presented to show the performance of the algorithms and to validate the theory developed in this paper.Concluding remarks are given in the last section.

BRR Theory and Multi-level Algorithm
In this section, we start with the abstract approximation theory of Brezzi, Rappaz and Raviart for nonlinear problems developed in [4,7,13].Then we introduce the multi-level algorithm analyzed in [22], which is a method to solve resulting nonlinear algebraic equations.
From here to the Remark 2.6, we will use the same parts of what we described in [22].Let X, X and Y be real Banach spaces with the norms ‖ ⋅ ‖ X , ‖ ⋅ ‖ X and ‖ ⋅ ‖ Y , respectively, and assume that X → X , continuously imbedding.Let Λ ⊂ ℝ be a compact interval and denote by L (X; Y) the set of all linear and continuous operators from X into Y.For a given C 2 map G : Λ × X → Y, we consider the following nonlinear problem: find where S ∈ L (Y; X ) is a linear operator independent of ν.The set {(ν, ϕ(ν)) : ν ∈ Λ} is called a branch of solutions of (2.1) if F(ν, ϕ(ν)) = 0 and the map ν → ϕ(ν) is continuous from Λ into X.If the Fréchet derivative D ϕ F(ν, ϕ(ν)) of F with respect to ϕ is an isomorphism from X onto X for all ν ∈ Λ, then the branch of solutions {(ν, ϕ(ν)) : ν ∈ Λ} is called nonsingular.To approximate the solution of problem (2.1), we introduce an operator S h ∈ L (Y; X ) intended to approximate the linear operator S, where h > 0 is a real parameter which will tend to zero.The approximation of nonlinear problem (2.1) is to find ϕ h ∈ X such that We assume that the following hypotheses.
(H1) G is a C 2 -mapping from Λ × X into Y, and the second-order Fréchet derivatives of G are bounded on all bounded subsets of Λ × X. (H2) There exists another Banach space Z → Y, continuously imbedding, such that Under the hypotheses, existence of a unique branch of nonsingular solutions of problem (2.2), its a priori error estimate and its convergence property are stated below.We skip the proof since it is very similar to that of [13,Theorem IV.3.3].We introduce the following quantities: where Theorem 2.1.Assume that (H1)-(H4) hold and that {(ν, ϕ(ν)) : ν ∈ Λ} is a branch of nonsingular solutions of (2.1).Then there exist ξ > 0 and h ( ξ ) > 0 such that, for all h ⩽ h ( ξ ), there is a C 2 -function ϕ h : Λ → X satisfying that {(ν, ϕ h (ν)) : ν ∈ Λ} is a branch of nonsingular solutions of (2.2), ϕ h (ν) is the only solution of (2.2) in the ball B(ϕ(ν), ξ ) for all ν ∈ Λ, and for K 1 := 4μ(ν) and a constant C S > 0, the following estimates hold: Note that, since Λ is compact, the function μ(ν) is bounded above on Λ.Thus, if needed, we can take K 1 as the constant independent of not only h but also ν, i.e., K 1 := 4 μ .You may refer to [22] for proofs of Lemma 2.2, Theorem 2.5 and Theorem 2.7 that follow in this section.
Lemma 2.2.Assume that (H1)-(H4) hold and that D ϕ F(ν, φ ) for φ ∈ X is an isomorphism from X onto X .Then there exists ξ ) is an isomorphism from X onto X with the bound, respectively, ) Remark 2.3.When F is replaced by F h for h small enough in Lemma 2.2, estimate (2.4) also holds.That is, if D ϕ F h (ν, φ ) for φ ∈ X is an isomorphism from X onto X , then there exists ξ * = ξ * ( φ ) > 0 such that, for We introduce the multi-level version of the two-grid algorithm.For j ⩾ 0, let h j ⩽ h given in Theorem 2.1 be mesh sizes, X j := X h j finite element spaces such that X j ⊂ X, and F j (ν, ϕ) := ϕ + S j G(ν, ϕ), where an operator S j := S h j is intended to approximate the linear operator S.
Step 2 (Linear solver) Update on each finer mesh level j with one Newton iteration.
For j = 1, 2, . . ., J, find ϕ j (ν) ∈ X j such that D ϕ F j (ν, Algorithm 1: Multi-level algorithm Remark 2.4.The two-grid scheme in [27,28] is based on passing information between finite element equations defined on two grids of different mesh sizes.In the first step, a nonlinear problem itself is solved in a coarse space, i.e., a finite-dimensional space with coarse grids.In the second step, the nonlinear problem is linearized locally at the solution obtained in the coarse space.Then the linearized problem is solved in a fine space.For better performance, this process can be iterated on a sequence of linearized problems with increasing dimensions.
In the second step of the multi-level algorithm in [14], the nonlinear problem is linearized locally at the solution ϕ 0 obtained in the fixed coarse space of Step 1.By contrast, in our multi-level algorithm, the nonlinear problem is linearized locally at the solution obtained on the very previous mesh.The advantages of our algorithm are discussed in detail in Remark 3.19.

Theorem 2.5 (A Priori Error Bound).
Assume that (H1)-(H4) hold.Let ϕ h j (ν) be nonsingular solutions of (2.2) on X j .And let ϕ j (ν) be solutions obtained from Step 2 of Algorithm 1. Then there exist ξ > 0 with ξ ⩽ ξ given in Theorem 2.1 and h * 0 > 0 such that, for h j ⩽ h * 0 , ϕ h j (ν) ∈ B(ϕ(ν), ξ/2) and ϕ j (ν) belongs to the ball B(ϕ h j (ν), ξ/2) for j ⩾ 1.Moreover, for the positive constant K 2 := 4 μ C S L(ϕ, ξ ) independent of mesh sizes h j and ν, we have the following quadratic relation: and for K 3 := max{K 1 , 2K 2 1 K 2 , 2K 2 }, we have an a priori estimate Remark 2.6.It should be noted that the above algorithm can be applied even if T j is not a uniform refinement of T j−1 and adaptive mesh refining can be exploited along with this algorithm.The multi-level algorithm has the quadratic convergence in the sense of equation (2.5) even on adaptive meshes.This is confirmed by numerical examples in Section 6.Now, we want to find an a posteriori error bound for solutions obtained from the multi-level algorithm.First, we introduce the following theorem presented and proved in [7, Theorem 2.1].Let F : X → X be a C 1 mapping and let υ ∈ X be such that DF(υ) ∈ L (X; X ) is an isomorphism.We introduce the notations Theorem 2.7.We assume that 2γL(2γϵ) ⩽ 1.Then F(x) = 0 has a unique solution u in the ball B(υ, 2γϵ) and DF(u) ∈ L (X; X ) is invertible with Moreover, for all υ  ∈ B(υ, 2γϵ), Theorem 2.8 (A Posteriori Error Bound).Assume that (H1)-(H4) hold.Let ϕ(ν) be a nonsingular solution of (2.1) and let ϕ j (ν) be the approximate solution produced by Algorithm 1. Then there exists 0 < h * 1 ⩽ h * 0 given in Theorem 2.5 such that, for all h j ⩽ h * 1 , we have ϕ h j (ν) ∈ B(ϕ(ν), ξ/2) and ϕ j (ν) ∈ B(ϕ h j (ν), ξ/2) for ξ given in Theo-rem 2.5, and we have a residual type a posteriori error bound with K 1 := 4μ(ν), Proof.The proof follows by applying Theorem 2.7.For this, we first notice that ϕ j (ν) ∈ B(ϕ(ν), ξ) for h j with h j ⩽ h * 0 given in Theorem 2.5 and ϕ j (ν) → ϕ(ν) strongly in X as h j → 0 thanks to Theorem 2.5 and (H3).By application of Lemma 2.2, D ϕ F(ν, ϕ j (ν)) is an isomorphism from X onto X with To apply Theorem 2.7, we set From the mean value theorem and (H1), we have, for which yields that there exist 0 < h  0 ⩽ h * 0 and 0 < α < ξ such that, for h j ⩽ h  0 , Applying Theorem 2.5, the mean value theorem and (H1) to the identity From (2.7)-(2.9),we have 2γ j L j (2γ j ϵ j ) ⩽ 2γL j (2γϵ j ) ⩽ 2γL j (α) ⩽ 1.

Application 1: Navier-Stokes Equations
In this section, we introduce the a priori error estimate and analyze an a posteriori error for multi-level algorithms for the pseudostress-velocity formulation of the NSE.The a priori error was analyzed in [21].The following details are only given for the reader's convenience.First, we introduce some notation and function spaces.Let div and ∇ denote the divergence and gradient operators, respectively.For υ = (υ 1 , . . ., υ d ) ∈ ℝ d and τ = (τ ij ) d×d ∈ ℝ d×d , let τ i = (τ i1 , . . ., τ id ) denote its ith row for i = 1, . . ., d and define Given data f and g, the stationary and incompressible NSE for the unknown velocity u and the pressure p reads −νΔu + u ⋅ ∇u + ∇p = f in Ω, (3.1) u = g on ∂Ω. (3. 3) The compatibility conditions read where n is the outward unit normal vector to the boundary.Let A : ℝ d×d → ℝ d×d be the deviatoric operator The pseudostress tensor σ := ν∇u − pδ allows the pseudostress-velocity formulation of NSE (3.1) ) ) Assume that Ω is a bounded, open, connected subset of ℝ d (d = 2 or 3) with the boundary ∂Ω of convex polygon/polyhedra.We use the standard notation and definitions for the Sobolev spaces W r,p (Ω) and W r,p (∂Ω) for an integer r ⩾ 0 and p ∈ [1, ∞].The standard associated norms are denoted by ‖ ⋅ ‖ r,p := ‖ ⋅ ‖ r,p,Ω and ‖ ⋅ ‖ r,p,∂Ω .For r = 0, W r,p (Ω) coincides with L p (Ω).Moreover, the space W r,2 (Ω) will be written in the form H r (Ω).Let We define spaces Ĥ(div; Then the mixed variational problem of the pseudostress-velocity formulation (3.6)-(3.7) is to find where Note that if we use the trial space Ŵ0,3 (div; Ω) d × L 2 (Ω) d , then we can guarantee the nonlinear term, u ⋅ Aσ, to be in L 2 (Ω) d from the imbedding theorem (see [5,Remark 2.1]).The well-posedness and uniqueness of system (3.8)-(3.9) is established by the following well-known theorem (see [5]).[5] or [13] for the specific value of ν 0 ), then the solution (σ, u) is unique.

BRR Framework for NSE
We let Here we used the triangle inequality and the Hölder inequality.So hypothesis (H2) holds.From the definition, G belongs to the C 2 -class clearly.So hypothesis (H1) holds.We let X := L 2 (Ω) d×d × L 2 (Ω) d and introduce the linear solution operator defined by the solution of The following result is well known (see [5,Lemma 5.1]).

Lemma 3.2. For any
Clearly, the pair if and only if the function ϕ = (σ, u) is a solution of (3.8)-(3.9).
We now consider mixed finite element approximations.Let {T h } be a family of shape-regular triangulations of Ω by triangles (d = 2) and tetrahedra (d = 3), respectively, of diameter h T , and let h := max T∈T h h T denote the mesh size.Associated with T h , we define a finite subspace of X by where , a product space of polynomials of degree k) and , a product space of Raviart-Thomas space of order k).For short notation on generic constants C independent of the local or global mesh sizes, for any expressions A and B, A ≲ B abbreviates A ⩽ CB.
We introduce a projection operator over the space W, where W := Ĥ1 (Ω) d×d or W := Ĥ(div; Ω) d ∩ L s (Ω) d×d for some s > 2. Let Πh denote the Raviart-Thomas projection operator (see [3]) associated with the degrees of freedom onto where |Ω| = ∫ Ω dx.Then we have ∫ Ω (tr Π h τ) dx = 0. We notice that the operator Πh approximates the normal components on element edges such that, for E ∈ E h , the set of all edges of T h , for But, in general, the above relation does not hold for the operator Π h .Let P h be the L 2 -projection onto P d k (T h ) with the well-known approximation property Then the following lemma holds (see [6] for details).
Lemma 3.4.The commutative property div Π h = P h div holds.Furthermore, we have The mixed finite element method for approximation of the solution of system (3.8)-(3.9) is defined by finding Now, the discrete nonlinear problem of (3.12) is to find Proof.See [21].
Theorem 3.5 implies that hypotheses (H3) and (H4) hold.We get the following theorem.

Error Estimates for Multi-level Algorithm for NSE
Now, we apply the multi-level algorithm to NSE based on pseudostress-velocity formulation.Note that the Fréchet derivative of for any ψ = (τ, υ) ∈ X.We generate a shape-regular triangulation T j+1 from T j (j ⩾ 1) with newest vertex bisection or longest edge bisection.For convenience, we omit ν from the functions ϕ j (ν).

.24)
Algorithm 2: Multi-level algorithm for NSE If T j with diameter h j is uniformly refined triangulations from the previous triangulations T j−1 with diameter h j−1 , from the results of Theorem 2.5 and Theorem 3.6, the following theorem for Algorithm 1 can be easily proved.
Theorem 2.8 indeed provides the reliability bound.For brevity, ν is dropped from the functions ϕ j (ν) and φ j (ν).First, observe that where is the unique solution for the pseudostress-velocity formulation of the Stokes problem (3.28) Remark 3.8.Equation (3.25) says that an a posteriori error estimate for the NSE is controlled by that for the Stokes equations with the additional term appearing in the right-hand side of (3.28).Thus, we can apply existing a posteriori error estimates for the pseudostress-velocity formulation of the Stokes equations to the NSE.Indeed, we apply the residual-based a posteriori error estimates for the two-dimensional Stokes equations analyzed by Carstensen, Kim and Park in [10].
Let us describe the well-known results used throughout this paper.Assume that Ω is a bounded Lipschitz domain with a polygonal boundary.Then there exists a constant C = C(Ω) such that Poincare's inequalities read ‖υ‖ 0,2 ⩽ C‖∇υ‖ 0,2 for all υ ∈ H 1 0 (Ω) d , where υ Ω is an average of υ over Ω, i.e., υ Ω = 1 |Ω| ∫ Ω υ(y) dy.The general interpolation operator cannot be applied to H 1 functions.But Clément's interpolation which defines a continuous interpolating polynomials using averages of υ instead of point values can be applied to H 1 functions.The nodal values of Clément's interpolating polynomials depend on the value of υ on the adjacent elements through an averaging process.We now consider the error estimate of Clément's interpolation operator.We define S 0 (T j ) d ⊂ L 2 (Ω) d as the piecewise constant and as continuous and piecewise affine vector functions.For each T j , let E j denote the set of all edges (faces) of T j .In what follows, h T and h E stand for the diameters of the triangle T ∈ T j and the edge E ∈ E j , respectively.For any T ∈ T j and E ∈ E j , we set Note that, since the triangulation T j is shape-regular, the maximal number of elements in ω T is h-independently bounded by C M .Clément's interpolation operator (see [11]) I j : H 1 (Ω) d → S 1 (T j ) d satisfies the following estimates: Finally, we state the usual inverse inequalities introduced in [25] which we will use in order to derive local lower bounds for the error.We need some modified results whose proof easily follows from the idea of the paper [25], and we omit the proof.Theorem 3.9.There exist constants C 1 , . . ., C 5 , which only depend on the shape of the elements, such that the following inequalities hold for all polynomials υ, w ∈ P k (T) d * and d ∈ P k (E) d * , d * = d or d × d, and for two real numbers p, q with 1 ⩽ p ⩽ ∞ and 1 p + 1 q = 1: ∫ T υb T wdx ‖w‖ 0,q,T ⩽ ‖υ‖ 0,p,T , (3.32) ) T ‖d‖ 0,q,E , (3.34) ) Let R j (f ) := f + div σ j − 1 ν u j ⋅ Aσ j , ē j = ū − u j and ε j = σ − σ j , and Σ := Ĥ(div; Ω) d and V := L 2 (Ω) d .For simplicity, we let P j := P h j be the L 2 -projection onto P d k (T h j ).Lemma 3.10.For the error ē j (1 ⩽ j ⩽ J), we have where g j is an approximation of g so that ⟨g j , n ⋅ (η − Πj η)⟩ = 0.
We notice that, for ψ ∈ H 1 (ω) d×d and r ∈ H 1 (ω) with the estimate ‖∇z‖ 0,2 where C is a positive constant independent of ε.
Proof.We consider the following Stokes problem: The solutions z and p satisfy −Δz + ∇ p = − div ε and div z = tr ε − s.
We define the jump [ ⋅ ] for E ∈ E j and for υ ∈ H 1 (⋃ T j ) d by and n E points from T + into its neighbor element T − .Lemma 3.12.For the error ε j (1 ⩽ j ⩽ J), we have Proof.
We find, from (3.23) and (3.26), Moreover, from the fact that (A ε , s δ) = (tr A ε , s ) = 0, we get (A ε , Curl(I j r)) = 0.By the identity A σ = ∇ ū and integration by parts (3.40), the second term (II) of (3.44) becomes Considering integration by parts (3.40) and estimate (3.30)-(3.31),recalling that the number of elements in ω T is bounded by the constant C M and using Poincare's inequality (3.29), we obtain )|r| where the error estimator η j is defined by η 2 j := ∑ T∈T j η j (T) 2 with the local error estimator η j (T) given by η j (T) Remark 3.14.We note that physical quantities such as pressure, gradient velocity and stress can be expressed as follows: From these identities, the approximation of pressure, gradient velocity and stress can be defined by Then the following relations hold: Now, we derive an efficient a posteriori error estimate.Let f j ∈ P k (T) d be a polynomial approximation of f .Lemma 3.15.For each T ∈ T j (1 ⩽ j ⩽ J), we have h T ‖R j (f )‖ 0,2,T ≲ h T ‖u − u j ‖ 0,6,T + ‖σ − σ j ‖ 0,3,T + h T ‖f − f j ‖ 0,2,T .
Proof.Using inverse inequality (3.32), the identity div σ = 1 ν u ⋅ Aσ − f and integration by parts, we have By the Hölder inequality, inverse inequality (3.32) and (3.35), and the mean value theorem, we get With an application of the triangle inequality to the identity we complete the proof by (3.47).
Lemma 3.16.For each T ∈ T j (1 ⩽ j ⩽ J), we have Proof.By inequality (3.32) and integration by parts, we have Lemma 3.17.For each T ∈ T j (1 ⩽ j ⩽ J), we have Proof.Inequality (3.32), the identity Curl(Aσ) = Curl(∇u) = 0 and integration by parts imply on each edge E in the weak sense and P ext d| E = d.Using inequality (3.33), integration by parts (3.40) and the identity Curl(Aσ) = 0, we have, for any edge E ⊈ Γ and for ‖Curl(Aσ j )‖ 0,2,T i ). (3.52) From Lemma 3.17, we have Remark 3.19.By the Hölder inequality and the triangle inequality, we have Assume that h 2 j−1 ⩽ h j for j > 0.Then, by the result of Theorem 3.7, Similarly, , which shows the asymptotic behavior of the additional term.Note that the additional term is of higher order compared to the other error estimates and the error ‖u j − u‖ 0,6,Ω + ‖σ j − σ‖ 0,3,Ω = O(h r j ) if we take h j = O(h 2−ϵ j−1 ) for some positive ϵ.In practice, adaptively refining mesh via newest vertex bisection or red refinement satisfies the relation h 2 j−1 ≪ 1 4 h j−1 ⩽ h j , which implies that the additional term can be negligible in the adaptive algorithm in contrast to the result in [14], where this additional term is not negligible as refinements proceed.The reason is that, in the second step of his two-grid algorithm, the nonlinear problem is linearized locally at the solution obtained in the fixed coarse space.When the linearized problem in the second step is solved in a fine mesh, as the mesh refinement proceeds, the additional term becomes relatively large.This issue can be resolved by using our multi-level algorithm, where the nonlinear problem is linearized locally at the solution obtained on the very previous mesh.
Combining , we arrive at the following efficiency estimates.

Theorem 3.20 (Efficient A Posteriori Error Estimate).
For the solution ϕ = (σ, u) of problem (3.8)-(3.9)and the approximate solution ϕ j = (σ j , u j ) to problem (3.23)-(3.24),we have Note that the derived a posteriori error estimator is efficient and reliable, i.e., the efficiency index EI := η j /‖ϕ − ϕ j ‖ X and its reciprocal are bounded for all triangulations.

Application 2: Semilinear Elliptic Equations
We consider the semilinear boundary value problems of the form where is a given mapping.This problem has been treated in [7] for d = 2. Multiplying by test functions and integrating by parts, we obtain the following weak formulation of problem (4.1):

Error Estimates for Multi-level Algorithm for SEE
We generate a shape-regular triangulation T j+1 from T j (j ⩾ 1) with newest vertex bisection or longest edge bisection.For simplicity, we let X j := X h j .Then, for problem (4.1) with G(λ, u) = u α − f , the multi-level algorithm reads as follows.
Step 1 (Nonlinear solver) Find u 0 ∈ X 0 such that Step 2 (Linear solver) For 1 ⩽ j ⩽ J, find u j ∈ X j such that Theorem 4.1 (A Priori Error Estimates).Let f ∈ Y. Let u be a solution of (4.1), u 0 the solution of (4.6), and u j the solution of (4.7).Assume that u ∈ H r (Ω) for 2 ⩽ r ⩽ k + 1.Then we have, for a constant C > 0 that depends on ‖u‖ 2,2 and independent of mesh sizes, [20,22].Remark 4.2.By referring to [20], we can prove that the result of (2.5) holds for X = L 2 (Ω) as well as X = H 1 0 (Ω).Moreover, we see that if the exact solution u is in W 2 p (Ω) with p > 2 for d = 2 and p = 12 5 for d = 3, then the following result holds: which shows improvement O(h 4 j−1 ) compared with the result O(h 2 j−1 ) via (2.5).
Proof.By integration by parts and (4.10), we get where P j is L 2 -projection onto X j .Using the error estimates of L 2 -projection, we conclude the proof.In fact, to estimate the term (w j , P j ē ), we use Hölder inequality and Sobolev imbedding theorem.Then we have, for α = 2, and for α = 3, Remark 4.4.Assume that h 2 j−1 ⩽ h j for j > 0.Then, by the triangle inequality, the Sobolev imbedding theorem, which is that H 1 (Ω) is continuously imbedded in L p (Ω) for p ∈ [1,6], and Theorem 4.1, we have, for q = 3 or 4, where a constant C depends on ‖u‖ 2,2 and is independent of mesh sizes.As in the case of the Navier-Stokes equations, the additional term w j is of higher order O(h j−1 ) compared to the other error estimates and the error ‖u j − u‖ 1,2 = O(h r−1 j ) if we take h j = O(h 2−ϵ j−1 ) for some positive ϵ.Efficient a posteriori error estimates can be derived in a way similar to those addressed in the NSE using the usual inverse inequalities introduced in [25].

Adaptive Multi-level Algorithm
We introduce an adaptive multi-level algorithm for nonlinear problems based on the local error estimator η j (T) and the global error estimator η j driven in Section 3 and 4. The adaptive procedure for the nonlinear problem goes as follows: find ϕ 0 from Step 1 of Algorithm 1 on a given initial (uniform) mesh T 0 , form a mesh T 1 uniformly refined from T 0 , and generate an adaptively refined mesh T j+1 from T j (j ⩾ 1) by repeating the following process: The step SOLVE means to perform Step 2 of Algorithm 1 on the mesh T j (j ⩾ 1).This process is iterated until the global error estimator η j becomes smaller than a given tolerance.We organize the adaptive multi-level algorithm for the nonlinear problems.
Input: Initial (uniform) mesh T 0 , parameter θ, and tolerance tol.Output: Final mesh T J and approximate solution ϕ J for J ⩾ 1. Find ϕ 0 from Step 1 of Algorithm 1 on T 0 .Form a mesh T 1 uniformly refined from T 0 .Obtain ϕ 1 from Step 2 of Algorithm 1 on T 1 .Compute the local error estimator η 1 (T) and global error estimator η 1 .j := 1. while η j ⩾ tol do Mark a set M T j ⊂ T j such that ∑ T∈M T j η j (T) 2 ⩾ θη 2 j (bulk strategy).j := j + 1. Generate a conforming triangulations T j from T j−1 by refining T ∈ M T j−1 .Obtain ϕ j from Step 2 of Algorithm 1 on T j .Compute the error estimators η j (T) and η j .

Algorithm 4: Adaptive multi-level algorithm
For the reader's convenience, we present Figure 1 in order to explain the outline of the (uniform/adaptive) multi-level algorithm.The details of the uniform multi-level algorithm can be found in [22].Figure 1 depicts the multi-level algorithm compared with the usual nonlinear solver (e.g., Newton-Raphson algorithm).The nonlinear solver finds discrete nonlinear solutions by iteration on a much higher-dimensional fixed space X J than X 0 .The multi-level algorithm finds the solution ϕ 0 from the nonlinear solver on the low-dimensional space X 0 and then repeatedly finds ϕ j by only one iteration on X j with starting value ϕ j−1 for j = 1, . . ., J. The notations

Multi-level algorithm
Nonlinear solver

Figure 1:
Multi-level process computing ϕ j ∈ X j ⊂ X, j = 0, 1, . . ., J, with J finest level k=0 denote the iteration sequence obtained from the nonlinear solver on the space X 0 and X J , respectively.The thickness of the arrows indicates the computational costs at the given level; thicker means higher costs.

Numerical Experiments
In this section, we perform various numerical experiments to test our adaptive multi-level algorithm and to validate the theory analyzed in this paper.We present the convergence of various errors and indicators on both uniform and adaptive meshes.In particular, we will numerically verify that the adaptive multi-level algorithm has the quadratic convergence in the sense of (2.5).

Semilinear Elliptic Equations
We approximate the semilinear elliptic problem (4.1) with G(λ, u) = u α − f , α = 2, 3, by using a conforming finite element space of piecewise linear polynomials.We introduce the following quantities: for 1 ⩽ j ⩽ J, We define the local error estimator η j (T) by η j (T) := (E j,1 (T) 2 + E j,2 (T) 2 + E j,3 (T) 2 ) 1/2 , the global error estimator η j by , and efficiency indices EI j by To show the performance, we consider two test examples: one with a boundary layer, the other with a corner singularity.In Table 1, we present convergence history of multi-level algorithm and nonlinear solver on the uniform meshes.
Here, N el means the number of elements at T j .The solution obtained from the multi-level algorithm is as accurate as the nonlinear solver at each level j.But we observe that various errors up to low levels (j ⩽ 2) exhibit non-optimal convergence behavior since meshes in low levels are in the pre-asymptotic range.We see that the errors have optimal order of convergence by carrying out more uniform refinement (j ⩾ 3).However, in this uniform refinement, for high levels j ⩾ 3, the number of elements at T j increases rapidly and the finite element space X j requires too many degrees of freedom, which necessitates the adaptive multi-level algorithm.In Table 2 and Figures 2 and 3, we present numerical results from the adaptive multi-level algorithm.Efficiency indices EI j = η j /|u − u j | 1,2 are presented in Table 2.The history of errors |u − u j | 1,2 and η j , and the efficiency indices EI j are plotted with a logarithmic scale in Figure 2, where the x-axis represents the number of elements.Numerical data for EI j show that the a posteriori error estimate is numerically reliable and efficient.The superiority of the adaptive multi-level algorithm can be clearly seen in Tables 1 and 2. Indeed, e.g., when α = 2, the accuracy on uniform mesh with N el = 262,144 in Table 1 is comparable to the accuracy on adaptive mesh with N el = 36,448 in Table 2 (see boldface in Tables 1 and 2). Figure 3 shows convergence orders of the indicators.The symbols in the legend of Figure 3 represent the indicators such as In Figure 3, we observe that the additional term E j,3 is of higher order than the other indicators for both uniform and adaptive meshes as mentioned in Remark 4.4.
To examine quadratic convergence behavior in the sense of (2.5) for X = L 2 (Ω) and in the sense of (4.8), we investigate the difference between u j−1 and u h j , and the difference between u j and u h j in the L p or H 1 -norm, and their ratios.In Figure 4, we present these quantities and the ratios RL j and RH j , where The ratios RL j and RH j are almost constant, which means that our adaptive multi-level algorithm retains quadratic convergence of Newton's method across different mesh levels.We present the adaptive meshes T j for j = 0, 1, 2, 10 and the approximate solution u 10 on the final mesh T 10 for α = 2 in Figure 5.We see that the adaptively designed meshes properly express the characteristics of the solution with a boundary layer in the upper right corner.We note that Δυ = 0 for υ := r 2/3 sin(2θ/3) ∈ H s (Ω) with s < 5/3, and for β > 1, w := r β sin(2θ/3) ∈ H 2 (Ω) and G(u) = u α − f = Δw ∈ L 2 (Ω).We take β = 1.01 for the numerical test.Note also that a priori error estimates for the Poisson model case on the L-shaped domain are well known (see [15,Remark 4.1, Example 4.5 and Problems 4.7]); for any ϵ > 0, To take the Neumann part of the boundary condition into account, we modify indicator E j,2 as follows: We only provide results for α = 2 as the case α = 3 shows similar behavior.We observe from Table 3 that convergence orders on the uniform meshes are in par with the Poisson problem (cf.(6.1)).In Table 4 and Figure 6, we display history of errors on the adaptive meshes for each level j, which shows the full rate of convergence.The superiority of the adaptive multi-level algorithm can be clearly seen in Tables 3 and 4. Indeed, the accuracy on uniform mesh with N el = 24,576 in Table 3 is comparable to the accuracy on adaptive mesh with N el = 3,705 in Table 4 (see boldface in Tables 3 and 4).In Figure 7 (a) and (b), we present indicators and their convergence orders on uniform and adaptive meshes, respectively.The convergence order of additional term E j,3 is higher than that of the other indicators E j,1 and E j,2 .The ratios RL j and RH j are almost constant in Figure 7 (c) and (d), which means that our adaptive multi-level algorithm retains quadratic convergence of Newton's method across different mesh levels.
We present adaptive meshes T j for j = 0, 1, 2, 6 and the approximate solution u 6 on the mesh T 6 in Figure 8.We see that the adaptively designed meshes well reflect the corner singularity of the solution.
We define the local error estimator η j (T) by Table 5 shows convergence history of multi-level algorithm and nonlinear solver on uniform meshes.We observe that the errors ‖ϕ − ϕ j ‖ X up to low levels (j ⩽ 2) exhibit non-optimal convergence behavior since meshes in low levels are in the pre-asymptotic range.We see that the errors have optimal order of convergence by carrying out more uniform refinement (j ⩾ 3).However, in this uniform refinement, for high levels j ⩾ 3, the number of elements at T j increases rapidly and the finite element space X j requires too many degrees of freedom as in Example 6.1; the adaptive multi-level algorithm is beneficial.

Conclusion
The adaptive multi-level algorithm analyzed in this paper can be used to efficiently solve nonlinear problems with singularity or layer.We start with a nonlinear system on a very coarse mesh and then solve linearized system on the successively refined meshes with adaptivity.The BRR framework facilitates to use the general theory established for a linear problem associated with given nonlinear equations.In particular, a posteriori error estimates for the linear problem can be utilized to find reliable error estimators for the nonlinear problem.
As applications, the Navier-Stokes equations and the semilinear elliptic equations are considered.Residual type reliable and efficient a posteriori error estimators are constructed for the adaptive multi-level algorithm for the considered problems.With few modifications we apply existing a posteriori error estimators for the linear Stokes equations and elliptic equations corresponding to the NSE and SEE, respectively.Note that the multi-level algorithm with uniform meshes has quadratic convergence in the sense of (2.5).We emphasize that the multilevel algorithm with adaptive meshes retains quadratic convergence of Newton's method across different mesh levels, which is validated from the numerical results presented in this paper.

. 27 )
The following error equations immediately follow from (3.23)-(3.24)and (3.26)-(3.27): .36) Here, for the barycenters x T of T and x E of E, two functions b T , b E ∈ C ∞ (T, ℝ) satisfy b T ∈ P 3 (T), b T (x T ) = 1 and b T = 0 on ∂T, and b E ∈ P 2 (T), b E (x E ) = 1 and b E = 0 on ∂T\E, respectively, and P ext is an extension operator P ext : C(E, ℝ d * ) → C(ω E , ℝ d * ) with P ext w| E = w.

3 Figure 2 :Figure 3 :Figure 4 :
Figure 2: The history of errors |u − u j | 1,2 and η j , and EI j on adaptive meshes for Example 6.1

Figure 5 :
Figure 5: Meshes and u 10 generated by the adaptive multi-level algorithm for Example 6.1

Figure 7 :
Figure 7: The convergence of indicators EI j and ratios RL j and RH j for Example 6.2

Table 1 :
Convergence orders (CO) of various errors on uniform meshes for Example 6.1

Table 2 :
Efficiency indices on adaptive meshes for Example 6.1

Table 3 :
Convergence orders (CO) of various errors on uniform meshes for Example 6.2