An accurate a posteriori error estimator for semilinear Neumann problem and its applications

Applied Mathematics and Computation 362 (2019) 124540

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

An accurate a posteriori error estimator for semilinear Neumann problem and its applications Fei Xu∗, Qiumei Huang∗ Beijing Institute for Scientific and Engineering Computing, College of Applied Sciences, Beijing University of Technology, Beijing 100124, China

a r t i c l e

i n f o

MSC: 65N30 35J61 65M55 65B99 Keywords: Semilinear Neumann problem A posteriori error estimate Cascadic multigrid method Adaptive finite element method Complementary method

a b s t r a c t In this paper, a type of accurate a posteriori error estimator is proposed for semilinear Neumann problem, which provides an asymptotic exact estimate for the finite element approximate solution. As its applications, we design two types of cascadic adaptive finite element methods for semilinear Neumann problem based on the proposed a posteriori error estimator. The first scheme is based on the Newton iteration, which needs to solve a linearized boundary value problem by some smoothing steps on each adaptive space. The second scheme is based on the multilevel correction method, which contains some smoothing steps for a linearized boundary value problem on each adaptive space and a solving step for semilinear Neumann equation on a low dimensional space. In addition, the proposed a posteriori error estimator provides the strategy to refine mesh and control the number of smoothing steps for both of the cascadic adaptive methods. Some numerical examples are presented to validate the efficiency of the proposed algorithms in this paper. © 2019 Elsevier Inc. All rights reserved.

1. Introduction This paper is to introduce a type of asymptotic exact a posteriori error estimator for semilinear Neumann problems and study its applications in adaptive finite element method. As we know, an efficient a posteriori error estimator is vital important in scientific and engineering computing. A posteriori error estimates are capable to quantify the size of the error. In addition, they play an irreplaceable role in mesh adaptation refinement (see, e.g., [6,7,14,19]). The approximate solution of a differential equation can be derived by many numerical methods such as finite element method, finite volume method and finite difference method. In the practical simulation, the numerical solutions of a certain partial differential equation alone is not enough to guarantee the convergence and accuracy of numerical algorithm. We need to derive the corresponding error estimate. But this target is not available since we do not know the exact solution in most cases. Although we have some a priori error estimates by the finite element theory, this estimate only can be used to verify the convergence rate of the numerical algorithm. That is why an accurate a posteriori error estimate is so important. In this paper, a type of a posteriori error estimator for semilinear Neumann problems is designed based on complementary energy method from Neittaanmaki and Repin, and Vejchodský [27,33]. Further, we can prove that this estimate is asymptotic exact when the mesh size is small enough. Based on the complementary a posteriori error estimator, we design two types of cascadic adaptive finite element methods for semilinear Neumann problem due to the more and more important roles of adaptive finite element method. The ∗

Corresponding authors. E-mail addresses: [email protected] (F. Xu), [email protected] (Q. Huang).

https://doi.org/10.1016/j.amc.2019.06.054 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

2

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

adaptive finite element method (AFEM) was proposed in [6,7] by Babuška and his cooperators for the partial differential equations with singularities. As we know, for the solution with singularity, we need to use a large mount of degrees of freedom of uniform mesh to reach the desired accuracy. Since the accuracy of the smooth part can reach the desired result easily, a lot of computational work is wasted by uniform mesh. Based on the idea of error equidistribution, the adaptive finite element method was developed to generate the optimal mesh. To date, the adaptive algorithm and corresponding theoretical analysis are well-developed for both linear and nonlinear elliptic equations. For detailed results, we refer to [14,19,25,26,28,32] and the reference cited therein. The key point of adaptive finite element method lies in elements selection. Through the proposed asymptotic exact error estimator in our paper, we can generate a high-quality mesh for numerical experiment. Another efficient method will be used in our paper is multigrid method. As we know, multigrid method [11,20,36] provides an optimal algorithm for elliptic boundary value problem. There have already existed many kinds of multigrid method such as full multigrid method and cascadic multigrid method, etc. In this paper, we resort to the cascadic multigrid method due to its simplicity and efficiency. The cascadic multigrid method was proposed by Bornemann and Deuhard [10]. The main feather of cascadic multigrid method lies in correction step free. In contrast to standard multigrid methods, cascadic multigrid methods do not revisit the coarser meshes. Instead, an approximation to the finite element solution on the finest level is computed by successively interpolating and smoothing from coarse to fine. Using an elaborate stopping criterion for the employed smoother, cascadic multigrid can be shown to provide approximations, whose error in energy norm is in the order of the discretization error on the finest level. Besides, we can still derive the optimality of computational work of cascadic algorithm due to the low dimensions of coarser levels. For more information about cascadic multigrid method, please refer to [10,30,31,34] and the references cited therein. Here, we build our cascadic adaptive method for semilinear Neumann problems by combing cascadic multigrid method, adaptive method. Two kinds of cascadic adaptive schemes will be proposed. In the first scheme, we will use the Newton iteration technique [22,29,37,38] to linearize the semilinear equation on each adaptive space. In this case, the bounded second order derivatives of the nonlinear terms are needed to guarantee the convergence. More precisely, we need to do some smoothing steps for the linearized boundary value problem derived from Newton iteration on each adaptive finite element space, and the approximate solution obtained from previous level acts as the initial iteration value. The number of smoothing step depends on the error estimate of the approximate solution. But we do not know the exact error estimate in most cases. This problem appears not only in the proposed cascadic adaptive method but also in the usual numerical simulation process. Fortunately, we have an asymptotic exact error estimator for the approximate solution which will play an important role in the smoothing process. Comparing with the standard AFEM, we do not need to solve the semilinear Neumann problem directly in each adaptive space. In addition, instead of solving the linearized boundary value problem directly in each adaptive space which is adopted by the standard AFEM based on Newton iteration, we only need to do some smoothing steps. Thus the efficiency will be further improved. This will be shown in our numerical results. The second scheme is based on the cascadic multigrid method, adaptive method and recent work on multilevel correction method [21,24,35]. The main idea is to transform the solution of the semilinear Neumann problem into a series of solutions of linear boundary value problems on the sequence of adaptive finite element spaces and then correct the approximate solution by solving a semilinear Neumann problem on a specially designed low dimensional space. For the linearized boundary value problem, it is unnecessary to solve the equation exactly. In the presented scheme, we only do some smoothing steps and the number is the same as the first scheme. Due to the low dimension of correction space, in this new version of cascadic adaptive algorithm, solving semilinear Neumann problem will not be significantly more expensive than the standard adaptive finite element method for the corresponding linear boundary value problems. Similarly, the second algorithm needs less computational work than the standard AFEM and multilevel correction method. This phenomenon will also be shown in numerical results. Moreover, comparing with the first type of cascadic adaptive method for semilinear Neumann equation, the second scheme only needs the bounded first order derivative of the nonlinear term. The rest of this paper is organized as follows. In Section 2, we review the finite element method for semilinear Neumann problem. A posteriori error estimator based on complementary approach is proposed in Section 3. Two types of cascadic adaptive finite element methods for the semilinear Neumann problem are given in Sections 4 and 5, respectively. In Section 6, some numerical examples are provided to validate the efficiency of the proposed numerical methods. Some concluding remarks are given in the last section. 2. Finite element method for semilinear Neumann problem In this paper, the letter C is used to denote a mesh independent constant. The standard notations of Sobolev spaces will be used including the norm and seminorm (see, e.g., [1]). For p = 2, we denote H s () = W s,2 (). In addition, we use the symbols  · s to denote  · s,2, and V to denote H1 () for simplicity in the rest of the paper. Here, we consider the following semilinear Neumann problem:



