# A D-N Alternating Algorithm for Solving 3D Exterior Helmholtz Problems

Authors: Qing Chen; Baoqing Liu; Qikui Du
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The nonoverlapping domain decomposition method, which is based on the natural boundary reduction, is applied to solve the exterior Helmholtz problem over a three-dimensional domain. The basic idea is to introduce a spherical artificial boundary; the original unbounded domain is changed into a bounded subdomain and a typical unbounded region; then, a Dirichlet-Nuemann (D-N) alternating method is presented; the finite element method and natural boundary element methods are alternately applied to solve the problems in the bounded subdomain and the typical unbounded subdomain. The convergence of the D-N alternating algorithm and its discretization are studied. Some numerical experiments are presented to show the performance of this method.


---

## Body

## 1. Introduction

Many scientific and engineering problems can be reduced to exterior boundary value problems of partial differential equations. Although the numerical methods for solving boundary value problems, such as the finite element method and the finite difference method, are very effective on a bounded domain, yet we often find it difficult to be applied to unbounded problem directly. To solve such problems in infinite region numerically, there are a variety of numerical methods, cf. [1–3] and references therein for more details.There have been studied for Dirichlet-Nuemann alternating algorithm, or D-N alternating algorithm on 2D unbounded domain problems extensively, such as Possion equation [4] and Helmholtz equation [5]. Wu and Yu [6] studied the natural integral equations of 3D Helmholtz problems. Jia et al. [7] investigated the coupled natural boundary element-finite element method for solving 3D exterior Helmholtz problem.In this paper, based on the natural boundary reduction, a D-N alternating algorithm is devised for the numerical solution of three dimensional Helmholtz problem in an infinite region with an inner spheroid boundary. Firstly, we turn the D-N alternating algorithm into the Richardson iterative algorithm which is equivalent to the original method. Secondly, we prove the convergence of algorithm in the general exterior domain. Then, we give the weak form and discretization of the original equations and prove the convergence of the discretization form. Finally, some numerical examples are presented to illustrate the feasibility and efficiency of this method.


## 2. D-N Alternating Algorithm

Let(
r
,
θ
,
φ
) denote the spherical coordinates; we consider the following 3D exterior Helmholtz problem:
(1)
-
Δ
u
-
k
2
u
=
f
,
in
Ω
c
,
∂
u
∂
n
=
g
,
on
⁡
Γ
R
,
where Ω is a bounded domain in R
3, Ω
c
=
R
3
∖
Ω
¯, and Γ
R
=
∂
Ω. k denotes the wave number related to the wavelength λ of the incident wave through k
=
2
π
/
λ. And assume that g is a known function. In order to assure the existence and uniqueness of the solution of (1), problem (1) satisfies the Sommerfeld radiation condition at infinity, which imposes that the scattered wave is outgoing:
(2)
lim
⁡
r
→
∞
r
(
∂
u
∂
r
-
i
k
u
)
=
0
,
where i
=
-
1, r
=
|
x
|
=
x
2
+
y
2
+
z
2. Take a sphere Γ
1
=
{
(
r
,
θ
,
φ
)
∣
r
=
R
1
,
0
⩽
θ
⩽
2
π
,
0
⩽
φ
⩽
π
} which surrounds Γ
R, such that dist
⁡
(
Γ
1
,
Γ
R
)
>
0. Then the artificial boundary Γ
1 divides the original unbounded region into two subregions Ω
ext and Ω
int
⁡, where Ω
ext is the outside subregion and Ω
int
⁡ is the bounded subregion surrounded by Γ
R and Γ
1, Ω
ext and Ω
int
⁡ are nonoverlaping domains. Then the original problem is turned into two subproblems over subdomains Ω
ext and Ω
int
⁡. For exterior boundary value problem (1), we construct the following D-N alternating algorithm:(i)
Step 1. Select the initialλ
0
∈
H
1
/
2
(
Γ
1
), n
=
0.
(ii)
Step 2. Solve the Dirichlet problem inΩ
ext:
(3)
-
Δ
u
2
n
-
k
2
u
2
n
=
0
,
in
Ω
ext
,
u
2
n
=
λ
n
,
on
⁡
Γ
1
,
lim
⁡
r
→
∞
r
(
∂
u
2
n
∂
r
-
i
k
u
2
n
)
=
0
.
(iii)
Step 3. Solve the Neumann problem inΩ
int
⁡:
(4)
-
Δ
u
1
n
-
k
2
u
1
n
=
f
,
in
Ω
int
⁡
,
∂
u
1
n
∂
n
1
=
g
,
on
⁡
Γ
R
,
∂
u
1
n
∂
n
1
=
-
∂
u
2
n
∂
n
2
,
on
⁡
Γ
1
.
(iv)
Step 4. Inputθ
n, and let
(5)
λ
n
+
1
=
θ
n
u
1
n
|
Γ
1
+
(
1
-
θ
n
)
λ
n
.
(v)
Step 5. Letn
=
n
+
1, and go to Step 2,where u
1
n and u
2
n are the nth approximate solutions in Ω
int
⁡ and Ω
ext, respectively. n
1 and n
2 denote the unit outward normals of Γ
1 with respect to two neighboring subdomains; θ
n denotes the nth relaxation factor and λ
0 is an arbitrary function in H
1
/
2
(
Γ
1
). Note that, on interface Γ
1, only the value of the normal derivative of the solution of (3) is needed in solving (4). So it is unnecessary to solve (3). Actually, we can obtain ∂
u
2
n
/
∂
n
2 directly from λ
n by making use of the following natural integral equation [6]:
(6)
∂
u
2
n
∂
n
2
=
K
λ
n
,
where K is the natural integral operator of Helmholtz equation in Ω
ext:
(7)
∂
u
2
n
∂
n
2
=
-
∂
u
2
n
∂
r
(
θ
,
φ
)
=
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
k
h
l
′
(
k
R
1
)
h
l
(
k
R
1
)
Λ
l
m
Y
l
m
(
θ
,
φ
)
,
h
l
(
·
) is the first-kind Hankel function of order l and
(8)
Λ
l
m
=
∫
0
2
π
‍
∫
0
π
‍
λ
n
(
θ
′
,
φ
′
)
Y
l
m
(
θ
′
,
φ
′
)
¯
sin
θ
′
d
θ
′
d
φ
′
,
Y
l
m
(
θ
,
φ
)
=
(
l
-
m
)
!
(
l
+
m
)
!
2
l
+
1
4
π
·
P
l
m
(
cos
⁡
θ
)
e
i
m
φ
,
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
(
1
+
l
2
)
1
/
2
|
Λ
l
m
|
2
<
+
∞
,
and P
l
m
(
·
) denotes the associated Legendre function of the first kind, Y
l
m
(
·
)
¯ denotes the complex conjugate of Y
l
m
(
·
).


