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

Authors: Haixiong Yu; Jinping Zeng
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We apply a lumped mass finite element to approximate Dirichlet problems for  nonsmooth elliptic equations. It is proved that the lumped mass FEM approximation error in energy norm is the same as that of standard piecewise linear finite element approximation. Under the  quasi-uniform mesh condition and the maximum angle condition, we show that the operator in the finite element problem is diagonally isotone and off-diagonally antitone. Therefore, some monotone convergent algorithms can be used. As an example, we prove that the nonsmooth Newton-like algorithm is convergent monotonically
if Gauss-Seidel iteration is used to solve the Newton's equations iteratively. Some numerical experiments are presented.


---

## Body

## 1. Introduction

In this paper, we consider a lumped mass finite element method (FEM) to the following Dirichlet problem for a nonsmooth elliptic equation:(1)
-
Δ
u
+
λ
u
+
=
f
in
Ω
,
u
=
0
on
⁡
∂
Ω
,
where Ω
⊂
R
2 is a bounded convex domain with a Lipschitz continuous boundary ∂
Ω, λ
>
0 is a constant, f is a given smooth function, and u
+
=
max
⁡
{
0
,
u
}. 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–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–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–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 spacesW
m
,
p
(
Ω
) on Ω with norm ∥
·
∥
m
,
p
,
Ω and seminorm |
·
|
m
,
p
,
Ω. We denote W
m
,
2
(
Ω
) by H
m
(
Ω
) and
(2)
∥
·
∥
m
,
2
,
Ω
=
∥
·
∥
m
,
Ω
,
|
·
|
m
,
2
,
Ω
=
|
·
|
m
,
Ω
,
and let H
0
1
(
Ω
) be the subspace of H
1
(
Ω
) with vanishing traces on ∂
Ω; that is,
(3)
H
0
1
(
Ω
)
≡
{
v
∈
H
1
(
Ω
)
:
v
|
∂
Ω
=
0
}
.We will use the letterc or C (with or without subscripts) to denote a general positive constant independent of the mesh size. When it is not important to keep track of these constants, we will conceal the letter c or C into Xu’s notations “≲” and “≅” for convenience. Here
(4)
a
≲
b
means
a
≤
C
b
,
a
≅
b
means
c
b
≤
a
≤
C
b
.The remainder of the paper is organized as follows. In Section2, we present a lumped mass finite element to problem (1) and estimate its error in energy norm (or equivalently H
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.


## 2. A Lumped Mass Finite Element Approximation to the Nonsmooth Dirichlet Problem

The weak form of (1) is to find u
∈
H
0
1
(
Ω
) such that
(5)
(
∇
u
,
∇
v
)
+
λ
(
u
+
,
v
)
=
(
f
,
v
)
,
∀
v
∈
H
0
1
(
Ω
)
,
where (
·
,
·
) is the inner product defined by (
u
,
v
)
≡
∫
Ω
‍
u
v.Let𝒯
h be a triangulation of Ω
h, where the mesh size h denotes the maximum diameter of its triangle elements and let Ω
h be a polygonal approximation to Ω with a boundary ∂
Ω
h. For simplicity, we assume that Ω
h
=
Ω. Let τ denote the triangle element of the triangulation 𝒯
h, and let X
i
∈
Ω, i
=
1,2
,
…
,
N
h, denote the interior nodes and X
i
∈
∂
Ω, i
=
N
h
+
1
,
…
,
M
h, the boundary nodes, respectively. Let ρ
τ and ρ
¯
τ be defined as follows:
(6)
ρ
τ
=
inf
⁡
{
diam
⁡
(
B
)
:
B
is
a
ball
that
contains
τ
}
,
ρ
¯
τ
=
sup
⁡
{
diam
⁡
(
B
)
:
B
is
a
ball
contained
in
τ
}
.
We assume that the triangulation 𝒯
h satisfies the following quasi-uniform mesh conditions:
(7)
h
≲
ρ
τ
≲
ρ
¯
τ
,
∀
τ
∈
𝒯
h
.
Associating with 𝒯
h, let V
h be the piecewise linear finite element space defined by
(8)
V
h
=
{
v
∈
C
(
Ω
¯
)
:
v
|
τ
∈
𝒫
1
|
τ
,
∀
τ
∈
𝒯
h
}
=
span
⁡
{
ϕ
1
,
ϕ
2
,
…
,
ϕ
M
h
}
,
where 𝒫
1 is the set of the polynomials of degree ≤1, and for each i
=
1,2
,
…
,
M
h, ϕ
i
≥
0 is the basis function corresponding to the node X
i. Let V
0
h
=
V
h
⋂
H
0
1
(
Ω
). Then V
0
h
=
span
⁡
{
ϕ
1
,
ϕ
2
,
…
,
ϕ
N
h
}, and the standard piecewise linear FEM approximation of (1) is to find u
h
∈
V
0
h such that
(9)
(
∇
u
h
,
∇
v
h
)
+
λ
(
u
h
+
,
v
h
)
=
(
f
,
v
h
)
,
∀
v
h
∈
V
0
h
.Letu
h
=
∑
i
=
1
N
h
μ
i
ϕ
i be the solution of problem (9). Then, problem (9) leads to the following system of nonsmooth equations: find μ
∈
R
N
h such that
(10)
F
h
(
μ
)
=
A
h
μ
+
B
h
(
μ
)
-
Q
h
f
=
0
,
where A
h
=
(
(
∇
ϕ
i
,
∇
ϕ
j
)
)
N
h
×
N
h is a symmetric positive definite (SPD) matrix, B
h
(
μ
)
=
λ
(
(
max
⁡
{
0
,
∑
i
=
1
N
h
μ
i
ϕ
i
}
,
ϕ
j
)
)
j
=
1
N
h, and Q
h
f
=
(
f
,
ϕ
j
)
j
=
1
N
h.Generally, the operatorF
h
(
μ
) 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
,
∑
i
=
1
N
h
μ
i
ϕ
i
}
,
ϕ
j
) 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 w
h
∈
V
0
h such that
(11)
(
∇
w
h
,
∇
v
h
)
+
λ
(
(
R
h
w
h
)
+
,
R
h
v
h
)
=
(
f
,
v
h
)
,
∀
v
h
∈
V
0
h
,
where R
h
:
C
0
(
Ω
)
→
K
h is defined by
(12)
R
h
v
=
∑
i
=
1
N
h
‍
v
(
X
i
)
χ
i
.
Here,
(13)
K
h
=
{
μ
h
:
μ
h
=
∑
i
=
1
N
h
‍
μ
i
χ
i
,
∀
μ
i
∈
R
,
i
=
1,2
,
…
,
N
h
}
with χ
i being the characteristic function on the lumping region Ω
i at node P
i
h. For any node P
i
h
∈
∑
h
0, its lumping region is the region by joining the centroids of the triangle elements which have P
i
h as a common vertex to the midpoints of the edges which have P
i
h as a common extremity [9, 10] (see an example in the shaded part in Figure 1).Figure 1
The lumped mass regionΩ
i
h at P
i
h.In the discrete problem (11), we do not have to calculate the nonsmooth term (
max
⁡
{
0
,
∑
i
=
1
N
h
μ
i
ϕ
i
}
,
ϕ
j
). Instead, we need only to calculate the following simple term
(14)
(
(
R
h
∑
i
=
1
N
h
‍
ω
i
ϕ
i
)
+
,
R
h
ϕ
j
)
=
∫
Ω
(
R
h
∑
i
=
1
N
h
‍
ω
i
ϕ
i
)
+
R
h
ϕ
j
=
∫
Ω
j
‍
ω
j
+
=
ω
j
+
meas
⁡
(
Ω
j
)
,
where w
h
=
∑
i
=
1
N
h
ω
i
ϕ
i 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 betweenu and w
h. To this end, we introduce two auxiliary problems of finding u
~
h and w
~
h
∈
V
0
h, such that
(15)
(
∇
u
~
h
,
∇
v
h
)
+
λ
(
u
+
,
v
h
)
=
(
f
,
v
h
)
,
∀
v
h
∈
V
0
h
,
(16)
(
∇
w
~
h
,
∇
v
h
)
+
λ
(
(
R
h
u
)
+
,
R
h
v
h
)
=
(
f
,
v
h
)
,
∀
v
h
∈
V
0
h
,
respectively. Between the solutions of problems (15) and (16) we have the following results.Lemma 1.
Letu
~
h and w
~
h be the solutions of problems (15) and (16), respectively. Then
(17)
∥
w
~
h
-
u
~
h
∥
1
,
Ω
≲
λ
h
.Proof.
By (15) and (16) we have
(18)
0
=
(
∇
(
w
~
h
-
u
~
h
)
,
∇
(
w
~
h
-
u
~
h
)
)
+
λ
(
(
R
h
u
)
+
,
R
h
(
w
~
h
-
u
~
h
)
)
-
λ
(
u
+
,
w
~
h
-
u
~
h
)
.
Then,
(19)
|
w
~
h
-
u
~
h
|
1
,
Ω
2
≤
λ
|
(
(
R
h
u
)
+
-
u
+
,
R
h
(
w
~
h
-
u
~
h
)
)
|
+
λ
|
(
u
+
,
R
h
(
w
~
h
-
u
~
h
)
-
(
w
~
h
-
u
~
h
)
)
|
≤
λ
∥
(
R
h
u
)
+
-
u
+
∥
0
,
Ω
∥
R
h
(
w
~
h
-
u
~
h
)
∥
0
,
Ω
+
λ
∥
u
+
∥
0
,
Ω
∥
R
h
(
w
~
h
-
u
~
h
)
-
(
w
~
h
-
u
~
h
)
∥
0
,
Ω
.
For any real constantsa and b, we have
(20)
|
a
+
-
b
+
|
≤
|
a
-
b
|
,
and
especially
,
|
a
+
|
≤
|
a
|
.
This combining with (19) implies
(21)
|
w
~
h
-
u
~
h
|
1
,
Ω
2
≤
λ
∥
R
h
u
-
u
∥
0
,
Ω
∥
R
h
(
w
~
h
-
u
~
h
)
∥
0
,
Ω
+
λ
∥
u
∥
0
,
Ω
∥
R
h
(
w
~
h
-
u
~
h
)
-
(
w
~
h
-
u
~
h
)
∥
0
,
Ω
≲
λ
h
∥
w
~
h
-
u
~
h
∥
1
,
Ω
.
Here, we used the following properties of R
h:
(22)
∥
R
h
v
-
v
∥
0
,
Ω
≲
h
|
v
|
1
,
Ω
,
∥
R
h
v
∥
0
,
Ω
≲
∥
v
∥
0
,
Ω
+
h
|
v
|
1
,
Ω
≲
∥
v
∥
1
,
Ω
.
Then, we immediately obtain (17) by Poincaré’s inequality.According to the finite element theory, the following lemma is obvious [15]. We give a simple proof for completion.Lemma 2.
Letu and u
~
h be the solutions of problems (5) and (15), respectively. Then
(23)
∥
u
-
u
~
h
∥
1
,
Ω
≲
h
.Proof.
We can prove the lemma by following standard steps of estimating the error in energy norm. Settingv
=
v
h in (5) and combining with (15), we have
(24)
(
∇
u
,
∇
v
h
)
=
-
λ
(
u
+
,
v
h
)
+
(
f
,
v
h
)
=
(
∇
u
~
h
,
∇
v
h
)
,
and then
(25)
(
∇
(
u
-
u
~
h
)
,
∇
v
h
)
=
0
,
∀
v
h
∈
V
0
h
.
Therefore,
(26)
|
u
-
u
~
h
|
1
,
Ω
2
=
(
∇
(
u
-
u
~
h
)
,
∇
(
u
-
u
~
h
)
)
=
(
∇
(
u
-
u
~
h
)
,
∇
(
u
-
u
I
)
)
+
(
∇
(
u
-
u
~
h
)
,
∇
(
u
I
-
u
~
h
)
)
=
(
∇
(
u
-
u
~
h
)
,
∇
(
u
-
u
I
)
)
≤
|
u
-
u
~
h
|
1
,
Ω
|
u
-
u
I
|
1
,
Ω
≲
h
|
u
-
u
~
h
|
1
,
Ω
,
where u
I is the nodal interpolation of u on V
h and we assume it in H
2
(
Ω
), while the “≲” comes from the following well known interpolation results:
(27)
∥
u
-
u
I
∥
1
,
Ω
≲
h
|
u
|
2
,
Ω
.
Therefore, by (26) and Poincaré’s inequality, we have ∥
u
-
u
~
h
∥
1
,
Ω
≲
|
u
-
u
~
h
|
1
,
Ω
≲
h, which implies (23).By Lemmas1 and 2, we have
(28)
∥
u
-
w
~
h
∥
1
,
Ω
≤
∥
u
-
u
~
h
∥
1
,
Ω
+
∥
u
~
h
-
w
~
h
∥
1
,
Ω
≲
max
⁡
{
1
,
λ
}
h
.Noting (11) and (16), we have
(29)
0
=
(
∇
(
w
~
h
-
w
h
)
,
∇
(
w
~
h
-
w
h
)
)
+
λ
(
(
R
h
u
)
+
,
R
h
(
w
~
h
-
w
h
)
)
-
λ
(
(
R
h
w
h
)
+
,
R
h
(
w
~
h
-
w
h
)
)
.
It then follows that
(30)
|
w
~
h
-
w
h
|
1
,
Ω
2
≤
(
∇
(
w
~
h
-
w
h
)
,
∇
(
w
~
h
-
w
h
)
)
+
λ
(
(
R
h
w
~
h
)
+
-
(
R
h
w
h
)
+
,
R
h
w
~
h
-
R
h
w
h
)
=
-
λ
(
(
R
h
u
)
+
-
(
R
h
w
~
h
)
+
,
R
h
(
w
~
h
-
w
h
)
)
≤
λ
∥
(
R
h
u
)
+
-
(
R
h
w
~
h
)
+
∥
0
,
Ω
∥
R
h
(
w
~
h
-
w
h
)
∥
0
,
Ω
≲
λ
∥
R
h
u
-
R
h
w
~
h
∥
0
,
Ω
∥
R
h
(
w
~
h
-
w
h
)
∥
0
,
Ω
≲
λ
∥
u
-
w
~
h
∥
1
,
Ω
∥
w
~
h
-
w
h
∥
1
,
Ω
≲
max
⁡
{
λ
,
λ
2
}
h
∥
w
~
h
-
w
h
∥
1
,
Ω
,
where the first “≲” is from (20), the second “≲” is from (22), and the last “≲” is from (28). Hence,
(31)
∥
w
~
h
-
w
h
∥
1
,
Ω
≲
max
⁡
{
λ
,
λ
2
}
h
.
This together with (28) obtains
(32)
∥
u
-
w
h
∥
1
,
Ω
≤
∥
u
-
w
~
h
∥
1
,
Ω
+
∥
w
~
h
-
w
h
∥
1
,
Ω
≲
max
⁡
{
1
,
λ
,
λ
2
}
h
.That is, the following theorem holds.Theorem 3.
Letu and w
h be the solutions of problems (5) and (11), respectively. Then we have the estimate ∥
u
-
w
h
∥
1
,
Ω
≲
h.