−∇ · (A∇ u ) + u + φ (x, u ) = f, (A∇ u ) · ν = g,

in , on ∂ ,

(2.1)

where  ⊂ Rd (d = 2, 3 ) is a bounded polygonal domain, the coefficient A = (ai, j )d×d denotes a symmetric positive definite matrix with ai, j ∈ W 1,∞ (i, j = 1, 2, . . . , d ), φ (x, u) is the nonlinear term, and ν is the outward normal to ∂ .

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

3

The weak form of (2.1) can be described as: Find u ∈ V such that

a(u, v ) + (φ (x, u ), v ) = ( f, v ) + (g, v )L2 (∂ ) ,

∀v ∈ V,

(2.2)

where

a(u, v ) := (A∇ u, ∇ v ) + (u, v ).

(2.3)

Obviously, a(u, v ) satisfies

a(u, v ) ≤ Ca u1, v1,

and

ca u21, ≤ a(u, u ),



∀u, v ∈ V.

(2.4)

Thus we can define the norm wa := a(w, w ) through the bilinear form. Further, we need to give some assumptions for the nonlinear term to guarantee the existence and uniqueness of the semilinear Eq. (2.2), which is presented in the following formula. Assumption A: The nonlinear term φ (x, · ) has a nonnegative derivative with respect to the second argument

0≤

∂φ (x, v ) ≤ Cφ , ∂v

∀x ∈  and ∀v ∈ V.

(2.5)

It is time for us to introduce the finite element method for semilinear Neumann problem (2.2). In the first step, we need to generate a conforming partition of computing domain . Then based on the partition, we can construct the corresponding conforming mesh Tk . Let hT denote the diameter of T for each element T ∈ Tk and hk := max{hT : T ∈ Tk }. Let Vk ⊂ V be the finite element space of continuous piecewise polynomials over the mesh Tk . The standard finite element scheme for semilinear Neumann problem (2.2) is: Find u¯ k ∈ Vk such that

a(u¯ k , vk ) + (φ (x, u¯ k ), vk ) = ( f, vk ) + (g, vk )L2 (∂ ) ,

∀vk ∈ Vk .

(2.6)

In order to deduce the global a priori error estimates, we define ηa (Vk , ) as

ηa (Vk , ) =

inf T1 f − va ,

sup

f ∈L2 (), f 0, =1 v∈Vk

(2.7)

where the operator T1 : L2 () → V is defined as

∀ f ∈ L2 () and ∀v ∈ V.

a(T1 f, v ) = ( f, v ),

(2.8)

Similarly, we define ηa (blueVk , ∂ ) as

ηa (Vk , ∂ ) =

sup

inf T2 g − va ,

g∈L2 (∂ ),g0,∂  =1 v∈Vk

(2.9)

where the operator T2 : L2 (∂ ) → V is defined as

a(T2 g, v ) = (g, v )L2 (∂ ) ,

∀g ∈ L2 (∂ ) and ∀v ∈ V.

(2.10)

It is easy to know that ηa (Vk , ) → 0 and ηa (Vk , ∂ ) → 0 as hk → 0 (cf. [8,12,16]). In order to describe the error estimate for the finite element approximations, we need to denote

δk (u ) = inf u − vk a . vk ∈Vk

From [29], we can give the following error estimates. Lemma 2.1. When Assumption A is satisfied, Eqs. (2.2) and (2.6) are uniquely solvable and the following estimates hold

u − u¯ k a  δk (u ),

(2.11)

u − u¯ k 0  ηa (Vk , )u − u¯ k a ,

(2.12)

u − u¯ k 0,∂   (ηa (Vk , ) + ηa (Vk , ∂ ))u − u¯ k a .

(2.13)

Proof. From Theorem 6.1 in [29], we can know that problems (2.2) and (2.6) are uniquely solvable. Now, it is time to prove the error estimates. To achieve the desired goal, we define the finite element projection operator Pk : V → Vk as follows

a(Pk w − w, vk )= 0,

∀w ∈ V, ∀vk ∈ Vk .

It is easy to know that

u − Pk ua = δk (u ),

(2.14)

u − Pk u0  ηa (Vk , )u − Pk ua ,

(2.15)

4

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

u − Pk u0,∂   ηa (Vk , ∂ )u − Pk ua .

(2.16)

Denoted by wk = Pk u − u¯ k . From (2.2), (2.5) and (2.6), we have

a(Pk u − u¯ k , wk ) ≤ a(Pk u − u¯ k , wk ) + (φ (x, Pk u ) − φ (x, u¯ k ), wk ) = a(Pk u, wk ) + (φ (x, Pk u ), wk ) − ( f, wk ) − (g, wk )L2 (∂ ) = a(Pk u − u, wk ) + (φ (x, Pk u ) − φ (x, u ), wk ) = (φ (x, Pk u ) − φ (x, u ), wk ) ≤ Cφ u − Pk u0 wk a . Then the following inequalities hold

Pk u − u¯ k a ≤ Cφ u − Pk u0 ≤ Cφ ηa (Vh , )u − Pk ua .

(2.17)

Combining (2.17) and triangle inequality leads to the following estimates

u − u¯ k a ≤ u − Pk ua + Pk u − u¯ k a ≤ δk (u ) + Cφ ηa (Vk , )u − Pk ua ≤ (1 + Cφ ηa (Vk , ))δk (u ),

(2.18)

which is the desired result (2.11). From (2.17) and triangle inequality, we have

u − u¯ k 0 ≤ u − Pk u0 + Pk u − u¯ k 0 ≤ u − Pk u0 + C Pk u − u¯ k a ≤ C ηa (Vk , )u − Pk ua + Cφ ηa (Vk , )u − Pk ua ≤ (C + Cφ )ηa (Vk , )u − Pk ua ≤ (C + Cφ )ηa (Vk , )u − u¯ k a and

u − u¯ k 0,∂  ≤ u − Pk u0,∂  + Pk u − u¯ k 0,∂  ≤ u − Pk u0,∂  + C Pk u − u¯ k a ≤ C ηa (Vk , ∂ )u − Pk ua + Cφ ηa (Vk , )u − Pk ua ≤ C (ηa (Vk , ∂ ) + ηa (Vk , ))u − u¯ k a , which is the desired result (2.12) and (2.13). The proof is completed.



3. Complementarity based error estimate In this section, we will propose a type of a posteriori error estimate for the finite element approximate solution. First, we recall the following Green’s theorem. Lemma 3.1. Let  ⊂ Rd be a bounded Lipschitz domain with unit outward normal ν to ∂ . Then the following Green’s formula holds





vdivpd +





p · ∇vd =



∂

v p · ν d ∂ ,

∀v ∈ H 1 () and ∀p ∈ W,

(3.1)

where W := {p ∈ (L2 ())d : divp ∈ L2 (), p · ν ∈ L2 (∂ )}. For the finite element approximate solution u¯ k , we have the following estimate. Theorem 3.1. Assume the mesh size hk is small enough, we have the following estimate for the finite element approximate solution u¯ k

u − u¯ k a ≤

1 + C (ηa (Vk , ) + ηa (Vk , ∂ )) min η (u¯ k , p ), p∈W 1 − C ηa (Vk , )

(3.2)

and a posteriori error estimator η (u¯ k , p ) is defined as

 1 η (u¯ k , p ) :=  f − φ (x, u¯ k ) − u¯ k + divp|20 + A−1/2 p − A1/2 ∇ u¯ k 20 + g − p · ν20,∂  2 .

