Numerical Solutions to Nonsmooth Dirichlet Problems Based on Lumped Mass Finite Element Discretization

and Applied Analysis 3 In the discrete problem (11), we do not have to calculate the nonsmooth term (max{0, ∑N i=1 μ i φ i }, φ j ). Instead, we need only to calculate the following simple term


Introduction
In this paper, we consider a lumped mass finite element method (FEM) to the following Dirichlet problem for a nonsmooth elliptic equation: where Ω ⊂  2 is a bounded convex domain with a Lipschitz continuous boundary Ω,  > 0 is a constant,  is a given smooth function, and  + = max{0, }.The above nonsmooth elliptic problem has many applications.For instance, it can arise from the MHD equilibria, thin stretched membranes problems, or reaction-diffusion problems (see, e.g., [1][2][3][4]).In order to solve problem (1), firstly, it is generally discretized by a finite element (volume) method or a finite difference method, and then various numerical algorithms are constructed to solve the corresponding discrete problems (see, e.g., [1,[4][5][6][7][8] and the references therein).In this paper, we apply a lumped mass finite element to approximate problem (1) by introducing the lumping domain for each node to deal with the nonsmooth term.We refer to [9][10][11][12] for such schemes of lumped mass type.For nonsmooth problem (1), one advantage of the lumped mass FEM is that one can calculate the nonsmooth terms in the finite element equations very easily and numerical quadrature algorithms are no longer needed.Another advantage of this method is that the operator in the finite element problem is diagonally isotone and off-diagonally antitone and thereby some monotone convergent algorithms can be applied (see, e.g., [13,14]).In this paper, we will prove that the FEM error in energy norm is the same as that of the standard FEM.Since the finite element problem is a nonsmooth equation, we will apply nonsmooth Newton-like algorithm to solve it and focus our attention on the monotone convergence property of the algorithm for some special inner iterators.Throughout this paper, we adopt the standard notations for Sobolev spaces  , (Ω) on Ω with norm ‖ ⋅ ‖ ,,Ω and seminorm | ⋅ | ,,Ω .We denote  ,2 (Ω) by   (Ω) and and let  1 0 (Ω) be the subspace of  1 (Ω) with vanishing traces on Ω; that is, We will use the letter  or  (with or without subscripts) to denote a general positive constant independent of the mesh The remainder of the paper is organized as follows.
In Section 2, we present a lumped mass finite element to problem (1) and estimate its error in energy norm (or equivalently  1 -norm).In Section 3, we present nonsmooth Newton-like algorithm and analyze its convergence, especially its monotone convergence.In Section 4, we illustrate some numerical results to confirm the theoretical results we obtained and the efficiency of the algorithm.And finally, in the last section, we make a simple conclusion.

A Lumped Mass Finite Element Approximation to the Nonsmooth Dirichlet Problem
The weak form of ( 1) is to find  ∈  1 0 (Ω) such that where (⋅, ⋅) is the inner product defined by (, V) ≡ ∫ Ω V.
Let T ℎ be a triangulation of Ω ℎ , where the mesh size ℎ denotes the maximum diameter of its triangle elements and let Ω ℎ be a polygonal approximation to Ω with a boundary Ω ℎ .For simplicity, we assume that Ω ℎ = Ω.Let  denote the triangle element of the triangulation T ℎ , and let   ∈ Ω,  = 1, 2, . . .,  ℎ , denote the interior nodes and   ∈ Ω,  =  ℎ + 1, . . .,  ℎ , the boundary nodes, respectively.Let   and   be defined as follows: We assume that the triangulation T ℎ satisfies the following quasi-uniform mesh conditions: Associating with T ℎ , let  ℎ be the piecewise linear finite element space defined by where P 1 is the set of the polynomials of degree ≤1, and for each  = 1, 2, . . .,  ℎ ,   ≥ 0 is the basis function corresponding to the node   .Let  ℎ 0 =  ℎ ⋂  1 0 (Ω).Then  ℎ 0 = span{ 1 ,  2 , . . .,   ℎ }, and the standard piecewise linear FEM approximation of ( 1) is to find  ℎ ∈  ℎ 0 such that Let =1     be the solution of problem (9).Then, problem (9) leads to the following system of nonsmooth equations: find  ∈   ℎ such that where Generally, the operator  ℎ () in problem (10) is not diagonally isotone and off-diagonally antitone.Moreover, from the computational point of view, it is tedious to calculate the nonsmooth term (max{0, ∑  ℎ =1     },   ) directly and some numerical quadrature algorithms have to be used.Instead of solving the finite element problem (9), the lumped mass finite element problem is to find  ℎ ∈  ℎ 0 such that where  ℎ :  0 (Ω) →  ℎ is defined by Here, with   being the characteristic function on the lumping region Ω  at node  ℎ  .For any node  ℎ  ∈ Σ 0 ℎ , its lumping region is the region by joining the centroids of the triangle elements which have  ℎ  as a common vertex to the midpoints of the edges which have  ℎ  as a common extremity [9,10] (see an example in the shaded part in Figure 1).
In the discrete problem (11), we do not have to calculate the nonsmooth term (max{0, ∑  ℎ =1     },   ).Instead, we need only to calculate the following simple term where =1     is the solution of problem (11), meas(Ω) denotes the area of Ω.Furthermore, it will be seen in Section 3 that the operator in problem ( 11) is a diagonally isotone and off-diagonally antitone operator if we assume that the triangulation satisfied the so-called maximum angle condition.
In the sequel, we estimate the error between  and  ℎ .To this end, we introduce two auxiliary problems of finding ũℎ and wℎ ∈  ℎ 0 , such that respectively.Between the solutions of problems ( 15) and ( 16) we have the following results.
According to the finite element theory, the following lemma is obvious [15].We give a simple proof for completion.Lemma 2. Let  and ũℎ be the solutions of problems ( 5) and (15), respectively.Then Proof.We can prove the lemma by following standard steps of estimating the error in energy norm.Setting V = V ℎ in ( 5) and combining with (15), we have and then Therefore, where   is the nodal interpolation of  on  ℎ and we assume it in  2 (Ω), while the "≲" comes from the following well known interpolation results: Therefore, by (26) and Poincaré's inequality, we have Noting ( 11) and ( 16), we have where the first "≲" is from (20), the second "≲" is from (22), and the last "≲" is from (28).Hence, This together with (28) obtains That is, the following theorem holds.

Nonsmooth Newton-Like Method
In this section, we consider some numerical solvers to the discrete problem (11).First, we give a definition [14,16].
Step 3. Compute the step length   , an approximation of the solution ξ of the Newton's equations such that Here, F = − Fℎ (  ),   ∈ co{  Fℎ (  )}; that is,   is a generalized derivative of  ℎ at   , and {  } is a sequence of positive numbers.
Noting that with  min ( ℎ ) > 0 being the smallest eigenvalue of matrix  ℎ , Fℎ is strongly monotone.On the other hand, noting that  ℎ () = ( +  meas(Ω  )) =1 is diagonal, that is, ( ℎ ())  is dependent only on the th entries   [14], Fℎ is a diagonal mapping.Assume furthermore that the triangulation satisfies the maximum angle condition, that is, for any  ∈ T ℎ , the maximum angle of  is less than /2.The matrix  ℎ = (  ) is then a symmetric -matrix.That is to say, the matrix  ℎ , with positive diagonals and nonpositive off-diagonals, has a nonnegative inverse (see, e.g., [18]).By the definition of Fℎ , where   is the th unit vector in   ℎ .Hence, for  1 ≥  2 , where we have used the inequalities   > 0 and   ≤ 0 ( ̸ = ).Equation (39) shows that Fℎ is diagonally isotone and offdiagonally antitone.
In Algorithm 5, we usually choose   → 0, so that the local convergence rate of Algorithm 5 is -superlinear (see, e.g., [19]).Furthermore, According to (39), some monotone convergent algorithms can be constructed.As an example, in the following, we will investigate the monotone convergence of nonsmooth GS-Newton algorithm, in which the subproblem (35) is solved iteratively by Gauss-Seidel iteration.In other words, in nonsmooth GS-Newton algorithm,   is generated by   -times Gauss-Seidel iteration with zero initial: where The following property is important to the monotone convergence.Lemma 6.Let  ℎ be an -matrix and ] satisfy Fℎ (]) ≥ 0.
Similar to Theorem 9, we have the following conclusion according to Remark 7. Proposition 10.Let  ℎ be an -matrix, and let {  } be generated by nonsmooth GS-Newton algorithm with the initial  0 satisfying Fℎ ( 0 ) ≤ 0. Then {  } is monotonically increasing and convergent to , where  is the solution of nonsmooth equation (34).
Remark 11.If we solve (35) by SOR iteration, that is, in the inner iteration, we use the following iteration scheme: Algorithm 5 becomes nonsmooth SOR-Newton algorithm, where  ∈ (0, 2) is a relaxation factor.It is not difficult to verify that the monotone convergence of the algorithm can be derived similar to the proof of Theorem 9.