## 3. Equivalent Iterative Method

To discuss the convergence of the D-N alternating algorithm in Section2, we first establish a law of deciding. Let u be the exact solution of (1). Denote λ
=
γ
0
u
=
u
|
Γ
1, u
1
=
u
|
Ω
int
⁡, u
2
=
u
|
Ω
ext; it is easy to know λ
∈
H
1
/
2
(
Γ
1
), and u
1 and u
2 are solutions of the following equations, respectively:
(9)
-
Δ
u
1
-
k
2
u
1
=
f
,
in
Ω
int
⁡
,
∂
u
1
∂
n
1
=
g
,
on
⁡
Γ
R
,
u
1
=
λ
,
on
⁡
Γ
1
,
(10)
-
Δ
u
2
-
k
2
u
2
=
0
,
in
Ω
ext
,
u
2
=
λ
,
on
⁡
Γ
1
,
lim
⁡
r
→
∞
r
(
∂
u
2
∂
r
-
i
k
u
2
)
=
0
.
By Green formula, we have
(11)
∂
u
1
∂
n
1
|
Γ
1
=
-
∂
u
2
∂
n
2
|
Γ
1
.
According to superimposition principle, the solutions of (9) and (10) can be denoted
(12)
u
s
=
R
s
λ
+
T
s
f
,
s
=
1,2
,
where T
s
f
(
s
=
1,2
) satisfy, respectively,(13)
-
Δ
T
1
f
-
k
2
T
1
f
=
f
,
in
Ω
int
⁡
,
∂
T
1
f
∂
n
1
=
g
,
on
⁡
Γ
R
,
T
1
f
=
0
,
on
⁡
Γ
1
,
-
Δ
T
2
f
-
k
2
T
2
f
=
0
,
in
Ω
ext
,
T
2
f
=
0
,
on
⁡
Γ
1
,
lim
⁡
r
→
∞
r
(
∂
T
2
f
∂
r
-
i
k
T
2
f
)
=
0
.
It is easy to verify that
(14)
∂
T
2
f
∂
n
1
=
0
.
Extend elements of Ω
1 and Ω
2 by zeros to R
1
λ and R
2
λ, respectively; therefore R
1
λ satisfies
(15)
-
Δ
R
1
λ
-
k
2
R
1
λ
=
0
,
in
Ω
int
⁡
,
∂
R
1
λ
∂
n
1
=
0
,
on
⁡
Γ
R
,
R
1
λ
=
λ
,
on
⁡
Γ
1
,
while R
2
λ satisfies
(16)
-
Δ
R
2
λ
-
k
2
R
2
λ
=
0
,
in
Ω
ext
,
R
2
λ
=
λ
,
on
⁡
Γ
1
,
lim
⁡
r
→
∞
r
(
∂
R
2
λ
∂
r
-
i
k
R
2
λ
)
=
0
,
where R
2 is the Possion integral operator of Helmholtz equation in the exterior spherical domain Ω
ext. Define
(17)
S
i
=
∂
R
i
∂
n
i
,
i
=
1,2
,
S
2 is K, which is the natural integral operator of Helmholtz equation in the exterior spherical domain Ω
ext. Substituting (12) in (11), we have
(18)
∂
(
R
1
-
R
2
)
λ
∂
n
1
=
∂
(
T
2
-
T
1
)
f
∂
n
1
=
-
∂
T
1
f
∂
n
1
.
Namely,
(19)
∂
(
R
1
-
R
2
)
λ
∂
n
1
=
-
∂
T
1
f
∂
n
1
.
Define the right-hand term in (19) χ
=
-
∂
T
1
f
/
∂
n
1; obviously, χ, a function independent of λ, can be solved independently by the problem in subdomains in advance. Then, (19) is equivalent to equation
(20)
S
λ
=
χ
,
where S
=
S
1
+
S
2
=
S
1
+
K is the Poincaré-Steklov operator [8] on Γ
1. Note that λ
∈
H
1
/
2
(
Γ
1
), S
=
S
1
+
S
2
=
S
1
+
K; it comes that
(21)
S
1