Proof. From (2.2), (2.6), (3.1) and Lemma 2.1, for any p ∈ W, we have

a(u − u¯ k , w ) = a(u, w ) − a(u¯ k , w ) = ( f, w ) + (g, w )L2 (∂ ) − (φ (x, u ), w ) − (A∇ u¯ k , ∇ w ) − (u¯ k , w ) = ( f, w ) + (g, w )L2 (∂ ) − (φ (x, u ), w ) − (A∇ u¯ k , ∇ w ) − (u¯ k , w ) +(divp, w ) + (p, ∇ w ) − (p · ν, w )L2 (∂ )

= ( f − φ (x, u¯ k ) − u¯ k + divp, w ) + (p − A∇ u¯ k , ∇ w ) + (g − p · ν, w )L2 (∂ )

(3.3)

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

5

+(φ (x, u¯ k ) − φ (x, u ), w )

 f − φ (x, u¯ k ) − u¯ k + divp0 w0 + A−1/2 p − A1/2 ∇ u¯ k 0 A1/2 ∇ w0 +g − p · ν0,∂  w0,∂  + C ηa (Vk , )u − u¯ k a wa  1/2 ≤  f − φ (x, u¯ k ) − u¯ k + divp20 + A−1/2 p − A1/2 ∇ u¯ k 20 + g − p · ν20,∂   1/2 × w20 + A1/2 ∇ w20 + w20,∂  + C ηa (Vk , )u − u¯ k a wa , ≤

(3.4)

where the Hölder inequality is used. Take w = u − u¯ k , we can derive

  u − u¯ k 2a ≤ η (u¯ k , p ) 1 + C (ηa (Vk , ) + ηa (Vk , ∂ )) u − u¯ k a + C ηa (Vk , )u − u¯ k 2a .

The desired result (3.2) immediately follows by the arbitrariness of p ∈ W.



Now, it is natural for us to discuss the optimization problem in (3.2): Find p∗ ∈ W such that

η (u¯ k , p∗ ) = min η (u¯ k , p ).

(3.5)

p∈W

In order to derive the solution of the optimization problem (3.5), we can transform the optimization problem into a partial differential equation, which can be solved by the standard finite element method. The corresponding results are presented in the following theorem. Theorem 3.2. The optimization problem (3.5) is equivalent to the following partial differential equation: Find p∗ ∈ W such that

a∗ ( p∗ , q ) = F ( q ), where

∀q ∈ W,

(3.6)

  a∗ (p, q ) := (divpdivq + A−1 pq )d + (p · ν )(q · ν )d∂ , ∂    F (q ) := (φ (x, u¯ k ) − f )divqd + (g + u¯ k )(q · ν )d∂ . 

∂

Proof. For any t ∈ R and q ∈ W, there holds

η2 (u¯ k , p∗ + tq ) ≥ η2 (u¯ k , p∗ ).

(3.7)

By expanding a posteriori error estimator according to the definition (3.3), we can simplify (3.7) into the following inequality

at 2 + 2bt ≥ 0,

(3.8)

where

a = divq20 + A−1/2 q20 + q · ν20,∂  and

b = ( f − φ (x, u¯ k ) − u¯ k + divp, divq ) + (A−1/2 p − A1/2 ∇ u¯ k , A−1/2 q ) + (p · ν − g, q · ν )L2 (∂ ) . From (3.7), we have b = 0, which leads to the desired result (3.6).



According to (3.3) and Theorem 3.2, we can get the following relationship. Theorem 3.3. Assume p∗ be the solution of problem (3.6). Then the following equality holds

η2 (u¯ k , p ) = η2 (u¯ k , p∗ ) + p∗ − p2W + (p − p∗ ) · ν20,∂  , ∀p ∈ W.

(3.9)

Proof. From the definition of a posteriori estimator (3.3) and dual problem (3.6), we can derive

η2 (u¯ k , p ) = η2 (u¯ k , p∗ + p − p∗ ) = η2 (u¯ k , p∗ ) + div(p − p∗ )20 + A−1/2 (p − p∗ )20 + (p − p∗ ) · ν20,∂  +2( f − φ (x, u¯ k ) − u¯ k + divp∗ , div(p − p∗ )) +2(A−1/2 p − A1/2 ∇ u¯ k , A−1/2 (p − p∗ )) + 2(p∗ · ν − g, (p − p∗ ) · ν ). Combing (3.6) and (3.10) leads to the desired result (3.9).

(3.10) 

Although we transform the optimization problem (3.5) into the partial differential equation (3.6), we can only derive an approximate solution since the exact solution is unknown. Once we get the approximations, we actually derive a computable error estimate which is shown in the following corollary.

6

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Corollary 3.1. Under the assumption of Theorem 3.1 and suppose pk ∈ W is an approximate solution for (3.6), then we get the following computable estimate

u − u¯ k a ≤

1 + C (ηa (Vk , ) + ηa (Vk , ∂ )) η (u¯ k , pk ). 1 − C ηa (Vk , )

(3.11)

Now, it is time for us to concentrate on the efficiency and reliability of a posteriori error estimator (3.3). Theorem 3.4. Let u¯ k be the approximate solution of the discrete problem (2.6) and p∗ be the solution of (3.6). Under the conditions of Theorem 3.1, we have

θ1 u − u¯ k a ≤ η (u¯ k , p∗ ) ≤ θ2 u − u¯ k a ,

(3.12)

where

θ1 :=

1 − C ηa (Vk , ) 1 + C (ηa (Vk , ) + ηa (Vk , ∂ ))

and

θ2 :=



1 + C ηa2 (Vk , ).

Furthermore, we have the following asymptotic exactness

η (u¯ k , p∗ ) = 1. ¯ k a hk →0 u − u lim

(3.13)

Proof. The left-hand side inequality of (3.12) can be obtained from (3.2) directly. Now, we prove the right-hand side one of (3.12). From (2.2), (3.3) and the fact A∇ u ∈ W, we have

η2 (u¯ k , A∇ u ) =  f − φ (x, u¯ k ) − u¯ k + ∇ · (A∇ u )20 + A1/2 ∇ (u − u¯ k )20 + g − A∇ u · ν20,∂  . =

φ (x, u ) − φ (x, u¯ k ) + u − u¯ k 20 + A1/2 ∇ (u − u¯ k )20 .

From (3.5), we can derive

η2 (u¯ k , p∗ ) ≤ η2 (u¯ k , A∇ u ) ≤ (1 + Cφ )2 u − u¯ k 20 + A1/2 ∇ (u − u¯ k )20 ≤ u − u¯ k 2a + (Cφ2 + 2Cφ )ηa2 (Vk , )u − u¯ k 2a ≤ u − u¯ k 2a + C ηa2 (Vk , )u − u¯ k 2a , which is the right-hand side inequality of (3.12). The asymptotically exact error estimate (3.13) can be deduced easily from the fact that ηa (Vk , ) → 0 as hk → 0.  Remark 3.1. If we solve the dual problem (3.6) approximately by high order finite element method to derive an approximate solution pk , which satisfies p∗ − pk W + (pk − p∗ ) · ν20,∂  ≤ γ u − u¯ k a for some constant γ > 0. Then the following error estimate holds

η (u¯ k , pk ) ≤



θ22 + γ 2 u − u¯ k a .

(3.14)

Furthermore, the error estimator η (u¯ k , pk ) is asymptotically exact if the following condition holds