## 3. Nonsmooth Newton-Like Method

In this section, we consider some numerical solvers to the discrete problem (11). First, we give a definition [14, 16].Definition 4.
A mappingF
:
R
n
→
R
n is said to be B-differential at a point x if there exists a function B
F
(
x
)
:
R
n
→
R
n called the B-derivative of F at x, which is positively homogeneous of degree 1 (i.e., B
F
(
x
)
(
t
h
)
=
t
B
F
(
x
)
h for all h
∈
R
n and all t
≥
0), such that
(33)
lim
⁡
∥
F
(
x
+
h
)
-
F
(
x
)
-
B
F
(
x
)
h
∥
t
=
0
.Letw
h
=
∑
i
=
1
N
h
ω
i
ϕ
i be the solution of problem (11). Then, problem (11) leads to the following system of nonsmooth equations: find ω
∈
R
N
h such that
(34)
F
~
h
(
ω
)
=
A
h
ω
+
B
h
(
ω
)
-
Q
h
f
=
0
,
where A
h
=
(
(
∇
ϕ
i
,
∇
ϕ
j
)
)
N
h
×
N
h is an SPD (Symmetric Positive Definite) matrix, B
h
(
ω
)
=
λ
(
ω
j
+
meas
⁡
(
Ω
j
)
)
j
=
1
N
h, and Q
h
f
=
(
f
,
ϕ
j
)
j
=
1
N
h. Due to the occurrence of ω
j
+, F
~
h is not differentiable but is semismooth. So, we use the following nonsmooth Newton-like method to solve problem (34) (see, e.g., [16, 17]).Algorithm 5 (nonsmooth Newton-like method).
Consider the following.Step 1.
Given an initial guessω
0 and a precision ϵ
>
0, let k
≔
0.Step 2.
If∥
F
~
h
(
ω
k
)
∥
∞
≤
ϵ, stop. Otherwise, turn to Step 3.Step 3.
Compute the step lengthξ
k, an approximation of the solution ξ
~
k of the Newton’s equations
(35)
G
k
ξ
~
k
=
F
~
k
,
such that
(36)
G
k
ξ
k
-
F
~
k
=
γ
k
,
with
∥
γ
k
∥
≤
η
k
∥
F
~
k
∥
.
Here, F
~
k
=
-
F
~
h
(
ω
k
), G
k
∈
co
⁡
{
∂
B
F
~
h
(
ω
k
)
}; that is, G
k is a generalized derivative of F
h at ω
k, and {
η
k
} is a sequence of positive numbers.Step 4.
Letω
k
+
1
=
ω
k
+
ξ
k,  k
≔
k
+
1, and turn to Step 2.Noting that(37)
(
F
~
h
(
μ
)
-
F
~
h
(
ν
)
)
T
(
μ
-
ν
)
=
(
μ
-
ν
)
T
A
h
(
μ
-
ν
)
+
λ
∑
i
=
1
N
h
‍
(
μ
j
+
-
ν
j
+
)
meas
⁡
(
Ω
j
)
(
μ
j
-
ν
j
)
≥
(
μ
-
ν
)
T
A
h
(
μ
-
ν
)
≥
λ
min
⁡
(
A
h
)
∥
μ
-
ν
∥
2
,
∀
μ
,
ν
∈
R
N
h
,
with λ
min
⁡
(
A
h
)
>
0 being the smallest eigenvalue of matrix A
h, F
~
h is strongly monotone. On the other hand, noting that B
h
(
ω
)
=
λ
(
ω
j
+
meas
⁡
(
Ω
j
)
)
j
=
1
N
h is diagonal, that is, (
B
h
(
ω
)
)
j is dependent only on the jth entries ω
j [14], F
~
h is a diagonal mapping. Assume furthermore that the triangulation satisfies the maximum angle condition, that is, for any τ
∈
𝒯
h, the maximum angle of τ is less than π
/
2. The matrix A
h
=
(
a
i
j
) is then a symmetric M-matrix. That is to say, the matrix A
h, with positive diagonals and nonpositive off-diagonals, has a nonnegative inverse (see, e.g., [18]). By the definition of F
~
h,
(38)
F
~
h
(
ν
+
t
1
e
i
)
-
F
~
h
(
ν
+
t
2
e
i
)
=
(
t
1
-
t
2
)
A
h
e
i
+
λ
(
(
ν
+
t
1
e
i
)
j
+
-
(
ν
+
t
2
e
i
)
j
+
)
meas
⁡
(
Ω
j
)
,
where e
i is the ith unit vector in R
N
h. Hence, for t
1
≥
t
2,
(39)
[
F
~
h
(
ν
+
t
1
e
i
)
-
F
~
h
(
ν
+
t
2
e
i
)
]
i
=
(
t
1
-
t
2
)
a
i
i
+
λ
[
(
ν
i
+
t
1
)
+
-
(
ν
i
+
t
2
)
+
]
meas
⁡
(
Ω
i
)
≥
0
,
[
F
~
h
(
ν
+
t
1
e
i
)
-
F
~
h
(
ν
+
t
2
e
i
)
]
j
=
(
t
1
-
t
2
)
a
j
i
≤
0
,
j
≠
i
,
where we have used the inequalities a
i
i
>
0 and a
i
j
≤
0 (j
≠
i). Equation (39) shows that F
~
h is diagonally isotone and off-diagonally antitone.In Algorithm5, we usually choose η
k
→
0, so that the local convergence rate of Algorithm 5 is Q-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, ξ
k is generated by t
k-times Gauss-Seidel iteration with zero initial:
(40)
ξ
k
,
0
=
0
,
(
D
k
-
L
k
)
ξ
k
,
l
+
1
=
U
k
ξ
k
,
l
+
F
~
k
,
l
=
0,1
,
…
,
t
k
-
1
,
ξ
k
=
ξ
k
,
t
k
,
where G
k
=
-
L
k
+
D
k
-
U
k.The following property is important to the monotone convergence.Lemma 6.
LetA
h be an M-matrix and ν satisfy F
~
h
(
ν
)
≥
0. Then, ν
≥
ω, where ω is the solution of nonsmooth equation F
~
h
(
ω
)
=
0. Similarly, let ξ satisfies G
k
ξ
-
F
~
k
≥
0. Then, ξ
≥
ξ
~
k, where ξ
~
k is the solution of (35).Proof.
Let(41)
θ
j
=
{
1
,
ν
j
≥
0
,
ω
j
≥
0
,
ν
j
ν
j
-
ω
j
,
ν
j
≥
0
>
ω
j
,
ω
j
ω
j
-
ν
j
,
ν
j
<
0
≤
ω
j
,
0
,
ν
j
<
0
,
ω
j
<
0
.
Then 0
≤
θ
j
≤
1, ν
j
+
-
ω
j
+
=
θ
j
(
ν
j
-
ω
j
), and B
h
(
ν
)
-
B
h
(
ω
)
=
D
ν
,
ω
(
ν
-
ω
), where D
ν
,
ω
=
λ
diag
⁡
(
θ
j
meas
⁡
(
Ω
j
)
). Therefore,
(42)
0
≤
F
~
h
(
ν
)
-
F
~
h
(
ω
)
=
[
A
h
+
D
ν
,
ω
]
(
ν
-
ω
)
.
Since that A
h is an M-matrix and D
ν
,
ω
≥
0 is a diagonal matrix, [
A
h
+
D
ν
,
ω
]
-
1
≥
0 (see, e.g., [14]). Multiplying (42) on the left by [
A
h
+
D
ν
,
ω
]
-
1, we obtain 0
≤
ν
-
ω. The rest part of the lemma can be proved in a similar way.Remark 7.
By Lemma6, for any ν satisfying F
~
h
(
ν
)
≥
0, we have ν
≥
ω. Therefore, ν is usually called an upper-solution of problem F
~
h
(
ω
)
=
0. Similarly, for any ν satisfying F
~
h
(
ν
)
≤
0, we can conclude ν
≤
ω and thereby ν is called a lower-solution of problem F
~
h
(
ω
)
=
0.Lemma 8.
LetA
h be an M-matrix and F
~
h
(
ω
k
)
≥
0. Then
(43)
G
k
ξ
k
,
l
-
F
~
k
≥
0
,
F
~
h
(
ω
k
+
ξ
k
,
l
)
≥
0
,
ξ
k
,
l
+
1
≤
ξ
k
,
l
≤
0
,
where ξ
k
,
l  (l
=
0,1
,
…
,
t
k) are generated by (40).Proof.
It is easy to check that(44)
G
k
ξ
k
,
0
-
F
~
k
=
-
F
~
k
=
F
~
h
(
ω
k
)
≥
0
,
F
~
h
(
w
k
+
ξ
k
,
0
)
=
F
~
h
(
w
k
)
≥
0
.
By the definition of G
k,
(45)
G
k
=
A
h
+
λ
diag
⁡
(
G
j
ω
j
+
∣
ω
j
k
meas
⁡
(
Ω
j
)
)
,
where G
j
ω
j
+
∣
ω
j
k is a generalized derivative of semismooth function ω
j
+ at point ω
j
k. Noting that
(46)
0
≤
G
j
ω
j
+
≤
1
,
G
k is still an M-matrix (see, e.g., [14]). Therefore, (
D
k
-
L
k
)
-
1
≥
0. On the other hand, by (40), we have
(47)
(
D
k
-
L
k
)
ξ
k
,
1
=
U
k
ξ
k
,
0
+
F
~
k
=
-
F
~
h
(
w
k
)
.
Hence
(48)
ξ
k
,
1
=
-
(
D
k
-
L
k
)
-
1
F
~
h
(
w
k
)
≤
0
=
ξ
k
,
0
.
We then conclude that (43) holds for l
=
0. Assume that (43) holds for l. By (40), we have
(49)
(
D
k
-
L
k
)
ξ
k
,
l
+
2
=
U
k
ξ
k
,
l
+
1
+
F
~
k
=
(
D
k
-
L
k
)
ξ
k
,
l
+
1
-
G
k
ξ
k
,
l
+
1
+
F
~
k
,
and then
(50)
ξ
k
,
l
+
2
=
ξ
k
,
l
+
1
-
(
D
k
-
L
k
)
-
1
(
G
k
ξ
k
,
l
+
1
-
F
~
k
)
≤
ξ
k
,
l
+
1
≤
0
,
where we used the inequalities (
D
k
-
L
k
)
-
1
≥
0 and ξ
k
,
l
+
1
≤
ξ
k
,
l
≤
0. By the use of (40) and ξ
k
,
l
+
1
≤
ξ
k
,
l
≤
0 again, we have
(51)
G
k
ξ
k
,
l
+
1
-
F
~
k
=
(
D
k
-
L
k
-
U
k
)
ξ
k
,
l
+
1
-
F
~
k
=
(
D
k
-
L
k
)
ξ
k
,
l
+
1
-
F
~
k
-
U
k
ξ
k
,
l
+
1
=
U
k
(
ξ
k
,
l
-
ξ
k
,
l
+
1
)
≥
0
.
On the other hand, it is easy to calculate that(52)
F
~
h
(
ω
k
+
ξ
k
,
l
+
1
)
=
A
h
(
ω
k
+
ξ
k
,
l
+
1
)
+
λ
(
(
ω
j
k
+
ξ
j
k
,
l
+
1
)
+
meas
⁡
(
Ω
j
)
)
-
Q
h
f
=
G
k
ξ
k
,
l
+
1
-
F
~
k
-
λ
diag
⁡
(
G
j
ω
j
+
∣
ω
j
k
meas
⁡
(
Ω
j
)
)
ξ
k
,
l
+
1
+
F
~
k
+
A
h
(
ω
k
)
+
λ
(
(
ω
j
k
+
ξ
j
k
,
l
+
1
)
+
meas
⁡
(
Ω
j
)
)
-
Q
h
f
=
(
D
k
-
L
k
)
ξ
k
,
l
+
1
-
F
~
k
-
U
k
ξ
k
,
l
+
1
-
λ
diag
⁡
(
G
j
ω
j
+
∣
ω
j
k
meas
⁡
(
Ω
j
)
)
ξ
k
,
l
+
1
-
A
h
ω
k
-
λ
(
(
ω
j
k
)
+
meas
⁡
(
Ω
j
)
)
+
Q
h
f
+
A
h
(
ω
k
)
+
λ
(
(
ω
j
k
+
ξ
j
k
,
l
+
1
)
+
meas
⁡
(
Ω
j
)
)
-
Q
h
f
=
U
k
(
ξ
k
,
l
-
ξ
k
,
l
+
1
)
+
λ
(

-
G
j
ω
j
+
∣
ω
j
k
ξ
j
k
,
l
+
1
]
meas
⁡
(
Ω
j
)
[
(
ω
j
k
+
ξ
j
k
,
l
+
1
)
+
-
(
ω
j
k
)
+
-
G
j
ω
j
+
∣
ω
j
k
ξ
j
k
,
l
+
1
]
meas
⁡
(
Ω
j
)
)
=
U
k
(
ξ
k
,
l
-
ξ
k
,
l
+
1
)
+
λ
(
(
θ
j
k
-
G
j
ω
j
+
∣
ω
j
k
)
ξ
j
k
,
l
+
1
meas
⁡
(
Ω
j
)
)
≥
λ
(
(
θ
j
k
-
G
j
ω
j
+
∣
ω
j
k
)
ξ
j
k
,
l
+
1
meas
⁡
(
Ω
j
)
)
,
where θ
j
k
∈
[
0,1
]  (
j
=
1,2
,
…
,
N
h
) are defined similar to (41). If ω
j
k
>
0, G
j
ω
j
+
∣
ω
j
k
=
1 and then θ
j
k
-
G
j
ω
j
+
∣
ω
j
k
=
θ
j
k
-
1
≤
0. If ω
j
k
≤
0, by ξ
k
,
l
+
1
≤
ξ
k
,
l
≤
0, we have (
ω
j
k
+
ξ
j
k
,
l
+
1
)
+
=
(
ω
j
k
+
ξ
j
k
,
l
)
+
=
0 and then θ
j
k
=
0. Therefore, by (46), we have θ
j
k
-
G
j
ω
j
+
∣
ω
j
k
=
-
G
j
ω
j
+
∣
ω
j
k
≤
0. In one word, θ
j
k
-
G
j
ω
j
+
∣
ω
j
k
≤
0 always holds. This together with (52) and (50) implies that
(53)
F
~
h
(
ω
k
+
ξ
k
,
l
+
1
)
≥
0
.
Combining above inequality with (50) and (51), we complete the proof by the principle of induction.Theorem 9.
LetA
h be an M-matrix, and let {
ω
k
} be generated by nonsmooth GS-Newton algorithm with the initial ω
0 satisfying F
~
h
(
w
0
)
≥
0. Then
(54)
F
~
h
(
ω
k
)
≥
0
,
ω
≤
ω
k
+
1
≤
ω
k
,
k
=
0,1
,
…
,
where ω is the solution of nonsmooth equation (34). Moreover, {
ω
k
} converges to ω monotonically.Proof.
Equation (54) can be obtained directly by Lemmas 6 and 8. Hence, lim
⁡
k
→
∞
ω
k
=
ω
* exists. On the other hand, it follows from (43) that ξ
k
,
l
+
1
-
ξ
k
,
l
≤
0 holds for each k
=
0,1
,
… and l
=
0
,
1
,
…
,
t
k
-
1. Therefore,
(55)
ω
*
⟵
ω
k
+
1
=
ω
k
+
ξ
k
,
t
k
≤
ω
k
+
ξ
k
,
l
+
1
≤
ω
k
+
ξ
k
,
l
≤
ω
k
→
ω
*
,
which implies
(56)
lim
⁡
k
→
∞
ξ
k
,
l
=
0
.
By (45) and (46), we have
(57)
∥
G
k
∥
≤
∥
A
h
∥
+
λ
∥
diag
⁡
(
meas
⁡
(
Ω
j
)
)
∥
and then by (40)
(58)
∥
F
~
h
(
ω
k
)
∥
=
∥
F
~
k
∥
=
∥
G
k
ξ
k
,
l
+
1
+
U
k
(
ξ
k
,
l
+
1
-
ξ
k
,
l
)
∥
≤
(
∥
A
h
∥
+
λ
∥
diag
⁡
(
meas
⁡
(
Ω
j
)
)
∥
)
∥
ξ
k
,
l
+
1
∥
+
∥
U
h
∥
∥
ξ
k
,
l
+
1
-
ξ
k
,
l
∥
→
0
,
k
→
∞
,
where U
h is the strictly upper-triangular part of A
h. Then, we conclude ∥
F
~
h
(
ω
*
)
∥
=
0, which implies ω
*
=
ω. The proof is then completed.Similar to Theorem9, we have the following conclusion according to Remark 7.Proposition 10.
LetA
h be an M-matrix, and let {
ω
k
} be generated by nonsmooth GS-Newton algorithm with the initial ω
0 satisfying F
~
h
(
w
0
)
≤
0. Then {
ω
k
} 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:
(59)
ξ
k
,
0
=
0
,
(
D
k
-
ρ
L
k
)
ξ
k
,
l
+
1
=
[
(
1
-
ρ
)
D
k
+
ρ
U
k
]
ξ
k
,
l
+
ρ
F
~
k
,
l
=
0,1
,
…
,
t
k
-
1
,
ξ
k
=
ξ
k
,
t
k
.
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.


