{'fraction_non_alphanumeric': 0.06414732399769807, 'fraction_numerical': 0.02451563399194322, 'mean_word_length': 3.638903719523047, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasi-Optimality of an Adaptive Finite Element Method for Cathodic Protection

In this work, we derive a reliable and efficient residual-typed error estimator for the finite element approximation of a 2d cathodic protection problem governed by a steady-state diffusion equation with a nonlinear boundary condition. We propose a standard adaptive finite element method involving the D\"{o}rfler marking and a minimal refinement without the interior node property. Furthermore, we establish the contraction property of this adaptive algorithm in terms of the sum of the energy error and the scaled estimator. This essentially allows for a quasi-optimal convergence rate in terms of the number of elements over the underlying triangulation. Numerical experiments are provided to confirm this quasi-optimality.


Introduction
Let Ω be a bounded polygonal domain in R 2 with its boundary Γ consisting of three mutually disjoint parts: Γ := Γ 0 ∪ Γ A ∪ Γ C , all of which are line segments. This work is concerned with the numerical treatment of the following problem: − ∇ · (σ∇u) = 0 in Ω, (1.1) σ ∂u ∂n = 0 on Γ 0 , σ ∂u ∂n = g on Γ A , σ ∂u ∂n where the conductivity σ is assumed to be a piecewise W 1,∞ function such that σ 1 ≤ σ ≤ σ 2 a.e. in Ω with two positive constants σ 1 and σ 2 , n is the unit outward normal on Γ and g ∈ L 2 (Γ A ). The system (1.1)-(1.2) arises in cathodic protection in electrochemistry. In a container Ω occupied by electrolyte, the first boundary condition in (1.2) describes insulation of the surface Γ 0 by painting. The second boundary equation in (1.2) reflects the fact that a current density g on anodes Γ A induces an electrical potential u in Ω governed by (1.1). The corrosion process on cathodes Γ C is slowed down through the nonlinear relation f , which depends on the electrode material and is given by [16] either f 1 (u) = C 1 u + C 2 u 3 (1. 3) or the Butler-Volmer function f 2 (u) = C 5 (e C 3 u − e −C 4 u ), (1.4) where C i , i = 1 · · · 5, are all positive constants. In problem (1.1)-(1.2), the sudden change of the boundary condition from Neumann type on anodes and the insulated part to a nonlinear one on cathodes gives rise to local solution singularities in these regions. Furthermore, internal layers may appear due to the discontinuity of the conductivity. Consequently, the computational efficiency will be compromised if a uniform mesh refinement is employed in the finite element discretization. One remedy in practice is to employ adaptivity techniques featuring local refinement so that numerical results can attain better accuracy with minimum degrees of freedom. The aim of this work is to investigate the computational complexity of an adaptive finite element method (AFEM) for problem (1.1)-(1.2).
A typical adaptive algorithm comprises successive iterations of the following loop: That is, SOLVE yields a finite element approximation on the current mesh; ESTIMATE computes the relevant a posteriori error estimator; MARK picks some elements to be subdivided; REFINE produces a new finer mesh. The module ESTIMATE, depending on some computable quantities, i.e., the discrete solution, local mesh size and given problem data, plays an indispensable role in (1.5). Since the seminal work [2], a posteriori error estimation for FEMs has been well understood in scientific computing and engineering [1,26]. As to the mathematical theory of AFEM, e.g., convergence and computational complexity, there have been great developments (see the overview [6,21] and the references therein) over the past thirty years. For linear elliptic problems, this issue has been investigated at depth [4,7,10,12,23]. Recently, the analysis has been extended to some nonlinear problems; see [3,9] for p-Laplacian and [14,15] for quasi-linear equations.
Recently, we [19] have proposed an AFEM of the form (1.5) for problem (1.1)-(1.2) and proved its plain convergence, namely the H 1 -norm error and the sequence of relevant estimators both go to zero as the loop (1.5) proceeds. This work is a continuation of [19], and it is devoted to the complexity of the algorithm. In the AFEM, ESTIMATE, MARK and REFINE use a residual-type a posteriori error estimator, Dörfler strategy and the bisection [18,24], respectively. The main contributions include a contraction property in Theorem 5.1 and the quasi-optimality computational complexity in terms of the number of elements associated with underlying triangulations in Theorem 6.1.
Our analysis is inspired by [12], to first obtain optimal marking for the error estimator, cf. Lemma 6.2 and then the optimal decay rate for the energy error plus an oscillation term so that the upper bound of the parameter in Dörfler strategy is independent of the efficiency constant. However, the nonlinear term f on the boundary Γ C requires a different treatment. First, due to the presence of the nonlinear term f on the boundary Γ C , the Galerkin orthogonality fails. We employ the energy functional instead of the energy norm as in [9,15]. By the equivalence of J (u T ) − J (u) and u − u T H 1 (Ω) , cf. Lemma 4.2, we prove that the adaptive algorithm reduces the sum of energy error and the scaled estimator for any two consecutive iterations. Second, instead of standard arguments for linear problems, we use the generalized Hölder inequality and the stability of solutions to establish a Céa-type lemma for complexity estimate, cf. Lemma 4.5.
The remainder of the paper is organized as follows. Section 2 is devoted to the a posteriori error analysis. An adaptive algorithm to approximate problem (1.1)-(1.2) is described in Section 3. We prove the convergence of this algorithm by a contraction property in Section 5 after presenting preliminary results in Section 4. Section 6 focuses on the quasi-optimal convergence rate. Throughout, we adopt standard notation for Sobolev spaces and related norms and seminorms. Moreover, any generic constant, with or without subscript, is independent of the mesh size and is not necessarily the same at each occurrence.

A posteriori error analysis
In this section, we shall derive a residual-type error estimator for the finite element approximation of problem (2.1), which forms the basis of our AFEM. To introduce the AFEM, we first recall the variational formulation of problem (1.1)-(1.2): find u ∈ H 1 (Ω) such that We refer to [16,17] Then by the application of the Poincaré inequality, the boundedness of σ and the trace theorem, we arrive at where β 1 and β 2 are positive constants depending on σ, α, Γ C and Ω. Utilizing (2.1) with v := u, together with the mean value theorem, (2.2) and (2.3), we can obtain Then the trace theorem yields the a priori estimate to problem (2.1): with C stb > 0 being a constant depending on σ, α, Γ C and Ω. Furthermore, the continuity of the imbedding H 1 (Ω) ֒→ H 1 2 (Γ) ֒→ L q (Γ) in 2d for all q < ∞ implies the existence of a constant C imb,q > 0 depending on Ω and q, satisfying (2.5) Next, we proceed to the discretization. Let T be a shape-regular conforming triangulation ofΩ into a set of disjoint closed triangles such that the coefficient σ is piecewise W 1,∞ over T 0 . For each element T ∈ T , we denote its mesh size h T := |T | 1 2 and ρ T the diameter of the largest inscribed ball. We associate each triangulation T with its shape regular parameter Over the mesh T , we consider the usual H 1 -conforming finite element space V m T , consisting of all piecewise polynomials of degree less than or equal to m ∈ N + , i.e., Then the discrete problem corresponding to (2.1) reads: Similar to the continuous case, the following stability estimate holds To describe the error estimator, we need a few notation and definitions. The collection of all edges (resp. all interior edges) in T is denoted by F T (resp. F T (Ω)) and its restriction on Γ (resp. Γ 0 , Γ A and Γ C ) by F T (Γ) (resp. F T (Γ 0 ), F T (Γ A ) and F T (Γ C )). The scalar h F := |F | stands for the diameter of F ∈ F T , which is associated with a fixed normal unit vector n F inΩ with n F = n on the boundary ∂Ω. For each T ∈ T , we denote ω T as the union of all elements in T with non-empty intersection with element T . For any F ∈ F T , ω F is the union of two elements that share F . Further, we let T be the Scott-Zhang quasi-interpolation operator over T [22]. Then for all v ∈ H 1 (Ω), T ∈ T and F ∈ ∂T ∩ F T , there holds with C I a constant depending on the shape regularity parameter C T . For any v T ∈ V m T , we define the residuals on each element T ∈ T and each edge F ∈ F T by where [·] denotes jumps across interior edges F : Then the local error indicator on any element T ∈ T is defined by The error estimator over the element patch M ⊆ T is Similarly, the oscillation term can be defined locally and globally by The following upper and lower bounds on the error estimator were given in [19]. Here we give a more precise estimate with respect to the occurring constant. For completeness, we provide the proof, since a related argument will be used in the proof of Lemma 4.7.
Theorem 2.1 (Reliability). Let u ∈ H 1 (Ω) and u T ∈ V m T be the solutions to problems (2.1) and (2.6), respectively. Then there exists a positive constant C rel > 0 depending on σ, α, Ω, Γ C , m and C T such that The discrete variational equation (2.6) implies These two estimates together yield By elementwise integration by parts and taking v T := I sz T v, we obtain Then a combination of (2.8) and Young's inequality completes the proof with C rel := 4C 2 I C ov /β 2 1 .

Theorem 2.2 (Efficiency).
Let u ∈ H 1 (Ω) and u T ∈ V m T be the solutions to problems (2.1) and (2.6), respectively. Then there exists a positive constant C ef f depending on σ, α, Ω, Γ C , m and

Adaptive algorithm
Now we present the AFEM for problem (2.1). Let T be the set of all possible conforming triangulations ofΩ obtained from some initial mesh by successive bisections [18,20,24]. The refinement process ensures that all constant depending on the shape regularity of T ∈ T are uniformly bounded by a constant only depending on the initial mesh [21,25]. T * is a called refinement of T for T ∈ T, if T * ∈ T is produced from T by a finite number of bisections.
The proposed adaptive algorithm is given below. For each triangulation Algorithm 3.1. Given an initial conforming mesh T 0 and a parameter θ ∈ (0, 1]. Set k := 0.

(REFINE)
Refine each T ∈ M k by bisection to get T k+1 .
The convergence and quasi optimality of Algorithm 3.1 will be analyzed in Sections 5 and 6. A key ingredient is the so-called closure estimate over the meshes {T k } k : with the constant C 0 depending on C T 0 and #T denoting the number of elements in T . This estimate was first proved in [

Auxiliary results
This section is devoted to several technical lemmas for the convergence analysis of Algorithm 3.1. As is well known, Galerkin orthogonality or Pythagoras property is key to the convergence analysis of the linear problems, which regretfully fails for the nonlinear case. Thus a new equivalent error has to be developed that can play the role of the Galerkin orthogonality property.
First, we introduce the associated functional to (2.1) by Then problem (2.1) is equivalent to the minimization problem [19]: Let u ∈ H 1 (Ω) and u T ∈ V T be solutions to problems (2.1) and (2.6), respectively. The nonnegative quantity is referred to as the equivalent error throughout this paper. The next lemma [17, Lemma 2.1] is useful to handle exponential nonlinearity in (1.4).
Next we show the equivalence between E(u T ) and u − u T 2 H 1 (Ω) .
Lemma 4.2. Let u ∈ H 1 (Ω) and u T ∈ V T be solutions to problems (2.1) and (2.6), respectively. Then there holds with positive constants c equ and C equ depending on σ, α, Ω, Γ C , m and C T .
Proof. Let y(t) := J (w(t)) for t ∈ [0, 1], with w(t) := (1 − t)u + tu T . Since w(0) = u is the minimizer of J over H 1 (Ω), consequently, we can obtain y ′ (0) = 0. In the meanwhile, we can obtain by Taylor's theorem that To finish the proof, we need to compute y ′′ (t). In view that F ′ (·) = f (·), an application of the chain rule implies which, together with the identity and the definition (4.1), yields Combining with (4.2), we arrive at After combining with (2.2) and (2.3), this derives the lower bound with c equ := 1 2 β 1 . In the following we will prove the upper bound. The convexity of f ′ implies We will discuss the cases f = f 1 and f = f 2 separately.
In the case of (1.4), f := f 2 . Then we can derive by application of the generalized Hölder inequality, together with Lemma 4.1, (2.4) and (2.5), that Finally, we can complete the proof after collecting these estimates above.
Remark 4.1. Note that let T * be a refinement of T and u T * be the solution to problem (2.6) over T * . Then the estimate in Lemma 4.2 still holds should u be replaced with u T * . Note also that when proving the contraction property of Algorithm 3.1 in Section 5, we will resort to the following identity instead of Pythagoras property in the form of energy norm, which is an important ingredient in relevant arguments for linear elliptic problems, but now fails for the nonlinear problem (2.1). The remaining of this section is concerned with the auxiliary results in the analysis: estimator reduction, Céa's lemma, oscillation perturbation and discrete reliability. Let T ∈ T, M ⊂ T and T * ∈ T be obtained from T by Algorithm 3.1 with M being the marked set. Let u T ∈ V T and u T * ∈ V T * be the solutions to problem (2.6) over T and T * respectively. Then there exists a constant C est depending only on σ, α, Ω, Γ C , C T , m and g L 2 (Γ C ) satisfying Before proceeding to its proof, we need an auxiliary result: Lemma 4.4 (Local perturbation of estimator). Let T ∈ T. Then there holds ∀δ ∈ (0, 1) : for all v T , w T ∈ V m T and T ∈ T . Here, the constant C sym,1 depends on C T , m, Ω, Γ C , σ, v T H 1 (Ω) and w T H 1 (Ω) .
Proof. Let e := v T − w T . Firstly, we can obtain by the definition (2.9) combining with an application of the Young's inequality that The inverse estimate indicates for some positive constant C inv depending only on C T and m, there holds Combining those two estimate, we can obtain ∀δ ∈ (0, 1) : Here, the constant C aux,1 depends on C T , m and σ. Therefore, it suffices to bound the last term.
For f := f 1 in (1.3), a direct calculation leads to The inverse estimate [5,Section 4.5] gives with the constant C inv,q depending on C T and m. The combination of these two inequalities with q := 6, together with the generalized Hölder's inequality and (2.5), yields ). Together with (4.5), this proves the desired result with C sym,1 := max{C aux,1 , h T C aux,2 }.
For f := f 2 in (1.4), note that the convexity of f ′ (·), together with the mean value theorem, implies for some s ∈ (0, 1), there holds Consequently, taking square on both sides and employing the definition (1.4), leads to Then integrating both sides over F and invoking the generalized Hölder's inequality, in combination with (4.6) with q := 4 and Lemma 4.1, we get Here, Therefore, noticing (4.5), the assertion follows with C sym,1 := max{C aux,1 , h Proof of Lemma 4.3. Let T ∈ T * . Utilizing Lemma 4.4 with T := T * and v T := u T , w T := u T * , we can derive

(4.9)
Summing over all elements T ∈ T * leads to

Now the Sobolev imbedding theorem (2.5) implies
with C est := C sym,1 (1 + 1 δ )(C ov + C 2 imb,2 ). Note that thanks to the stability estimate (2.7) for u T and u T * , the positive constant C est depends only on σ, α, Ω, Γ C , C T , m and g L 2 (Γ C ) . Note also that In view that for each element T ∈ T * \T , there exists a uniqueT ∈ T \T * , s.t., T ⊂T and h T ≤ 1 √ 2 hT . Then by definition (2.9), we arrive at Finally, the combination of (4.10), (4.11), (4.12) and M ⊆ T \T * yields the assertion.
Remark 4.2. If we exchange u T and u T * in (4.9) and then sum up it over T ∩ T * instead, a similar argument leads to Proof. We derive from (2.12) and the Galerkin Orthogonality that for any v T ∈ V T , there holds We focus on the second term, and claim For therefore, an application of the generalized Hölder's inequality implies Thus the claim follows from (2.5), (2.4) and (2.7) with C aux,4 : Multiplying by u − v T , integrating over Γ C and applying generalized Hölder's inequality and Lemma 4.1 yield Thereafter, we establish the claim by appealing to (2.5), (2.4) and (2.7) with C aux,4 : = C 5 C 2 imb,4 C 3 2 + 2|Γ C | 1/2 + 2|Γ C | 1/2 e 2C 2 3 CexpC 2 stb g 2 .
To establish the oscillation perturbation estimate, we may follow the proof of Lemma 4.4 to obtain the local perturbation of oscillation. However, this leads to the issue that the related constant depends on v T and w T as C sym,1 in Lemma 4.4. As a result, the constant in the subsequent convergence rate involves the finite element function v T , which should be avoided in order to establish the quasi-optimality estimate. To this end, we keep the nonlinear function in the following estimate.
Lemma 4.6 (Oscillation perturbation). Let T , T * ∈ T with T * being a refinement of T . If v T ∈ V T and v T * ∈ V T * , then where C op depends on C T , m, σ and C ov , and (4.14) Proof. Let T ∈ T ∩ T * and denote e := v T − v T * . We can obtain from the definition (2.10) and the Young's inequality that Together with the inverse estimate (4.4), we arrive at where C aux,5 depends on C T , m and σ.
Let f = f 1 . Note that since J F (v) ∈ P 3m (∂T ∩ Γ C ) for all v ∈ V T and T ∈ T , and sincē J F (v) is the L 2 -projection of J F (v) on P 3m (F ), therefore, we can obtain ∀v ∈ V T and F ∈ ∂T ∩ Γ C : Plugging those estimates above into (4.15) leads to In the meanwhile, note that osc T * (v T , T ) = osc T (v T , T ) and v T ∈ V T * , then summing over T ∈ T ∩ T * in (4.16) leads to The desired result follows by taking δ := 1 and C op := 2 × max{C aux,5 C ov , 2}.
Lemma 4.7 (Discrete reliability). Let T , T * ∈ T with T * being a refinement of T and let u T ∈ V T , u T * ∈ V T * be solutions to problem (2.6) over T and T * , respectively. Then there exists C drel > 0 depending only on σ α, Ω, Γ C , m and C T such that Proof. Using the operator I sz T (2.8) to v := u T * − u T and noting v = I sz T v on unrefined elements in T ∩ T * , the argument of Theorem 2.1 completes the proof with C drel := C rel .

Convergence
Now we show each iteration of Algorithm 3.1 reduces the sum of the equivalent error and the scaled estimator, which implies the convergence of the algorithm.

Quasi-optimality
Now we give a quasi-optimal convergence rate for Algorithm 3.1. We begin with a generalization of Cea's lemma in Lemma 4.5.
Lemma 6.1. Let u and u T be solutions to problems (2.1) and (2.6) over some mesh T ∈ T.
Here, the positive constant C qs depends on C T , m, σ, g and C ov .
Proof. Given v T ∈ V T . An application of Lemma 4.6 with v T := u T , T * := T yields Then an application of the triangle inequality yields Combining with Lemma 4.5, we arrive at Let f := f 1 . Thenf = 0 by (4.14). This proves the assertion by taking C qs := max{2, 2C op + 2C op C cea + C cea }.
Remark 6.1. We refer to the square root of the left hand side of (6.1) as the total error. Therefore, Lemma 6.1 establishes the quasi-optimality of the solution u T in the sense of the total error. Note that the right hand side of (6.1) can also be controlled by the total error: This estimate can be derived from (6.3) directly.
Next, we introduce the approximation class. Let T N ⊂ T be a subset consisting of all triangulation T ∈ T satisfying #T − #T 0 ≤ N . The approximation class A s for 0 < s ≤ m/2 is defined by The upper bound m/2 is attained for the uniform refinement. Then we give the fundamental ingredients in the analysis, i.e., the optimal marking and cardinality of M k . We follow [12] to derive the optimal marking that relates a strict error estimator reduction to Dörfler marking. This type of estimate was first given in [23] for the Poisson equation with an H 1 -norm reduction, and then extended in [7] to the total error for linear elliptic problems. Below, we present a version in terms of the error estimator as in [12], the proof of which does not require the efficiency estimate in Theorem 2.2. Lemma 6.2 (Optimal marking). Suppose that the marking parameter θ in (3.1) satisfies θ ∈ (0, 1/(1 + C est C drel )). (6.4) Let T * ∈ T be any refinement of T ∈ T and let u T ∈ V T , u T * ∈ V T * be solutions to problem (2.6) over T and T * , respectively. For any δ > 0, define with a fixed k ∈ N + , there exist a triangulation mesh T ε ∈ T and v ε ∈ V Tε such that (6.7) Let T * be the smallest common refinement of T ε and T k , i.e., T * := T ε ⊕ T k . Then T * is a refinement of T ε . An application of Theorem 2.2 and (6.1), combining with the inequality osc 2 T * (v ε ) ≤ osc 2 Tε (v ε ) and (6.7), leads to This, together with Theorem 2.1, the inequality osc 2 k (u k ) ≤ η 2 k (u k ) and (6.6), gives η 2 T * (u T * ) ≤ λη 2 k (u k ). Consequently, the subset T k \ T * satisfies the Dörfler marking strategy owing to Lemma 6.2. But the module MARK in Algorithm 3.1 selects a subset M k ⊂ T k with minimal cardinality such that the same property holds, which, together with Lemma 3.7 in [7], implies Therefore, the assertion readily follows from (6.7) and (6.8).
Now we can establish the quasi-optimality of Algorithm 3.1.
Theorem 6.1. Let condition (6.4) hold. Let u be the solution to problem (2.1) and {T k , V k , u k } be the sequence of meshes, finite element spaces and discrete solutions generated by Algorithm 3.1. If (u) ∈ A s , then there holds where C qopt depends on C T 0 , m, µ and β in Theorem 5.1 but is independent of s or u.
Proof. Let M := (C qs (C rel + 1)/λC eff ) 1/2s |(u)| 1/s s . By (3.2) and Lemma 6.3, we deduce Since the oscillation term (2.11) is dominated by the global error, we can obtain from Theorem 2.2 and Lemma 4.2 that Now collecting the last three estimates, we arrive at with C θ := µ 1/2s /(1 − µ 1/2s ) bounding the geometric series. Raising this to the s-th power and noting C s give the desired estimate.

Numerical results
Now we present two numerical tests using Algorithm 3.1 with affine elements. The implementation of the algorithm is based on [13]. In the experiments, Ω is an L-shaped domain Ω = [−1, 1] 2 \([0, 1] × [−1, 0]), and σ = 1. The initial mesh T 0 is a uniform triangulation of the domain, cf. Fig. 1. At the k th adaptive iteration with triangulation T k and k = 0, 1, 2, · · · , we employ the Newton method to obtain the corresponding solution u k (of the nonlinear system). Specifically, we take the initial guess to be linear interpolation from the previous mesh, i.e., with I k : V k−1 → V k being the linear interpolation operator. The stopping criterion for the Newton iteration is u (Ω) ≤ ǫ with n denoting the Newton iteration number, and ǫ the prescribed accuracy. In the adaptive algorithm, we take ǫ = 10 −7 . Algorithm 3.1 is terminated once the sum of error indicators T ∈T k η k (u k , T ) 2 falls below a pre-specified threshold tolerance τ . We take τ = 10 −3 for both examples below. After obtaining the adaptive solution u k , we check whether the stopping condition is satisfied. If not, a refinement is carried out for those with large error indicators.
In our simulation, the Dörfler bulk criterion is used to mark elements for refinement, i.e., given θ ∈ (0, 1), we look for the minimal set M k ⊂ T k satisfying (3.1). In the experiment we show results with θ = 0.1 and θ = 0.3. To refine the mesh, we apply the newest vertex bisection (NVB) refinement [13] and bisect all edges of the elements in M k to get a finer mesh T k+1 . Alternatively, we can mark only one edge of those elements in M k for a refined mesh. Note that a smaller θ yields a more adaptive mesh and a larger iteration number.
Example 7.1. In the first example, denote Γ 0 = ∅, Γ C as the left boundary, and Γ A as the rest of the boundary. We take g(x, y) = x 2 + y 2 and f (u) = u + u 3 .
Due to the nonlinearity of the problem, the exact solution is not available. Hence, we use the solution on a very fine uniform mesh with a mesh size h = 1/2000 as the reference solution (and analogously, the Newton method is employed with a much smaller tolerance ǫ = 10 −11 ). The solution on an adaptive mesh is shown in Fig. 2(a). Since the solution singularity is localized around the re-entrant corner of the domain and the corners where the boundary condition changes, the adaptive algorithm properly refines these regions. In Fig. 3(a), we observe a convergence rate O(N −0.51 ) for the error estimator, which agrees well with the convergence rate O(N −0.54 ) in the H 1 (Ω)-norm error of the adaptive solution from Theorem 6.1, numerically verifying the reliability of the estimator. The adaptive algorithm is more efficient than the uniform refinement. Fig. 3(b) displays the convergence history with a larger parameter θ = 0.3. We obtain a smaller iteration number but a larger degrees of freedom over each refinement. The convergence rate of the error estimator and the H 1 (Ω)-error are O(N −0.50 ) and O(N −0.53 ), respectively.
In the second example, we consider the oscillatory boundary condition.     where the boundary condition changes. On a very coarse mesh, the oscillatory boundary data is not properly resolved, which leads to a slower decay at the beginning. Nonetheless, as the adaptive procedure proceeds, the convergence of the algorithm is fairly steady, with the estimator decay rate O(N −0.50 ) and the H 1 (Ω) convergence rate O(N −0.56 ) for θ = 0.1. We observe similar convergence rates for θ = 0.3 from Fig. 4(b).

Concluding remark
In this paper, for a 2D variational problem governed by a linear diffusion equation and a nonlinear boundary condition, we have analyzed an adaptive finite element method based on a residual-typed a posteriori error estimator and the Dörfler marking. We established a quasioptimal decay rate in terms of the number of elements for the algorithm, which is confirmed by the numerical experiments. One natural question is to extend the analysis to the 3D case.