p∗ − pk W + (pk − p∗ ) · ν20,∂  = 0. u − u¯ k a hk →0 lim

(3.15)

Proof. From (3.9) and (3.12), we have

η2 (u¯ k , pk ) = η2 (u¯ k , p∗ ) + pk − p∗ 2W + (pk − p∗ ) · ν20,∂  θ22 u − u¯ k 2a + γ 2 u − u¯ k 2a ≤ (θ22 + γ 2 )u − u¯ k 2a , ≤

(3.16)

which is the desired result (3.14) and the asymptotical exactness of the estimator follows immediately from the condition (3.15). 4. Cascadic adaptive finite element method based on Newton iteration In this section, we give a type of cascadic adaptive finite element method based on Newton iteration and complementarity based error estimate. As we know, for the cascadic multigrid finite element method, we need to do some smoothing steps on the sequence of finite element spaces and the number of smoothing steps increase gradually on coarser levels. Following the idea of cascadic multigrid method, we propose a cascadic adaptive finite element method for semilinear Neumann problem in this section.

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

7

First, we introduce the smoothing operator Sk : Vk → Vk of cascadic algorithm which satisfies

⎧ ⎨

Skm wk a ≤ mCα h1k wk 0 , Skm wk a ≤ wk a , ⎩ m Sk (wk + vk )a ≤ Skm wk a + Skm vk a ,

(4.1)

where the constant α is dependent on the type of smoothing operator. For the Richardson iteration and symmetric Gauss– Seidel, we have α = 1/2. For the conjugate-gradient iteration, we have α = 1 (cf. [30,31]). Then we use the following equation

wk = Smooth(Vk , f, g, ξk , m, Sk )

(4.2)

to represent the smoothing operate for the following linearized elliptic problem derived from Newton iteration

a(uk , vk ) + (φ (x, ξk )uk , vk ) = ( f, vk ) + (g, vk )L2 (∂ ) ,

∀vk ∈ Vk ,

(4.3)

where wk is the final result of the smoothing operate, ξ k is the iteration initial value, Sk is the smoothing operator, m is the iteration number, φ is the derivative of nonlinear term φ with respect to the second argument. 4.1. Cascadic adaptive finite element method The standard adaptive finite element method can be written as loop of the form

Solve → Estimate → Mark → Refine. In order to refine mesh, we need to use a local error estimator. The complementary a posteriori error estimator η(uk , pk ) is a competitive candidate and will be used to describe the algorithms in this paper. Besides, other efficient a posteriori error estimators such as the residual type a posteriori error estimator (see (6.7)) can also be adopted. In our numerical experiments, we have compared the results of the complementary a posteriori error estimator and residual type a posteriori error estimator in mesh refinement, and the former works better. In order to simplify the notation, ηT (uk , pk ) := ( f − φ (x, uk ) − uk + divpk 20,T + A−1/2 pk − A1/2 ∇ uk 20,T + g − pk · ν20,∂ ∩∂ T )1/2 is used to denote the restriction of η(uk , pk ) on each element T ∈ Tk . Based on the local error indicator, we use the Dörfler’s Marking Strategy [19] presented as follows to mark the elements for local refinement.

Dörfler’s Marking Strategy 1. Given a parameter θ ∈ (0, 1). 2. Construct a minimal subset Mk from Tk by selecting some elements in Tk such that



ηT2 (uk , pk ) ≥ θ η2 (uk , pk ).

T ∈Mk

3. Mark all the elements in Mk .

Based on the Dörfler’s Marking Strategy, the first type of cascadic adaptive finite element method is defined as follows.

Algorithm 1. Cascadic adaptive finite element method. 1. Generate an initial mesh T1 for the computing domain . Then build the initial finite element space V1 on the triangulation T1 and solve the following semilinear Neumann equation: Find u1 ∈ V1 such that

a(u1 , v1 ) + (φ (x, u1 ), v1 ) = ( f, v1 ) + (g, v1 )L2 (∂ ) ,

∀v1 ∈ V1 .

(4.4)

2. Set k = 1. 3. Compute the local error indicator ηT (uk , pk ) on each element T ∈ Tk . 4. Construct the submesh Mk ⊂ Tk by Dörfler’s Marking Strategy and refine Tk to generate a new conforming mesh Tk+1 , then construct the corresponding space Vk+1 .

8

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

5. Derive the new approximate solution uk+1 ∈ Vk+1 by the following smoothing step

uk+1 = Smooth(Vk+1 , f − φ (x, uk ) + φ (x, uk )uk , g, uk , mk+1 , Sk+1 ).

(4.5)

6. Let k = k + 1 and go to step 3.

Since only smoothing steps are needed on each mesh level, the key point of the cascadic adaptive finite element method lies in the adaption of smoothing number mk on each coarse space such that the final error of approximate solution on fine space has already been eliminated on the coarse space. In the standard cascadic multigrid method, this goal is achieved by increasing the number of smoothing steps gradually on coarser levels. Requiring the number of smoothing operations is proportional to the number of unknowns on the final space, the error of the final approximation is of the same order as the error of standard finite element method. The normal choice of the iteration number mk on level k is

mk = m

 n d1α

nk

,

(4.6)

where nk = dimVk , d denotes the space dimension and denotes the index of the final level. The choice of (4.6) comes from the theoretical analysis of standard cascadic multigrid method, which is used to derive the convergence property. For more details, please refer to [10,18]. However, at the intermediate level k, we have not derived the number n yet for the adaptive mesh, which means that the strategy (4.6) is not implementable. To make it implementable, we define the final level as the first level on which the approximate error is below some user given tolerance TOL. Hence we have the following relation

u − uk a TOL

≈

 n 1d

nk

,

(4.7)

which leads us to replace (4.6) by

mk = m

u − uk a

 α1 .

TOL

(4.8)

Herein the approximate error u − uk is not known since the approximate solution uk is not derived yet, but can be replaced by the estimate

u − uk a ≈ u − uk−1 a

n

k−1

1d

nk

.

(4.9)

For the linear boundary value problems, there exists some methods to derive the discretization error u − uk−1 a (see, e.g., [9,17]). But unfortunately, these methods do not work for the nonlinear case. Here, we use the complementarity based error estimate η (uk−1 , pk−1 ) to control u − uk−1 a since a posteriori error estimator η (uk−1 , pk−1 ) is an excellent error estimate for u − uk−1 a . Hence, the number of smoothing steps are chosen by

mk = m

η (uk−1 , pk−1 )  nk−1 d

1

TOL

 α1 ,

nk

(4.10)

and the good performance of the scheme is examined by a variety of numerical examples. Now we come to estimate the computational work for Algorithm 1. Here we have to use additionally, that the sequence of unknowns belongs to a geometric progression:

nk < σ0 nk ≤ nk+1 < σ1 nk ,

k = 1, 2, . . .

(4.11)

Theorem 4.1. Assume the computational work of solving  process for the semilinear Neumann problem in the initial space V1 is M1 . Then the computational work of Algorithm 1 is O M1 + n ) if d > 1/α , and the computational work is O(n ) if M1 is small









enough. Further, the computational work will change to O M1 + n log(n ) and O n log(n ) if d = 1/α . Proof. Let wk denote the computational work on the space Vk , then we have w1 = M1 and wk = mk nk , for k = 2, . . . , . Then from Algorithm 1, (4.6) and (4.11), there holds

Computational Work =





w k = O M1 +

k=1



= O M1 + m n



mk nk

k=2



1 k=2

σ0



