Quadratic convergence of monotone iterates for semilinear elliptic obstacle problems

In this paper, we consider the numerical solution for the discretization of semilinear elliptic complementarity problems. A monotone algorithm is established based on the upper and lower solutions of the problem. It is proved that iterates, generated by the algorithm, are a pair of upper and lower solution iterates and converge monotonically from above and below, respectively, to the solution of the problem. Moreover, we investigate the convergence rate for the monotone algorithm and prove quadratic convergence of the algorithm. The monotone and quadratic convergence results are also extended to the discrete problems of the two-sided obstacle problems with a semilinear elliptic operator. We also present some simple numerical experiments.

To solve problem (.), we generally apply a finite difference or finite element approximation to obtain a discrete problem. If we use standard finite difference or lumping finite element approximation with a Delaunay triangle triangulation (see, e.g., in [-]), the discrete problem can be a finite dimensional nonlinear complementarity problem of finding u ∈ K = {v ∈ R n : v ≥ φ} such that where A = (a ij ) ∈ R n×n is an M-matrix, f : R n → R n is a diagonal and nondecreasing monotone mapping. In other words, matrix A has a nonnegative inverse A - as well as nonpositive off-diagonals, and mapping f has a form of In this paper, we also assume that A is a positive definite matrix. That is, there exists a positive constant σ such that, for any vector v ∈ R n , These conditions can be satisfied for suitable small meshsize, see, e.g., in [, , ]. The numerical algorithms for solving problem (.) have been developing rapidly. Among these algorithms, there are a kind of monotone algorithms based on the upper or lower solutions of the problem. We refer to [, , , ] and the references therein for details. These algorithms can be regarded as extensions of the monotone algorithms for solving elliptic boundary value problems or their discretizations (see, e.g., in [, -]). In these algorithms, any generated iterate is an upper (or lower) solution sequence which then converges to the solution monotonically. In this paper, we extend the monotone iterative approach for elliptic boundary value equation, presented in [], to complementarity as well as two-sided obstacle problems. By using a pair of upper and lower solutions as two initial iterates, one can construct two monotone sequences which converge monotonically from above and below, respectively, to the solutions of the problems. Especially, the initial iterative solutions in the monotone iterative algorithms can be obtained directly by solving two discrete linear complementarity problems without any knowledge of the exact solution. Quadratic convergence is also proved for the algorithms.
The structure of the paper is as follows. In Section , we provide two procedures, in which only a pair of linear complementarity problems are needed to be solved. Following these procedures, we can obtain a pair of upper and lower solutions of the nonlinear complementarity problem we are concerned with. In Sections  and , we propose a monotone iterative algorithm and deal with the quadratic convergence of the monotone iterates, respectively. In Section , we extend the results obtained in Sections - to the two-sided obstacle problem. In Section , we present some simple numerical experiments.

Upper and lower solutions and their initializations
The approach presented in this paper is based on the upper and lower solutions of the problems. In this section, we introduce the definitions of upper and lower solutions of problem (.) and discuss their properties.
Nonlinear complementarity problem (.) is equivalent to the following system of nonsmooth equations: According to (.), let Then S  and S - are called the upper and lower solution sets of problem (.), respectively. Besides, any element of S  (S - ) is called an upper (lower) solution of problem (.). It is well known that the solution of the problem is the minimal (maximal) element of upper (lower) solution set (see, e.g., in [, , ]). It is obvious that v ∈ S  is equivalent to That is,S - consists of lower solutions in K . Since the solution of problem (.) is in K , in the sequel, we will consider the lower solution inS - . Obviously, φ is a candidate of such lower solutions.
Lemma . For any w ∈ R n , letū be a solution of the following linear complementarity problem (LCP) of findingū ∈ K such that Proof Assuming thatū satisfies (.) and |Aw + f (w)| + Aw + f (w) ≥ , where z =ūw. It follows from (.) again that Since A is an M-matrix, so is A + * (see, e.g., in [, ]). That is, (A + * ) - ≥ . This together with (.) implies z ≥ , and hence by (.), For the lower solution, we have a similar result as follows.
Lemma . For any w ∈ R n , let u be a solution of the following LCP of finding u ∈ K such that Then u ∈S - .
Problems (.) and (.) can be regarded as lower and upper obstacle problems with variablesū and -u, respectively. According to Lemmas . and ., we can obtain a pair of upper and lower solutions of nonlinear complementarity problem (.) by solving two linear complementarity problems (.) and (.). As for linear complementarity problems, there are many classic and efficient iterative or direct algorithms to solve them. We refer to [, -] for further discussions.