Some Numerical Examples
In this section, we carry out some numerical experiments to confirm the theoretical results we have obtained.The experiments are performed under Windows XP and MATLAB v7.10 (R2010a) running on a personal computer with an Intel Core 2 Duo CPU at 2.20 GHz and 1.00 GB of memory.
It is shown in Table 1 that the convergent error order in energy norm is (ℎ), which is consistent with Theorem 3.
The monotone convergence of nonsmooth SOR-Newton algorithm can be investigated in the experiments if we choose the initial suitably.But we do not list the details here.In Table 2, we listed CPU times spent in nonsmooth SOR-Newton algorithm for  = 10, 100, and 1000, respectively.In the experiments, we let  = 10 −8 ,  0 = 0,   = min{10 −2 /, ‖ Fℎ (  )‖}, and  = 1.93.As a comparison, we also present corresponding results in Table 3 for the active set method (ASM), which is presented in [5] based on primaldual active set method for constrained optimal control problem (in Table 3, "-" indicates that the corresponding CPU time is more than 6000 seconds).
It follows from Tables 2 and 3 that nonsmooth SOR-Newton algorithm performs much better than ASM.The reason may be that a pair of discrete PDEs have to be solved at each step for ASM.
It can be seen in Tables 2 and 3 that the smaller ℎ is, the more the CPU time becomes.Indeed, by (46), we have for any where the second equality is from ( ℎ   ,  ℎ   ) = 0 for  ̸ =  and the "≲" is from (22) and the inverse inequality, while the "≅" is from the following equivalent property between the  2 -norm and the mesh-dependent norm ‖| ⋅ |‖ ℎ under the condition of quasi-uniform grid (see, e.g., [20]): Therefore, the condition number of   satisfies Cond(  ) ≅ ℎ −2 .This indicates that the conditioner of the linear system (35) deteriorates as the mesh ℎ becomes finer.Roughly speaking, for the large scale problem, much time may be spent on solving linear subproblem (35).This stimulates us to construct numerical solvers combining nonsmooth Newton with multigrid technique.A simple idea is to apply linear multigrid directly to the linear systems produced in Newtonlike method for the correction term at each iteration.We refer to [21,22] and the references therein for this kind of methods.
In Table 4, we present some numerical results for this kind of nonsmooth multigrid-Newton algorithm, in which V-cycle multigrid is used to solve Newton's equations (35) and Gauss-Seidel iteration is taken as the smoother (with one time presmooth and two times postsmooth at each cycle).It follows from the Tables 2-4 that compared with nonsmooth SOR-Newton algorithm and ASM, nonsmooth multigrid-Newton algorithm performs excellently for small ℎ.

Concluding Remarks
In this paper, we apply a lumped mass finite element to approximate Dirichlet problems for nonsmooth elliptic equations.The operator in the relating discrete system of nonsmooth equations is diagonal, which leads to a diagonal generalized gradient matrix for the nonsmooth term.It is proved that the FEM error in energy norm is (ℎ), which is optimal.For the discrete finite element problem, we use nonsmooth Newton-like algorithm to solve it and locally -superlinear convergence of the algorithm can be obtained directly by the existing literatures.Especially, since the operator in the discrete system of nonsmooth equations is diagonally isotone and off-diagonally antitone, the monotone convergence is obtained for the algorithm if Gauss-Seidel iteration (or more generally SOR ) is used to solve the subproblem iteratively.In other words, the iterate generated by the algorithm converges to the solution monotonically from an upper-solution or a lower-solution of the problem.In the numerical experiments, we adopt SOR or multigrid as the inner iterator.Numerical results indicate that the algorithm is efficient.

Table 1 :
The errors in energy norm.

Table 2 :
The CPU times for nonsmooth SOR-Newton algorithm.