( −k)(1− d1α ) 

.

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540





9



Thus, the total computational work is O M1 + n ) if d > 1/α and is O M1 + n log(n ) if d = 1/α , and moreover, the total computational work changes to O (n ) and O (n log(n )) if M1 is small enough.  5. Cascadic adaptive finite element method based on multilevel correction scheme In this section, we give another type of cascadic adaptive finite element method based on the multilevel correction scheme and complementarity based error estimate. Similarly, we use the following equation

wk = Smooth(Vk , f, g, ξk , m, Sk )

(5.1)

to represent the smoothing operate for the following linear equation

∀vk ∈ Vk ,

a(uk , vk ) = ( f, vk ) + (g, vk )L2 (∂ ) ,

(5.2)

where wk is the final result of the smoothing operate, ξ k is the iteration initial value, Sk is the smoothing operator, m is the iteration number. 5.1. One correction step In order to describe the cascadic adaptive algorithm clearly, a type of correction step will be introduced first in this subsection. This correction step describes the way to refine the mesh and derive the corresponding approximate solution. First, assume that we have derived an approximate solution uk on the current finite element space Vk . Then we need to refine the current mesh Tk in the adaptive way based on a posteriori error estimator and Dörfler’s Marking Strategy (cf. [19]). Let us denote the new mesh and finite element space by Tk+1 and Vk+1 . The correction process includes two steps. In the first step, we only need to solve a linearized problem by inserting the approximate solution from the previous level into the nonlinear term. Since there exists many mature numerical methods for the linear problem, the efficiency will be improved compared with the nonlinear equation which need to be solved by iteration computation. In the second step, we correct the approximate solution of the linearized problem by solving a low dimensional semilinear Neumann problem on a special space VH,k . The special space VH,k is constructed by combing a coarse space VH and approximate solution of the linearized problem. The details of the correction step are presented in the following algorithm.

Algorithm 2. One correction step. 1. Compute the local error indicator ηT (uk , pk ) on each element T ∈ Tk . Construct the submesh Mk ⊂ Tk by Dörfler’s Marking Strategy and refine Tk to generate a new conforming mesh Tk+1 , then construct the corresponding finite element space Vk+1 . 2. Define a linearized boundary value problem in the space Vk+1 : Find  uk+1 ∈ Vk+1 such that

a ( uk+1 , vk+1 ) = ( f, vk+1 ) + (g, vk+1 )L2 (∂ ) − (φ (x, uk ), vk+1 ),

∀vk+1 ∈ Vk+1 .

(5.3)

Perform the smoothing process (5.1) to the equation (5.3) and obtain a new approximate solution  uk+1 ∈ Vk+1 by

 uk+1 = Smooth(Vk+1 , f − φ (x, uk ), g, uk , mk+1 , Sk+1 ).

(5.4)

3. Define an enriched finite element space VH,k+1 = VH + span{ uk+1 } and solve the following semilinear Neumann equation: Find uk+1 ∈ VH,k+1 such that

a(uk+1 , vk+1 ) + (φ (x, uk+1 ), vk+1 ) = ( f, vk+1 ) + (g, vk+1 )L2 (∂ ) ,

∀vk+1 ∈ VH,k+1 .

In order to simplify the notation, we summarize above steps into

uk+1 = SmoothCorrection(VH , uk , mk+1 , Sk+1 , Vk+1 ).

The number of smoothing step in Algorithm 2 is chosen by the same way as Algorithm 1:

mk = m

 η (u

k−1 , pk−1

TOL

)  nk−1  1d α 1

nk

.

(5.5)

Remark 5.1. We want to emphasize that the dimension of the enriched space VH,k remains unchange during the adaptive process (In fact, dimVH,k = dimVH +1). So we only need to solve a low dimensional semilinear Neumann problem in the third step of Algorithm 2 and the computational time can be ignored with the refinement of mesh.

10

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

5.2. Cascadic adaptive finite element method In this subsection, we will design a type of cascadic adaptive finite element method based on one correction step defined in Algorithm 2. First, we need to generate a coarse mesh TH on the computing domain and construct the corresponding finite element space VH . Then we generate an initial mesh T1 by refining TH several times by uniform way and construct the initial space V1 . On the initial space V1 , we solve the semilinear Neumann problem and derive the approximate solution u1 . Then next we execute the one correction step defined in Algorithm 2 to generate a series of adaptive solutions. Based on the cascadic adaptive method described above, and one correction step defined by Algorithm 2, the second type of cascadic adaptive method is given below.

Algorithm 3. Cascadic adaptive finite element method. 1. Given a mesh refinement parameter 0 < θ < 1. Construct a coarse space VH and the initial space V1 . Then solve the semilinear Neumann equation in the initial space: Find u1 ∈ V1 such that

a(u1 , v1 ) + (φ (x, u1 ), v1 ) = ( f, v1 ) + (g, v1 )L2 (∂ ) ,

∀v1 ∈ V1 .

2. Set k = 1. 3. Obtain a new approximate solution uk+1 ∈ Vk+1 by Algorithm 2:

uk+1 = SmoothCorrection(VH , uk , mk+1 , Sk+1 , Vk+1 ). 4. Let k = k + 1 and go to step 3.

Remark 5.2. In the standard adaptive finite element method for semilinear Neumann problem, the semilinear equation is solved on each adaptive space Vk by iteration method. Here, we can regard Algorithm 3 as the adaptive algorithm on a series of spaces VH,k . The special space VH,k is constructed by inserting an approximate solution  uk into the coarse space VH . That is, the solution of the linearized boundary value problem is regarded as a new basis function of the enriched space. Since the accuracy of approximate solution  uk improves gradually, the space VH,k will become more and more accurate to approximate the semilinear Neumann problem. Meanwhile, the dimension of VH,k is small and fixed during the adaptive process and the efficiency of Algorithm 3 would be benefit from it. Now we come to estimate the computational work for Algorithm 3. Theorem 5.1. Assume the computational work of solving process for the semilinear Neumann problem in the coarse space  VH and initial space V1 are MH and M1 . Then the computational work of Algorithm 3 is O M1 + MH log(n ) + n ) if d > 1/α , and  the computational work is O(n ) if MH and M1 are small enough. Further, the computational work will change to O M1 + MH log(n ) + n log(n ) and O n log(n ) if d = 1/α . Proof. Let wk denote the computational work on the space Vk , then we have w1 = M1 and wk = mk nk + MH , 2, . . . , . Then from Algorithm 3, (4.6) and (4.11), there holds

Computational work =





w k = O M1 +

k=1



= O M1 +





k=



( mk n k + MH )

k=2

for



mk n k + MH ( − 1 )

k=2

= O M1 + MH log(n ) + m n





1 k=2

σ0

( −k)(1− d1α )

 .





Thus, the total computational work is O M1 + MH log(n ) + n ) if d > 1/α and is O M1 + MH log(n ) + n log(n ) if d = 1/α , and moreover, the total computational work changes to O (n ) and O (n log(n )) if MH and M1 are small enough.  6. Numerical results In this section, we present four examples to illustrate the efficiency of the cascadic adaptive algorithms. Different examples have different nonlinear terms including polynomial function, exponential function and the nonlinear term that has no

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann problem

11

Errors for semilinear elliptic equation

η (u ,Ω)

η (u ,Ω)

η(uk,p0k )

η(u ,p )

k

k

k

||u−u ||

k a

slope=−1/2

0

10

Errors

Errors

10

1 k

||u−u ||

k a

slope=−1/2