## 4. 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.In our numerical examples, we consider problem (1) with the following data (see [5]): Ω
=
(
0,1
)
×
(
0,1
) and
(60)
f
=
2
[
(
x
-
x
2
)
+
(
y
-
y
2
)
+
π
2
sin
(
2
π
x
)
sin
(
2
π
y
)
]
+

λ
max
⁡
(

0
,
(
x
-
x
2
)
(
y
-
y
2
)

+

1
4
sin
(
2
π
x
)
sin
(
2
π
y
)

)
.It is easy to verify that the exact solution of (1) is u
=
(
x
-
x
2
)
(
y
-
y
2
)
+
(
1
/
4
)
sin
(
2
π
x
)
sin
(
2
π
y
).It is shown in Table1 that the convergent error order in energy norm is O
(
h
), which is consistent with Theorem 3.Table 1
The errors in energy norm.
h
λ
=
10
λ
=
100
λ
=
1000
|
w
h
-
u
|
1
,
Ω
Order
|
w
h
-
u
|
1
,
Ω
Order
|
w
h
-
u
|
1
,
Ω
Order
2
-
3
4.347769
E
-
1
∖
3.619760
E
-
1
∖
3.262923
E
-
1
∖
2
-
4
2.170599
E
-
1
1.0022
1.825936
E
-
1
0.9873
1.621645
E
-
1
1.0087
2
-
5
1.084135
E
-
1
1.0015
9.068785
E
-
2
1.0097
7.902224
E
-
2
1.0371
2
-
6
5.419942
E
-
2
1.0002
4.528526
E
-
2
1.0019
3.914398
E
-
2
1.0135
2
-
7
2.709797
E
-
2
1.0001
2.261904
E
-
2
1.0015
1.950529
E
-
2
1.0049
2
-
8
1.354897
E
-
2
1.0000
1.131196
E
-
2
0.9997
9.757465
E
-
3
0.9993
2
-
9
6.774435
E
-
3
1.0000
5.656166
E
-
3
1.0000
4.879369
E
-
3
0.9998The 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 Table2, we listed CPU times spent in nonsmooth SOR-Newton algorithm for λ
=
10, 100, and 1000, respectively. In the experiments, we let ϵ
=
1
0
-
8, ω
0
=
0, η
k
=
min
⁡
{
1
0
-
2
/
k
,
∥
F
~
h
(
ω
k
)
∥
}, 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 primal-dual active set method for constrained optimal control problem (in Table 3, “—” indicates that the corresponding CPU time is more than 6000 seconds).Table 2
The CPU times for nonsmooth SOR-Newton algorithm.
h
λ
=
10
λ
=
100
λ
=
1000
2
-
3
0.006091
0.009314
0.010742
2
-
4
0.027895
0.040443
0.047889
2
-
5
0.101071
0.137251
0.179251
2
-
6
0.467282
0.711851
0.717870
2
-
7
4.551622
4.099358
4.886615
2
-
8
92.435878
74.219925
72.793171
2
-
9
1520.338999
1277.017092
1217.694818Table 3
The CPU times for active set method in [5].
h
λ
=
10
λ
=
100
λ
=
1000
2
-
3
0.045277
0.045805
0.046686
2
-
4
0.598815
0.606200
0.825561
2
-
5
7.772854
10.562347
11.056002
2
-
6
165.160697
218.744505
179.321274
2
-
7
—
—
—
2
-
8
—
—
—
2
-
9
—
—
—It follows from Tables2 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 Tables2 and 3 that the smaller h is, the more the CPU time becomes. Indeed, by (46), we have for any v
=
∑
i
=
1
N
n
ν
i
ϕ
i that(61)
0
≤
ν
T
diag
⁡
(
G
j
ω
j
+
|
ω
j
k
meas
⁡
(
Ω
j
)
)
ν
≤
ν
T
diag
⁡
(
meas
⁡
(
Ω
j
)
)
ν
=
ν
T
diag
⁡
(
(
R
h
ϕ
i
,
R
h
ϕ
i
)
)
ν
=
ν
T
(
(
R
h
ϕ
i
,
R
h
ϕ
j
)
)
ν
=
(
R
h
∑
i
=
1
N
h
‍
ν
i
ϕ
i
,
R
h
∑
j
=
1
N
h
‍
ν
j
ϕ
j
)
=
∥
R
h
ν
∥
0
,
Ω
2
≲
∥
ν
∥
0
,
Ω
2
≅
h
2
∥
ν
∥
2
,
where the second equality is from (
R
h
ϕ
i
,
R
h
ϕ
j
)
=
0 for i
≠
j and the “≲” is from (22) and the inverse inequality, while the “≅” is from the following equivalent property between the L
2-norm and the mesh-dependent norm ∥
|
·
|
∥
h under the condition of quasi-uniform grid (see, e.g., [20]):
(62)
∥
v
∥
0
,
Ω
≅
∥
|
ν
|
∥
h
=
Δ
h
∥
ν
∥
.
Estimate (46) together with (61) implies that G
k is an SPD matrix and satisfies
(63)
h
2
≅
λ
min
⁡
(
A
h
)
≤
λ
min
⁡
(
G
k
)
≲
λ
min
⁡
(
A
h
)
+
λ
h
2
≅
(
1
+
λ
)
h
2
,
1
≅
λ
max
⁡
(
A
h
)
≤
λ
max
⁡
(
G
k
)
≲
λ
max
⁡
(
A
h
)
+
λ
h
2
≅
1
+
λ
h
2
.
Therefore, the condition number of G
k satisfies Cond
(
G
k
)
≅
h
-
2. This indicates that the conditioner of the linear system (35) deteriorates as the mesh h 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 Newton-like 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).Table 4
The CPU times for nonsmooth multigrid-Newton algorithm.
h
λ
=
10
λ
=
100
λ
=
1000
2
-
3
0.011277
0.009646
0.013446
2
-
4
0.043442
0.042609
0.040126
2
-
5
0.149155
0.140483
0.163493
2
-
6
0.589568
0.620956
0.806871
2
-
7
2.732603
2.860542
3.987436
2
-
8
15.872829
17.257198
23.665717
2
-
9
93.825685
108.019914
153.901755It follows from the Tables2–4 that compared with nonsmooth SOR-Newton algorithm and ASM, nonsmooth multigrid-Newton algorithm performs excellently for small h.


## 5. 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 isO
(
h
), which is optimal. For the discrete finite element problem, we use nonsmooth Newton-like algorithm to solve it and locally Q-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.



---

*Source: 549305-2014-01-19.xml*