Monotone iterative algorithm for complementarity problem
In this section, we propose an algorithm for solving the nonlinear complementarity problem (.) and discuss its monotone convergence. Firstly, we present the algorithm as follows.
For k ≥ , we calculate a pair of iterates by solving the following linear complementarity problems of finding u (k) α ∈ K , α = , -, respectively, such that Direct algorithms with a polynomial computational complexity are fruitful for subproblems (.). Especially, there are a lot of polynomial algorithms for linear complementarity problems with an M-matrix. We refer to [-] for more details.
By Lemmas . and ., for any w ∈ R n ,ū and u, generated by (.) and (.), are in S  andS - , respectively. So, we may let the initials in Algorithm . be obtained by solving (.) and (.). That is, is also an M-matrix. Therefore, similar to S  andS - , we can define the upper and lower solution sets of problem (.), respectively, as follows: It is sufficient to prove that (z (k) ) + = . By (.) and (z (k) ) + ≥ , we have Thereby, This together with (.) gives By the use of (.), we have then Hence, This together with (.) implies that the right-hand side of the equality in (.) is nonpositive, and then where the last equality is from ((z (k) ) + ) i ((z (k) ) -) i = , and the last inequality is from ((z (k) ) + ) i ((z (k) ) -) j ≤  and a ij ≤  for i = j. By the positive definite assumption of matrix A, we conclude (z (k) ) + = , and the proof is then complete.
The following theorem gives the monotone convergence of Algorithm ..

Quadratic convergence rate of the monotone algorithm
Introduce the notations Here, we have assumed that f is twice continuously differentiable. The following theorem gives the quadratic convergence of Algorithm ..
Theorem . Let the sequences {u (k) α }, α = , -, be generated by Algorithm .. Then the following estimate holds: Proof The linear complementarity problems (.) can be reformulated as the following variational inequality problems of finding u (k) α ∈ K , α = , -, respectively, such that where (v, w) = v T w. Letting v  = u (k) - and v - = u (k)  in above variational forms, we get Denoting u (k) -u (k)  by z (k) and adding above two inequalities, we obtain Then, by (.), we have That is, Taking into account (.), we conclude that Therefore, by (.), This together with (.), (.) as well as (.) and (.) implies That is, (.) holds. The proof is then completed.
By estimate (.), the quadratic convergence holds true for the difference u (k) u (k) - between the upper and lower iterative solutions. This is not like the usual quadratic convergence estimate. Nevertheless, similar to the proof of Theorem ., we can get and And then, Noting that u (k) - ≤ u ≤ u (k)  , we can obtain the following standard estimates in a similar way.

Theorem . Let the sequences {u (k)
α }, α = , -, be generated by (.). Then the following estimates hold: Estimates (.) and (.) indicate that when f i (v i ) (i = , , . . . , n) have convex (or concavity ) property in a neighborhood of the solution of problem (.), the maximal (or minimal) sequence, generated by (.), converges quadratically to the solution of the problem.

Extensions to a two-sided obstacle problem
In this section, we extend the results obtained in the previous sections to the case of a two-sided obstacle problem. Consider the discrete two-sided obstacle problem of finding u ∈ K such that where A and f are defined as before, Problem (.) is equivalent to the following system of nonsmooth equations: max min Au + f (u), uφ , uψ =  (.) and the following variational inequality problem of finding u ∈ K such that According to (.), we define the upper and lower solution sets of problem (.), respectively, as follows: which is equivalent to and which is equivalent to Similar to the case of complementarity problem, the solution of problem (.) is the minimal (maximal) element of S  (S - ).
Obviously, φ and ψ are candidates of S  and S - , respectively. In the following, we present two schemes, which can produce an upper or a lower solution of problem (.) from any w ∈ R n . Scheme . Let w ∈ R n .
Step . Solve the following LCP of findingū ≥ φ such that where * andḡ are the same as those given in Section .
Step . Let where e ∈ R n is the vector of ones.
By Scheme .,ũ ≤ū andũ ∈ K . Moreover, ifũ i < ψ i , where the first inequality is from Lemma ., τ i =  andū i =ũ i , the second inequality is from the definition of M-matrix and τ j ≤  for each j. Thereby, we obtain the following result.