0

k

k

−1

−1

10

10

2

3

10

10

4

10

5

2

10

10

Number of elements

3

10

4

10

5

10

Number of elements

Fig. 1. Errors for Example 1 when the semilinear Neumann problem (6.2) is solved by the linear finite element method. η (uk , p0k ) and η (uk , p1k ) denote η(uk , pk ) when the dual problem is solved by Wk0 and Wk1 . ηk (uk , ) denotes the residual type a posteriori error estimator.

bounded second order derivatives. Here, the numerical examples are implemented on the machine PowerEdge R720 with the linux system, and the memory is 24G. In all the benchmark examples, we choose conjugate gradient iteration as the smoother. And the number of smoothing times are settled by (4.10) with m = 2. In order to compute the complementary based a posteriori error estimator, we need to solve the dual problem (3.6). Here, the dual problem (3.6) and semilinear Neumann problem (2.2) are solved using the same mesh. The corresponding p finite element space for dual problem is the H(div, ) space Wk defined as follows [13]:

Wkp = {w ∈ W : w|K ∈ RTp , ∀K ∈ Tk }, where RTp = (P p

)d

(6.1)

+ xP˜p .

6.1. Example 1 In the first example, we will show the asymptotic exactness of the complementary error estimate. Let us consider the following semilinear Neumann problem:



−u + u + u3 = f, ∇ u · ν = g,

in , on ∂ ,

(6.2)

where  = (0, 1 )2 . We choose the right hand side term f and Neumann boundary condition g such that the exact solution is given by

u = sin(π x ) cos(π y ).

(6.3)

In order to verify the asymptotic exactness of the complementary a posteriori error estimator, we solve semilinear Neumann equation (6.2) and corresponding dual problem on a series of uniform mesh sequence. In our numerical experiment, linear finite element space is constructed for the semilinear Neumann problem (6.2), Wk0 and Wk1 are constructed for dual problem (3.6). We can find the corresponding numerical results in Fig. 1, from which we can see that η(uk , pk ) is asymptotic exact when we solve the dual problem in the space Wk1 . Moreover, we also compare the results between the a complementary error estimator η(uk , pk ) and residual type a posteriori error estimator defined as follows. First, we construct the element residual RT (uk ) and the jump residual Je (uk ) for the approximate solution uk as follows:

RT (uk ) := −∇ · (A∇ uk ) + uk + φ (x, uk ) − f,

 Je (uk ) :=

1 2

( A∇ u k | T + · ν + + A∇ u k | T − · ν − ) , A∇ uk · ν − g,

in T ∈ Tk , for for

e ∈ Ek , e ∈ E∂  ,

(6.4)

(6.5)

where e is the common side of elements T + and T − with the unit outward normals ν + and ν − , respectively. Let Ek be the set of interior faces of Tk , and E∂  be the set of boundary faces. For T ∈ Tk , we define the local error indicator ηk2 (uk , T )

12

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

Mesh after 10 iterations

η (u ,Ω) k

k

η(u ,p1) k

k

||u−u ||

k a

0

slope=−1/2

Errors

10

−1

10

−2

10

2

3

10

4

10

5

10

10

Number of elements Fig. 2. Left: Errors of Algorithm 1 for Example 2 when the complementary a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

by

ηk2 (uk , T ) := h2T RT (uk )20,T +

e∈Ek ,e⊂∂ T

he Je (uk )20,e .

(6.6)

Given a subset ω ⊂ , we define the error estimator ηk2 (uh , ω ) by

ηk2 (uk , ω ) =



T ∈Tk ,T ⊂ω

ηk2 (uk , T ).

(6.7)

The corresponding numerical results of the residual type a posteriori error estimator are also presented in Fig. 1 which show that the complementary a posteriori error estimator has a better efficiency than the residual type a posteriori error estimator. 6.2. Example 2 In the second example, we solve the following semilinear Neumann problem:



−u + u + u3 = 4, ∇ u · ν = 1,

in , on ∂ ,

(6.8)

where  = (−1, 1 )2 /[0, 1]2 . Due to the reentrant corner of , the exact solution with singularity is expected. We give the numerical results for the approximate solutions by Algorithms 1 and 3 with the same parameter θ = 0.4. The number of smoothing step is chosen as (4.10). In each algorithm, both the residual type a posteriori error estimator and the complementary a posteriori error estimator are used in mesh refinement to test the efficiency. The linear conforming element space is used to solve the semilinear elliptic problem (6.8) and finite element space Wk1 is used to solve the dual problem (3.6). Since the exact solution is not known, we choose an adequately accurate approximation on a fine enough mesh as the exact solution. The corresponding numerical results are presented in Figs. 2–5. From Figs. 2 to 5, we can find that both of the two proposed cascadic adaptive finite element methods are able to obtain the optimal accuracy. In order to show the efficiency of the complementary a posteriori error estimator in mesh refinement, we compare the computational time of Algorithms 1 and 3 with the complementary a posteriori error estimator and residual type a posteriori error estimator for mesh refinement. In our test, we set the same accuracy of tolerance for all examples (TOL=1E2). The CPU time (in seconds) is provided in Table 1 and we can find that the complementary a posteriori error estimator has a better efficiency than the residual type a posteriori error estimator for both of the two algorithms. Besides, we can find that Algorithm 1 has a better efficiency than Algorithm 3 when the nonlinear term has a bounded second order derivative. In order to further show the efficiency of Algorithms 1 and 3, we propose the computational time of direct AFEM (i.e., solve the semilinear Neumann problem directly on each adaptive space), AFEMs based on Newton iteration and multilevel correction method (i.e., solve the linearized boundary value problem involved in Algorithms 1 and 3 directly by conjugate gradient method rather than executing smoothing steps). The corresponding computational times are presented in Table 2. From Table 2, we can find that our algorithms have a better efficiency.

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

13

Mesh after 10 iterations

η (u ,Ω) k

k

η(u ,p1) k

k

||u−u ||

k a

0

slope=−1/2

Errors

10

−1

10

−2

10

2

10

3

4

10

5

10

10

Number of elements Fig. 3. Left: Errors of Algorithm 1 for Example 2 when the residual type a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Errors for semilinear Neumann equation

Mesh after 10 iterations

η (u ,Ω) k

k

η(u ,p1) k

k

||u−u ||

k a

0

slope=−1/2

Errors

10

−1

10

−2

10

2

10

3

10

4

10

5

10

Number of elements Fig. 4. Left: Errors of Algorithm 3 for Example 2 when the complementary a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements. Table 1 The computational time (in seconds) of the proposed cascadic adaptive algorithms with the same tolerance TOL for Example 2. Algorithm Time Algorithm 1 with complementary a posteriori error estimator

Computational time 529.84

Error estimate 9.37E-3

Algorithm 1 with residual type a posteriori error estimator Algorithm 3 with complementary a posteriori error estimator Algorithm 3 with residual type a posteriori error estimator

1936.55 572.90 2154.75

9.10E−3 9.34E−3 9.12E−3

Table 2 The computational time (in seconds) of the proposed cascadic adaptive algorithms, direct AFEM, AFEM based on Newton iteration and multilevel correction method for Example 2. Algorithm Time

Computational time

Error estimate

Direct AFEM Algorithm 1 with complementary a posteriori error estimator AFEM based on Newton iteration Algorithm 3 with complementary a posteriori error estimator AFEM based on multilevel correction method

2016.72 529.84 660.72 572.90 721.29

9.32E−3 9.37E−3 9.36E−3 9.34E−3 9.35E−3