:
H
1
/
2
(
Γ
1
)
⟶
H
-
1
/
2
(
Γ
1
)
.
Then,
(22)
S


:
H
1
/
2
(
Γ
1
)
⟶
H
-
1
/
2
(
Γ
1
)
.Theorem 1.
S
i
:
H
1
/
2
(
Γ
1
)
→
H
-
1
/
2
(
Γ
1
)
,
i
=
1,2, is symmetric positive operator.Proof.
We note thatH
1
/
2
(
Γ
1
)
⊂
L
2
(
Γ
1
)
⊂
H
-
1
/
2
(
Γ
1
); therefore, for any λ
,
μ
∈
H
1
/
2
(
Γ
1
)
⊂
L
2
(
Γ
1
). By Green formula, we know
(23)
〈
S
i
λ
,
μ
〉
=
∫
Γ
1
‍
S
i
λ
μ
d
s
=
∫
Γ
1
‍
∂
R
i
λ
∂
n
i
μ
d
s
=
a
i
(
R
i
λ
,
R
i
μ
)
,
where
(24)
a
i
(
u
,
v
)
=
∑
s
=
1
3
‍
∫
Ω
ext
‍
∂
u
∂
x
s
·
∂
v
∂
x
s
d
x
+
∫
Ω
ext
‍
k
2
u
v
d
x
.
This proves S
i is symmetric positive operator.Corollary 2.
There existsS
i
-
1, which is the inverse operator of S
i.Corollary 3.
The bilinear forma
i
(
·
,
·
) on a normed linear space V
k is said to be coercive on H
1
/
2
(
Γ
1
); that is, there exists a positive constant C such that
(25)
〈
S
k
λ
,
λ
〉
⩾
C
∥
λ
∥
H
1
/
2
(
Γ
1
)
2
,
where V
k
=
{
v
∣
v
∈
H
0
1
(
Ω
k
)
,
v
|
∂
Ω
∩
∂
Ω
k
=
0
,
k
=
1,2
,
}.By Poincaré-Steklov operator, solving the original problem (1) is rewritten as solving the operator equation (20). Generally speaking, it is easier to solve S
1
-
1 than to solve the inverse of S. In fact,
(26)
S
1
λ
=
h
,
h
∈
H
-
1
/
2
(
Γ
1
)
is equivalent to finding u
1
∈
V
1, such that
(27)
a
1
(
u
1
,
v
)
=
∫
Γ
1
‍
h
v
d
s
,
∀
v
∈
V
1
.
Let
(28)
λ
=
u
1
|
Γ
1
.
It is easy to verify that
(29)
λ
=
S
1
-
1
h
∈
H
1
/
2
(
Γ
1
)
.
Now we introduce S
1 as preprocessor and preconditioned Richardson iterative method. Given λ
0
∈
H
1
/
2
(
Γ
1
), the procedure of iterative
(30)
S
1
(
λ
n
+
1
-
λ
n
)
=
θ
n
(
χ
-
S
λ
n
)
,
where
(31)
χ
=
-
∂
∂
n
1
(
T
1
f
)
.
Following [9], to analyze the convergence of algorithm (30), what we need to do is just to estimate the eigenvalues of operator S
1
-
1
S, namely, to estimate the upper bound and the lower bound of 〈
S
λ
,
λ
〉
/
〈
S
1
λ
,
λ
〉.Lemma 4.
There exist two positive constantsσ and τ, such that
(32)
σ
=
inf
⁡
λ
∈
H
1
/
2
(
Γ
1
)
〈
S
λ
,
λ
〉
〈
S
1
λ
,
λ
〉
,
τ
=
sup
⁡
λ
∈
H
1
/
2
(
Γ
1
)
〈
S
λ
,
λ
〉
〈
S
1
λ
,
λ
〉
.Lemma 5.
For anyλ
∈
H
1
/
2
(
Γ
1
), one has σ
⩽
〈
S
λ
,
λ
〉
/
〈
S
1
λ
,
λ
〉
⩽
τ.Theorem 6.
If one putsɛ
=
2
/
τ and 0
<
θ
n
<
ɛ, preconditioned Richardson iterative algorithm (30) is convergent. Particularly, optimal relaxation factor is θ
o
p
t
=
2
/
(
σ
+
τ
). Correspondingly, optimal compression ratio of the iteration is δ
o
p
t
=
(
τ
-
σ
)
/
(
σ
+
τ
).Theorem 7.
The D-N alternating algorithm (3)–(5) is equivalent to the preconditioned Richardson iterative method (30).Proof.
Lete
1
n
=
u
1
-
u
1
n, let e
2
n
=
u
2
-
u
2
n, and let v
n
=
λ
-
λ
n. Following (3)–(10), we have
(33)
-
Δ
e
2
n
-
k
2
e
2
n
=
0
,
in
Ω
ext
,
e
2
n
=
v
n
,
on
⁡
Γ
1
,
lim
⁡
r
→
∞
r
(
∂
e
2
n
∂
r
-
i
k
e
2
n
)
=
0
,
(34)
-
Δ
e
1
n
-
k
2
e
1
n
=
0
,
on
⁡
Ω
int
⁡
,
e
1
n
=
0
,
on
⁡
Γ
R
,
∂
e
1
n
∂
n
1
=
-
∂
e
2
n
∂
n
2
,
on
⁡
Γ
1
,
(35)
v
n
+
1
=
θ
n
e
1
n
|
Γ
1
+
(
1
-
θ
n
)
v
n
.
Meanwhile,
(36)
e
2
n
=
R
2
v
n
,
e
1
n
=
R
1
(
e
1
n
)
.
Substituting the above equation to the third equation of (34), we can obtain
(37)
S
1
(
e
1
n
)
=
-
S
2
(
v
n
)
=
-
K
(
v
n
)
.
Therefore,
(38)
S
1
(
v
n
+
1
-
v
n
)
=
S
1
(
θ
n
e
1
n
|
Γ
1
-
θ
n
v
n
)
=
-
θ
n
(
K
v
n
)
-
θ
n
(
S
1
v
n
)
=
-
θ
n
S
v
n
;
thus,
(39)
S
1
(
λ
n
+
1
-
λ
n
)
=
θ
n
S
(
λ
-
λ
n
)
.
It comes that
(40)
S
λ
=
∂
∂
n
1
(
R
1
λ
)
+
∂
∂
n
2
(
R
2
λ
)
=
∂
∂
n
1
(
u
1
-
T
1
f
)
+
∂
u
2
∂
n
2
=
-
∂
T
1
f
∂
n
1
=
χ
.
Substituting the above equation in (39), we have
(41)
S
1
(
λ
n
+
1
-
λ
n
)
=
θ
n
(
χ
-
S
λ
n
)
.