Lemma . Letũ be produced by Scheme
Similarly, we can obtain a lower solution of the problem by the following scheme.
Step . Solve the following LCP of finding u ≤ ψ such that where * and g are the same as those given in Section .
Step . Let Similar to Lemma ., we have the following result.
By Schemes . and ., we can obtain a pair of upper and lower solutions of two-sided obstacle problem (.) by solving two affine upper and lower obstacle problems (.) and (.), instead of solving one two-sided obstacle problem. To our knowledge, as for the twosided obstacle problem, direct algorithms with polynomial computational complexity are few.
Algorithm . Let the initials u ()  ∈ S  and u () - ∈ S - . For k ≥ , we calculate a pair of iterates by solving the following affine two-sided obstacle problems of finding u (k) α ∈ K , α = , -, respectively, such that where matrix A (k) and mapping r are the same as those in Algorithm .. and Problem  We discuss the following nonlinear complementarity problem: Here F(u) = Au + D(u) + f , with A being the same as in Problem . The vector f is generated from a uniform distribution in the interval (-, ). Let D(u) = (D i ) : R n → R n be a given diagonal mapping with D i : R → R, for i = , , . . . , n. The components of D(x) are D i (u) = arctan(u i ).
We compare different algorithms from the point of view of iteration numbers and cpu times (seconds). Here, we consider three algorithms: Algorithm ., denoted by (AL); the semismooth equation approach proposed in [], denoted by SSN; and primal-dual algorithm proposed in [], denoted by PDA.
In the algorithm AL, we choose initial point u () - =  and u ()  to be obtained by Lemma . for all problems. All subproblems are solved by PSOR and the tolerance in PSOR is chosen to be equal to  - in ·  norm. In order to determine the relaxation parameter ω, we consider the following nonlinear complementarity problem: Here F(u) = Au + D(u) + f , where A is the same in Problems  and , D(u) = (D i ) : R n → R n being a given diagonal mapping like in Problem . Set D i (u i ) = u  i , and let f i = -. Fix h =   , and use AL to solve the above problem. We change the ω, and present the result in Table . From the table, we can see that the relaxation parameter ω = . is a good choice, and we use this in all problems.
The termination criterion of the algorithm AL is chosen to be u (k) u (k) - ≤  - . In the algorithm SSN, we choose initial point u  = , tolerance =  - , p = , ρ = ., β = . and H k ∈ ∂ B is defined by the procedure proposed in []. In the algorithm PDA, we fix c =  and choose the initial point u  = . The stopping criterion is that the active sets  are not changed between two iterations. In this algorithm, the subproblems are systems of nonlinear equations, which are solved by Newton iteration. The numerical results are listed in Tables -, where 'iter' denotes the iteration number for the algorithm converging to the solution and 'cpu' denotes the execution time.
From Tables -, we can easily see that the iteration numbers of Algorithm . are stable, which may also mean that the initials obtained by (.) may be a good solution guess. While for SSN and PDA, the iteration numbers increase when the dimensions of Problem  become large. As we can see from Table , the algorithm we proposed seems to be more effective for solving large-scale problems. The main reason may be as follows. Algorithm . takes only a few iterations for all problems, and in each iteration it only needs to solve two linear complementarity problems. These subproblems are solved by PSOR rapidly. For SSN, it takes a lot of time to solve the system of linear equations to obtain directions, especially for large-scale problems. For PDA, by using the active set strategy, in each iteration, it only needs to solve a reduced system of linear equations, where the dimension is much less than the original one. Therefore, the execution time for PDA is also less compared with SSN.

Conclusions
In this paper, we have considered the numerical solution for the discretization of semilinear elliptic complementarity problems. Based on the upper and lower solutions of the problem, we have proposed a monotone algorithm and proved that iterates are a pair of upper and lower solution iterates and converge monotonically from above and below, respectively, to the solution of the problem. Moreover, we have established quadratic convergence of the algorithm. The limited numerical results showed that the proposed algorithm is effective.