14

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

Mesh after 10 iterations

ηk(uk,Ω) η(u ,p1)

0

k

Errors

10

k

slope=−1/2

−1

10

−2

10

3

10

4

5

10

10

Number of elements Fig. 5. Left: Errors of Algorithm 3 for Example 2 when the residual type a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Errors for semilinear Neumann equation

Mesh after 10 iterations

η (u ,Ω) k

k

1

η(uk,pk )

0

Errors

10

slope=−1/2

−1

10

−2

10

2

3

10

10

4

10

5

10

Number of elements Fig. 6. Left: Errors of Algorithm 1 for Example 3 when the complementary a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

6.3. Example 3 In the third example, we consider the following semilinear Neumann problem:



−∇ · (A∇ u ) + u − e−u = 1, ( A∇ u ) · ν = 1 ,



in , on ∂ ,

(6.9)

2 0 ,  = (−1, 1 )2 /[0, 1]2 . Due to the reentrant corner of , the exact solution with singularity is expected. 0 9 We give the numerical results for the approximate solutions by Algorithms 1 and 3 with the same parameter θ = 0.4. The number of smoothing step is chosen as (4.10). In each algorithm, both the residual type a posteriori error estimator and the complementary a posteriori error estimator are used in mesh refinement to test the efficiency. The linear conforming element space is used to solve the semilinear elliptic problem (6.9) and finite element space Wk1 is used to solve the dual problem (3.6). The corresponding numerical results are presented in Figs. 6–9. From Figs. 6 to 9, we can also find that the proposed two types of cascadic adaptive finite element methods are able to obtain the optimal accuracy. where A =

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

15

Mesh after 10 iterations

η (u ,Ω) k

k

10

Errors

k

η(u ,p1)

0

k

slope=−1/2

−1

10

−2

10

3

4

10

5

10

10

Number of elements Fig. 7. Left: Errors of Algorithm 1 for Example 3 when the residual type a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Errors for semilinear Neumann equation

Mesh after 10 iterations

η (u ,Ω) k

k

10

Errors

k

η(u ,p1)

0

k

slope=−1/2

−1

10

−2

10

3

4

10

10

5

10

Number of elements Fig. 8. Left: Errors of Algorithm 3 for Example 3 when the complementary a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Table 3 The computational time (in seconds) of the proposed cascadic adaptive algorithms with the same tolerance TOL for Example 3. Algorithm Time Algorithm Algorithm Algorithm Algorithm

1 1 3 3

with with with with

complementary a posteriori error estimator residual type a posteriori error estimator complementary a posteriori error estimator residual type a posteriori error estimator

Computational time

Error estimate

110.34 146.58 124.55 169.36

8.26E−3 8.95E−3 8.21E−3 8.97E−3

In order to show the efficiency of the complementary a posteriori error estimator in mesh refinement, we also compare the computational time of Algorithms 1 and 3 with the complementary a posteriori error estimator and residual type a posteriori error estimator for mesh refinement. In our test, we set the same accuracy of tolerance for all examples (TOL=1E−2). The CPU time (in seconds) is provided in Table 3 and we can find that the complementary a posteriori error estimator has a better efficiency than the residual type a posteriori error estimator for both of the two algorithms. Besides, we can also find that Algorithm 1 has a better efficiency than Algorithm 3 when the nonlinear term has a bounded second order derivative.

16

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

Mesh after 10 iterations

ηk(uk,Ω) η(u ,p1)

0

k

Errors

10

k

slope=−1/2

−1

10

−2

10

3

4

10

10

5

10

Number of elements Fig. 9. Left: Errors of Algorithm 3 for Example 3 when the residual type a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements. Table 4 The computational time (in seconds) of the proposed cascadic adaptive algorithms, direct AFEM, AFEM based on Newton iteration and multilevel correction method for Example 3. Algorithm time

Computational time

Error estimate

Direct AFEM Algorithm 1 with complementary a posteriori error estimator AFEM based on Newton iteration Algorithm 3 with complementary a posteriori error estimator AFEM based on multilevel correction method

458.75 110.34 138.53 124.55 159.23

8.20E−3 8.26E−3 8.25E−3 8.21E−3 8.23E−3

In order to show the efficiency of Algorithms 1 and 3, we also propose the computational time of direct AFEM, AFEMs based on Newton iteration and multilevel correction method. The corresponding computational times are presented in Table 4. From Table 4, we can also find that our algorithms have a better efficiency. 6.4. Example 4 In the last example, we solve the following semilinear Neumann problem:



−u + u + φ (x, u ) = f, ∇ u · ν = g,

with

 φ (x, u ) =

u 3/2 ,

(−u )3/2 ,

if if

in , on ∂ ,

(6.10)

u ≥ 0, u < 0,

(6.11)

where  = (−1, 1 )2 /[0, 1]2 . We choose the right hand side term f and Neumann boundary condition g such that the exact solution is given by

u = sin(10π x ) cos(10π y ).

(6.12)

As we can see from equation (6.10), the nonlinear term φ (x, u) only has bounded first order derivative and has no bounded second order derivative. Then the numerical algorithm given in Algorithm 1 can not be used for this example, which means that Algorithm 3 has a wider scope of application. We give the numerical results for the approximate solutions by Algorithm 3 with the parameter θ = 0.4. Both the residual type a posteriori error estimator and the complementary a posteriori error estimator are used for mesh refinement. The linear conforming finite element space is used to solve the semilinear Neumann problem (6.10) and the finite element space Wk1 is used to solve the dual problem (3.6). The corresponding numerical results are presented in Figs. 10 and 11. The computational time of Algorithm 3 with the accuracy of tolerance 1E-2 is propose in Table 5. From Figs. 10 and 11, we can find the proposed cascadic adaptive finite element method is able to obtain the optimal accuracy. From Table 5, we can also find that the complementary a posteriori error estimator has a better efficiency than the residual type a posteriori error estimator.

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

Errors for semilinear Neumann equation

17

Mesh after 10 iterations

η (u ,Ω) k

1

10

k

η(u ,p1) k

k

slope=−1/2 0

Errors

10

−1

10

−2

10

3

10

4

10

5

10

Number of elements Fig. 10. Left: Errors of Algorithm 3 for Example 4 when the complementary a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Errors for semilinear Neumann equation

Mesh after 10 iterations

ηk(uk,Ω)

1

10

η(u ,p1) k

k

slope=−1/2 0

Errors

10

−1

10

−2

10

3

10

4

10

5

10

Number of elements Fig. 11. Left: Errors of Algorithm 3 for Example 4 when the residual type a posteriori error estimator is adopted in mesh refinement. Right: The corresponding mesh after 10 adaptive refinements.

Table 5 The computational time (in seconds) of the proposed cascadic adaptive algorithm with the same tolerance TOL for Example 4. Algorithm time

Computational time

Error estimate

Algorithm 3 with complementary a posteriori error estimator Algorithm 3 with residual type a posteriori error estimator

57.33 96.45

8.76E−3 9.01E−3

Table 6 The computational time (in seconds) of the proposed cascadic adaptive algorithms, direct AFEM, AFEM based on multilevel correction method for Example 4. Algorithm time

Computational time

Error estimate

Direct AFEM Algorithm 3 with complementary a posteriori error estimator AFEM based on multilevel correction method

274.51 57.33 78.34

8.72E−3 8.76E−3 8.79E−3

18

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