## 4. Convergence of the Algorithm

Considering a special case in whichΓ
0
=
{
(
r
,
θ
,
φ
)
∣
r
=
R
0
,
0
⩽
θ
⩽
2
π
,
0
⩽
φ
⩽
π
} and the interface Γ
1
=
{
(
r
,
θ
,
φ
)
∣
r
=
R
1
,
0
⩽
θ
⩽
2
π
,
0
⩽
φ
⩽
π
}, following [10], we have the following.Lemma 8.
Suppose(42)
h
h
h
h
f
n
(
x
)
=
n
2
n
+
1
+
n
+
1
2
n
+
1
x
n
+
1
,
x
∈
[
0,1
]
,
n
∈
Z
+
=
{
0,1
,
2
,
…
}
,
and set
(43)
U
(
x
)
=
sup
⁡
n
∈
Z
+
{
f
n
(
x
)
}
,
L
(
x
)
=
inf
⁡
n
∈
Z
+
{
f
n
(
x
)
}
;
then, we have
(44)
U
(
x
)
=
{
1
2
,
0
⩽
x
⩽
1
2
,
x
,
1
2
<
x
⩽
1
,
L
(
x
)
=
{
f
0
(
x
)
,
0
⩽
x
⩽
x
0
,
f
k
+
1
(
x
)
,
x
k
<
x
⩽
x
k
+
1
,
where x
k is the only root of equation f
k
+
1
(
x
)
-
f
k
(
x
)
=
0 in interval (
0,1
) and x
k
<
x
k
+
1, k is a nonnegative integer.Forλ
(
θ
,
φ
)
∈
H
1
/
2
(
Γ
1
), in the L
2-sense, the following expansion holds
(45)
λ
(
θ
,
φ
)
=
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
Λ
l
m
Y
l
m
(
θ
,
φ
)
,
where P
l
m
(
x
) denotes the associated Legendre function of the first kind; Y
l
m
(
θ
,
φ
)
¯ denotes the complex conjugate of Y
l
m. On the one hand,
(46)
〈
K
c
λ
,
λ
〉
=
R
1
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
(
1
+
l
)
|
Λ
l
m
|
2
satisfies
(47)
R
1
{
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
(
1
+
l
2
)
1
/
2
|
Λ
l
m
|
2
}
1
/
2
<
+
∞
;
Λ
l
m
¯
=
(
-
1
)
|
m
|
Λ
l
-
m
,
on the other hand, set
(48)
H
1
N
λ
=
∑
l
=
0
N
‍
∑
m
=
-
l
l
‍
R
1
l
+
1
Λ
l
m
R
1
2
l
+
1
-
R
0
2
l
+
1
(
r
l
-
R
0
2
l
+
1
r
l
+
1
)
Y
l
m
(
θ
,
φ
)
,
h
h
h
h
h
h
h
h
h
h
h
h
f
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
R
0
<
r
<
R
1
.
Following [11] and (48), we have
(49)
lim
⁡
N
→
∞
∥
H
1
λ
-
H
1
N
λ
∥
H
1
(
Ω
1
)
=
0
.
According to Green formula and orthogonality of Y
l
m
(
θ
,
φ
), it comes that
(50)
〈
S
1
λ
,
λ
〉
=
∫
Ω
1
‍
∇
H
1
λ
·
∇
H
1
λ
d
x
=
lim
⁡
M
→
∞
,
N
→
∞
∫
Ω
1
‍
∇
H
1
M
λ
·
∇
H
1
M
λ
d
x
=
lim
⁡
M
→
∞
,
N
→
∞
∫
Γ
1
‍
∂
H
1
M
λ
∂
r
H
1
N
λ
d
s
=
lim
⁡
M
→
∞
,
N
→
∞
∑
l
=
0
min
⁡
{
M
,
N
}
‍
∑
m
=
-
l
l
‍
l
R
1
2
l
+
1
+
(
l
+
1
)
R
0
2
l
+
1
R
1
2
l
+
1
-
R
0
2
l
+
1
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
R
1
Λ
l
m
·
Λ
l
m
¯
.Lemma 9.
Suppose(51)
〈
S
1
λ
,
λ
〉
=
R
1
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
l
R
1
2
l
+
1
+
(
l
+
1
)
R
0
2
l
+
1
R
1
2
l
+
1
-
R
0
2
l
+
1
|
Λ
l
m
|
2
;
〈
K
λ
,
λ
〉
=
R
1
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
(
l
+
1
)
|
Λ
l
m
|
2
,
then,
(52)
1
U
(
R
0
/
R
1
)
⩽
〈
S
λ
,
λ
〉
〈
S
1
λ
,
λ
〉
⩽
1
L
(
R
0
/
R
1
)
.Theorem 10.
If relaxation factorθ
n satisfies
(53)
0
<
θ
n
<
2
L
(
R
0
R
1
)
,
then the D-N alternating algorithm (3)–(5) converges. Particularly, the optimal relaxation factor is
(54)
θ
o
p
t
=
2
U
(
R
0
/
R
1
)
L
(
R
0
/
R
1
)
U
(
R
0
/
R
1
)
+
L
(
R
0
/
R
1
)
,
while the optimal compression ratio of the iteration is
(55)
δ
o
p
t
=
U
(
R
0
/
R
1
)
-
L
(
R
0
/
R
1
)
U
(
R
0
/
R
1
)
+
L
(
R
0
/
R
1
)
.Remark 11.
We have analyzed the convergence for exterior spherical domain. Following [10], the convergence analysis for the general exterior domain can be extended from above content. Obviously, the optimal relaxation factor and optimal compression ratio of the iteration depend on R
1 and the geometry of Ω
c in some way, which is confirmed by our numerical results.


## 5. The Weak Form and Discretization

The weak form of problem (4) is Find u
1
n
∈
H
g
1
(
Ω
int
⁡
), such that
(56)
D
Ω
int
⁡
(
u
1
n
,
v
)
=
∫
Ω
int
⁡
‍
f
v
d
x
-
∫
Γ
1
‍
v
∂
u
2
n
∂
n
2
d
s
,
∀
v
∈
H
0
1
(
Ω
int
⁡
)
,
where
(57)
D
Ω
int
⁡
(
u
1
n
,
v
)
=
∫
Ω
int
⁡
‍
(
∇
u
1
n
·
∇
v
+
k
2
u
1
n
v
)
d
x
,
H
g
1
(
Ω
int
⁡
)
=
{
v
∣
v
∈
H
1
(
Ω
int
⁡
)
,
v
|
Γ
R
=
g
}
,
H
0
1
(
Ω
int
⁡
)
=
{
v
∣
v
∈
H
1
(
Ω
int
⁡
)
,
v
|
Γ
0
=
0
}
.Obviously,v
|
Γ
1
∈
H
1
/
2
(
Γ
1
). Following (6), the weak form (56) can be rewritten as Find u
1
n
∈
H
g
1
(
Ω
int
⁡
), such that
(58)
D
Ω
int
⁡
(
u
1
n
,
v
)
=
∫
Ω
int
⁡
‍
f
v
d
x
-
D
^
c
(
λ
n
,
v
)
,
∀
v
∈
H
0
1
(
Ω
int
⁡
)
,
i
=
1,2
,
where
(59)
D
^
c
(
λ
n
,
v
)
=
〈
K
λ
n
,
v
〉
=
-
∑
l
=
0
∞
‍
∑
m
=
-
l
l
‍
k
R
1
2
h
l
′
(
k
R
)
h
l
(
k
R
)
(
λ
n
)
l
m
V
l
m
¯
,
V
l
m
=
∫
0
2
π
‍
∫
0
π
‍
v
(
θ
′
,
φ
′
)
Y
l
m
(
θ
′
,
φ
′
)
¯
sin
θ
′
d
θ
′
d
φ
′
,
(
λ
n
)
l
m
=
∫
0
2
π
‍
∫
0
π
‍
λ
n
(
θ
′
,
φ
′
)
Y
l
m
(
θ
′
,
φ
′
)
sin
θ
′
d
θ
′
d
φ
′
.
(58) and (5) are the weak form of D-N method (3)–(5). In the following, we consider the discretization of this weak form. Discretize Ω
int
⁡ into a finite number of element domains. Let S
h
(
Ω
int
⁡
) denote the linear subspace of H
0
1
(
Ω
int
⁡
) corresponding to this partition. Define M
0, M
1, and M
i to be the sets of all nodes belonging to Γ
R, Γ
1, and Ω
int
⁡, respectively. Denote U
A
n the approximate value and L
A
(
x
) a basis function at node A. Obviously,
(60)
S
h
(
Ω
int
⁡
)
=
span
⁡
{
L
A
i
(
x
)
,
L
A
1
(
x
)
∣
A
i
∈
M
i
,
A
1
∈
M
1
}
.
The finite element approximation of u
1
n can be expressed as
(61)
u
1
h
n
=
∑
A
i
∈
M
i
‍
U
A
i
n
L
A
i
(
x
)
+
∑
A
1
∈
M
1
‍
U
A
1
n
L
A
1
(
x
)
+
∑
A
0
∈
M
0
‍
g
(
x
A
0
)
L
A
0
(
x
)
.
By (61) and (5), we can obtain the following discrete D-N alternating algorithm:
(62)
(
Q
11
Q
1
i
Q
10
Q
i
1
Q
i
i
Q
i
0
Q
01
Q
0
i
Q
00
)
(
U
1
n
U
i
n
U
0
n
)
=
(
-
B
Λ
n
0
F
0
)
,
(63)
Λ
n
+
1
=
θ
n
U
A
1
n
+
(
1
-
θ
n
)
Λ
n
,
where
(64)
Λ
n
=
[
λ
n
(
A
1
)
]
T
,
U
s
n
=
[
U
A
s
n
]
T
,
g
g
g
g
g
g
g
g
Q
s
t
=
[
D
Ω
int
⁡
(
L
A
s
,
L
B
t
)
]
,
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
s
,
t
=
0,1
,
i
,
Q
01
=
Q
10
=
0
,
B
=
[
b
(
γ
0
N
A
1
,
γ
0
N
B
1
)
]
,
g
g
g
g
g
g
g
g
g
F
0
=
[
∫
Γ
R
‍
g
N
B
0
d
s
]
T
,
g
g
g
g
g
g
g
A
s
,
B
s
∈
M
s
,
s
=
0,1
,
i
.


## 6. Convergence Analysis for the Discretization Form

In the following, we consider the convergence of the discretization form.Theorem 12.
The discrete D-N alternating algorithms (62) and (63) are equivalent to the following preconditioned Richardson iteration:
(65)
S
h
(
1
)
(
Λ
n
+
1
-
Λ
n
)
=
θ
n
(
F
1
-
S
h
Λ
n
)
,
where
(66)
S
h
(
1
)
=
Q
11
-
Q
1
i
(
Q
i
i
-
Q
i
0
(
Q
00
)
-
1
Q
0
i
)
-
1
Q
i
1
,
S
h
=
S
h
(
1
)
+
B
,
F
1
=
Q
1
i
(
Q
i
i
-
Q
i
0
(
Q
00
)
-
1
Q
0
i
)
-
1
Q
i
0
(
Q
00
)
-
1
F
0
.Proof.
Let[
U
1
,
U
i
,
U
0
]
T be the solution of the equation
(67)
(
Q
11
+
B
Q
1
i
Q
10
Q
i
1
Q
i
i
Q
i
0
Q
01
Q
0
i
Q
00
)
(
U
1
U
i
U
0
)
=
(
0
0
F
0
)
.
Simplified to
(68)
[
Q
11
-
Q
1
i
(
Q
i
i
-
Q
i
0
(
Q
00
)
-
1
Q
0
i
)
-
1
Q
i
1
+
B
]
U
1
=
Q
1
i
(
Q
i
i
-
Q
i
0
(
Q
00
)
-
1
Q
0
i
)
-
1
Q
i
0
(
Q
00
)
-
1
F
0
,
namely,
(69)
S
h
U
1
=
F
1
.
Then, (67) can be rewritten as
(70)
(
Q
11
Q
1
i
Q
10
Q
i
1
Q
i
i
Q
i
0
Q
01
Q
0
i
Q
00
)
(
U
1
U
i
U
0
)
=
(
-
B
U
1
0
F
0
)
.
By (62) and (70),
(71)
(
Q
11
Q
1
i
Q
10
Q
i
1
Q
i
i
Q
i
0
Q
01
Q
0
i
Q
00
)
(
U
1
n
-
U
1
U
i
n
-
U
i
U
0
n
-
U
0
)
=
(
-
B
(
Λ
n
-
U
1
)
0
0
)
.
We have
(72)
S
h
(
1
)
(
U
1
n
-
U
1
)
=
B
(
U
1
-
Λ
n
)
.
Therefore,
(73)
S
h
(
1
)
(
Λ
n
+
1
-
Λ
n
)
=
θ
n
S
h
(
1
)
(
U
1
n
-
Λ
n
)
=
θ
n
[
S
h
(
1
)
(
U
1
n
-
U
1
)
+
S
h
(
1
)
(
U
1
-
Λ
n
)
]
=
θ
n
(
B
+
S
h
(
1
)
)
(
U
1
-
Λ
n
)
=
θ
n
S
h
(
U
1
-
Λ
n
)
=
θ
n
(
F
1
-
S
h
Λ
n
)
.Theorem 13.
Letα be spectral radius of (
S
h
(
1
)
)
-
1
S
h, which is iterative matrix of preconditioned Richardson iteration. Then, there is a positive constant β, which is independent of finite element mesh parameter h of subdomain Ω
int
⁡, such that α
⩽
β.Theorem 14.
Eigenvalues of iterative matrix(
S
h
(
1
)
)
-
1
S
h are real and greater than 1.Theorem 15.
Putθ
n
=
η
(
n
=
0,1
,
…
); then, there exists a constant δ
(
0
<
δ
<
1
), which is independent of finite element mesh parameter h of subdomain Ω
int
⁡. For 0
<
η
<
δ, the preconditioned Richardson iteration converges, the discrete D-N alternating algorithm (62)-(63) converges and the convergence rate is independent of mesh parameter h of subdomain Ω
int
⁡.Proof.
By (73), we have
(74)
U
1
-
Λ
n
+
1
=
[
I
-
η
(
S
h
(
1
)
)
-
1
S
h
]
(
U
1
-
Λ
n
)
=
[
I
-
η
(
S
h
(
1
)
)
-
1
S
h
]
n
+
1
(
U
1
-
Λ
0
)
;
it comes that
(75)
∥
U
1
-
Λ
n
+
1
∥
2
⩽
δ
n
+
1
∥
U
1
-
Λ
0
∥
2
,
where
(76)
δ
=
∥
I
-
η
(
S
h
(
1
)
)
-
1
S
h
∥
2
.
Following Theorem 13 and Theorem 14, there exists a constant δ
(
0
<
δ
<
1
), which is independent of h. For 0
<
η
<
σ, spectral radius of I
-
η
(
S
h
(
1
)
)
-
1
S
h is less than 1, and spectral norm δ
<
1; therefore,
(77)
lim
⁡
n
→
+
∞
∥
U
1
-
Λ
n
+
1
∥
2
=
0
.
It follows that the preconditioned Richardson iteration converges; then, the discrete D-N alternating algorithm (62)-(63) converges and the convergence rate is independent of mesh parameter h of subdomain Ω
int
⁡.


## 7. Numerical Examples

To test the effectiveness of the method in this paper, we give two numerical examples, using the discrete D-N alternating algorithm in Section2. In the two examples, the exact solutions are known. The purpose of these examples is to check the convergence in terms of iteration and mesh size h.Suppose thatΩ
c
=
{
(
x
,
y
,
z
)
∣
|
x
|
>
1
,
|
y
|
>
1
,
|
z
|
>
1
} is the exterior domain of cube [
-
1,1
]
×
[
-
1,1
]
×
[
-
1,1
]. The artificial boundary is Γ
1
=
{
(
r
,
θ
,
φ
)
∣
r
=
R
1
,
R
1
>
3
}, r
=
x
2
+
y
2
+
z
2. Split the bounded domain Ω
int
⁡ into finite element mesh as follows. Firstly, make every inner boundary edge into N
1 equivalent parts and connect nodes by using segments which are parallel to the corresponding coordinate axis; then, the grid located in the inner boundary Γ
1 is generated. Secondly, along the radial direction, rays starting from each node of located in the inner boundary intersect artificial boundaries at some nodes, and all the other nodes are produced by dividing segments along the radial direction into N
2 equivalent parts. At last, we get an eight-node trilinear isoparametric finite elements of Ω
int
⁡.Denotee the maximum node-error on Ω
¯
i:
(78)
e
(
n
)
=
sup
⁡
P
i
∈
Ω
¯
i
|
u
(
P
i
)
-
u
1
h
n
(
P
i
)
|
;
e
h denotes the maximum node-error of the abjacent two-steps on nodes
(79)
e
h
(
n
)
=
sup
⁡
P
i
∈
Ω
¯
i
|
u
1
h
n
+
1
(
P
i
)
-
u
1
h
n
(
P
i
)
|
;
q
h is the convergence rate
(80)
q
h
(
n
)
=
e
h
(
n
-
1
)
e
h
(
n
)
.We substitute∑
l
=
0
M
‍ for ∑
l
=
0
∞
‍ in the computing of the entries of K
h. Denote by m the total number of nodes on Ω
int
⁡
¯. By computing, the results are as follows.Example 16.
In the problem (1), we take
(81)
g
=
3
4
π
(
k
2
i
-
2
i
-
2
k
k
2
)
e
i
k
cos
⁡
θ
;
the exact solution is
(82)
u
(
r
,
θ
)
=
-
3
4
π
k
r
+
i
k
2
r
2
e
i
k
r
cos
⁡
θ
.
Let k
=
2.0, N
1
=
N
2
=
2, m
=
78, M
=
4, the numerical results are shown in Tables 1 and 2.
Tables1 and 2 show that the D-N alternating algorithm (62)-(63) is convergent. Table 1 shows that the speed of convergence is faster while θ
n
=
0.5. Conversely, Table 2 shows the different phenomenon that the speed of convergence is faster while θ
n
=
0.3. It proves Theorem 10 and shows relaxation factor θ
n is dependent on R
1. But it is valuable to know that θ
n
=
0.5 unsatisfied the requirement of Theorem 10  (
2
L
(
R
0
/
R
1
)
=
0.5
), but the D-N alternating algorithm (62)-(63) is convergent, which shows (53) is the sufficient condition of Theorem 10.Table 1
Maximum node-error onΩ
¯
i for R
1
=
2.0.
θ
Iterative times
0
1
2
3
4
12
0.3
0.33062
0.18917
0.11558
0.07835
0.05963
0.04858
0.4
0.33062
0.14107
0.07438
0.05198
0.04806
0.04713
0.5
0.33062
0.08362
0.05926
0.04816
0.04823
0.04818
0.6
0.33062
0.06345
0.09187
0.06537
0.06823
0.09684
0.7
0.33062
0.11998
0.15164
0.06356
0.08451
0.06529Table 2
Maximum node-error onΩ
¯
i for R
1
=
6.0.
θ
n
Iterative times
0
1
2
3
4
12
0.3
0.13381
0.05442
0.04968
0.04930
0.04927
0.04927
0.4
0.13381
0.03435
0.05269
0.04849
0.49448
0.04927
0.5
0.13381
0.04332
0.06801
0.03922
0.05465
0.04939
0.6
0.13381
0.11119
0.14317
0.07156
0.11630
0.08039
0.7
0.13381
0.16126
0.24338
0.22273
0.31370
0.63435Example 17.
In the problem (1), we take
(83)
g
=
15
32
π
(
k
3
+
4
k
2
i
-
9
k
-
9
i
k
3
)
e
i
(
k
+
ϕ
)
sin
2
θ
;
the exact solution is
(84)
u
(
r
,
θ
)
=
-
15
32
π
(
3
i
k
3
r
3
+
3
k
2
r
2
-
i
k
r
)
e
i
(
k
r
+
φ
)
sin
2
θ
.
Let k
=
2.0, R
1
=
6.0, N
1
=
N
2
=
2, m
=
78, M
=
4, the numerical results are shown in Table 3.
Table3 shows that in the scope of relaxation factor, the bigger the relaxation factor θ
n is, the faster the convergence rate is.
While we putθ
n
=
0.4 and R
1
=
6.0 and denote Mesh I: N
1
=
N
2
=
1, m
=
16, M
=
4; Mesh II: N
1
=
N
2
=
2, m
=
78, M
=
4; Mesh III: N
1
=
N
2
=
4, m
=
490, M
=
4. The numerical results are shown in Table 4.
Table4 shows that the convergence rate of the discrete D-N alternating algorithm is independent of mesh parameter h. At the same time, the closer we refine the mesh, the smaller the maximum node-error will be. The maximum node-error is roughly of O
(
h
) order.Table 3
The relation between convergence rate and relaxation factor.
θ
n
Iterative times
0
1
2
3
4
12
0.2
e
0.13381
0.07449
0.05900
0.05302
0.05072
0.04928
e
h
0.02358
0.05223
0.01157
0.00256
0.00057
q
h
4.51432
4.51329
4.51206
4.51101
0.3
e
0.13381
0.54416
0.04968
0.04330
0.04927
0.04927
e
h
0
0.50323
0.0806
0.01294
0.00209
0.00031
q
h
0
0
6.2432
6.2301
6.1952
6.1809
0.4
e
0.13381
0.03434
0.05268
0.04849
0.04945
0.04927
e
h
0
0.01512
0.00176
0.00021
0.00002
0
q
h
0
0
8.5961
8.5532
8.3201
7.9326Table 4
The relation between convergence rate and mesh.
Mesh
Iterative times
0
1
2
3
4
12
I
e
0.16651
0.12025
0.11342
0.11354
0.11339
0.11327
e
h
0.11547
0.10300
0.08651
0.07381
0.06525
q
h
1.12105
1.19058
1.17205
1.13114
II
e
0.13381
0.03434
0.05268
0.04849
0.04945
0.04927
e
h
0.02834
0.01814
0.01162
0.00745
0.00478
q
h
1.56223
1.56134
1.56001
1.56018
III
e
0.08874
0.01239
0.01753
0.01623
0.01648
0.01644
e
h
0.01056
0.00449
0.00191
0.00083
0.00036
q
h
2.35073
2.35327
2.30124
2.29801



---

*Source: 418426-2014-07-20.xml*