In order to show the efficiency of Algorithm 3, we also propose the computational time of direct AFEM and AFEM based on multilevel correction method. The corresponding computational times are presented in Table 6. From Table 6, we can also find that our algorithm have a better efficiency. 7. Concluding remarks In this paper, a type of asymptotic exact error estimator is proposed for the semilinear Neumann equation. In addition, two types of cascadic adaptive methods for semilinear Neumann equation are proposed based on the complementary error estimate. The cascadic multigrid method is simple and efficient. But it is not an easy job to extend it to adaptive case since we need to know the number of smoothing steps in advance. The current literatures dealing with adaptive cascadic multigrid method for semilinear problem is rarely found. Here, the complementary based error estimate provides a candidate and we make an attempt for studying the adaptive cascadic algorithm. The first scheme is based on the Newton iteration, which needs to do some smoothing steps for a linearized boundary value problem on each adaptive finite element space. The second scheme is based on the multilevel correction method, which contains some smoothing steps on each adaptive finite element space and some solving steps to semilinear Neumann equations on a very coarse space. The number of smoothing steps is determined by the complementary error estimate and the number of degrees of freedom. The efficiency of the presented cascadic algorithms is validated by a variety of examples. In this paper, we only consider a simple semilinear Neumann problem. The algorithms presented in this paper can also be extended to some other more complicated nonlinear problems, such as the models presented in [2–5,15,23], etc. By similar technique, solving these nonlinear models will not be significantly more complicated than the linear boundary value problem. These will be our future work. Acknowledgments This work was supported by National Natural Science Foundation of China (Grant Nos. 11801021, 11571027), Foundation for Fundamental Research of Beijing University of Technology (0 060 0 0546318504), International Research Cooperation Seed Fund of Beijing University of Technology (No. 2018B32). References [1] R.A. Adams, Sobolev spaces, Academic Press, New York, 1975. [2] P. Assari, M. Dehghan, A meshless discrete collocation method for the numerical solution of singular-logarithmic boundary integral equations utilizing radial basis functions, Appl. Math. Comput. 315 (2017) 424–444. [3] P. Assari, M. Dehghan, Solving a class of nonlinear boundary integral equations based on the meshless local discrete Galerkin (MLDG) method, Appl. Numer. Math. 123 (2018) 137–158. [4] P. Assari, M. Dehghan, Application of thin plate splines for solving a class of boundary integral equations arisen from Laplaces equations with nonlinear boundary conditions, Int. J. Comput. Math. 96 (2019) 170–198. [5] P. Assari, M. Dehghan, A meshless Galerkin scheme for the approximate solution of nonlinear logarithmic boundary integral equations utilizing radial basis functions, J. Comput. Appl. Math. 333 (2018) 362–381. [6] I. Babuška, W. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (1978) 736–754. [7] I. Babuška, W. Rheinboldt, A-posteriori error estimates for the finite element method, Int. J. Numer. Methods Eng. 12 (1978) 1597–1615. [8] A. Bermudez, R. Rodriguez, D. Santamarina, A finite element solution of an added mass formulation for coupled fluid-solid vibrations, Numer. Math. 87 (20 0 0) 201–227. [9] F. Bornemann, B. Erdmann, R. Kornhuber, A posteriori error estimates for elliptic problems in two and three space dimensions, SIAM J. Numer. Anal. 33 (3) (1996) 1188–1204. [10] F. Bornemann, P. Deuhard, The cascadic multigrid method for elliptic problems, Numer. Math. 75 (1996) 135–152. [11] J.H. Bramble, X. Zhang, The analysis of multigrid methods, Handb. Numer. Anal. 7 (20 0 0) 173–415. [12] S. Brenner, L. Scott, The mathematical theory of finite element methods, Springer-Verlag, New York, 1994. [13] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, SpringerVerlag, New York, 1991. [14] J. Cascon, C. Kreuzer, R. Nochetto, K. Siebert, Quasi-optimal convergence rate for an adaptive finite element method, SIAM J. Numer. Anal. 46 (5) (2008) 2524–2550. [15] H. Chen, L. He, A. Zhou, Finite element approximations of nonlinear eigenvalue problems in quantum physics, Comput. Methods Appl. Mech. Eng. 200 (21) (2011) 1846–1865. [16] P.G. Ciarlet, in: The finite element method for elliptic problems, North-Holland, New York, 1978. [17] P. Deuflhard, P. Leinen, H. Yserentant, Concepts of an adaptive hierarchical finite element code, IMPACT Compt. Sci. Engrg. 1 (1989) 3–35. [18] P. Deuflhard, M. Weiser, Adaptive numerical solution of PDEs, De Gruyter, Berlin, 2012. [19] W. Dórfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (3) (1996) 1106–1124. [20] W. Hackbusch, Multi-grid methods and applications, Springer-Verlag, Berlin, 1985. [21] X. Han, H. Xie, F. Xu, A cascadic multigrid method for eigenvalue problem, J. Comput. Math. 322 (1) (2017) 747–759. [22] Y. Huang, Z. Shi, T. Tang, W. Xue, A multilevel successive iteration method for nonlinear elliptic problem, Math. Comp. 73 (2004) 525–539. [23] W. Layton, L. Toblska, A two level method with backtracking for the Navier–Stokes equation, SIAM J. Numer. Anal. 35 (5) (1998) 2035–2054. [24] Q. Lin, H. Xie, A multi-level correction scheme for eigenvalue problems, Math. Comp. 84 (291) (2015) 71–88. [25] K. Mekchay, R. Nochetto, Convergence of adaptive finite element methods for general second order linear elliptic PDEs, SIAM J. Numer. Anal. 43 (2005) 1803–1827. [26] P. Morin, R. Nochetto, K. Siebert, Convergence of adaptive finite element methods, SIAM Rev. 44 (4) (2002) 631–658. [27] P. Neittaanmaki, S. Repin, Reliable methods for computer simulation, error control and a posteriori estimates, in: Studies in Mathematics and its Applications, 33, Elsevier Science B.V., Amsterdam, 2004. [28] R. Nochetto, K. Siebert, A. Veeser, Theory of adaptive finite element methods: an introduction, in: R.A. DeVore, A. Kunoth (Eds.), Multiscale, nonlinear and adaptive approximation, Springer, Berlin, 2009. [29] V. Shaidurov, Multigrid methods for finite elements, Springer, 1995. [30] V. Shaidurov, Some estimates of the rate of convergence for the cascadic conjugate-gradient method, Comput. Math. Appl. 31 (1996) 161–171.

F. Xu and Q. Huang / Applied Mathematics and Computation 362 (2019) 124540

19

[31] V. Shaidurov, L. Tobiska, The convergence of the cascadic conjugate-gradient method applied to elliptic problems in domains with re-entrant corners, Math. Comput. 69 (20 0 0) 501–520. [32] R. Stevension, Optimality of a standard adaptive finite element method, Found. Math. Comput. 7 (2) (2007) 245–269. [33] T. Vejchodský, Complementarity based a posteriori error estimates and their properties, Math. Comput. Simulation 82 (2012) 2033–2046. [34] L. Wang, X. Xu, The basic mathematical theory of finite element methods, Science Press, Beijing, 2004. [35] H. Xie, A multigrid method for eigenvalue problem, J. Comput. Phys. 274 (2014) 550–561. [36] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Rev. 34 (4) (1992) 581–613. [37] J. Xu, Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal. 33 (5) (1996) 1759–1777. [38] H. Yu, J. Zeng, A cascadic multigrid method for a kind of semilinear elliptic problem, Numer. Algorithms 58 (2) (2011) 143–162.