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

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.

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 (𝑟,𝜃,𝜑) denote the spherical coordinates; we consider the following 3D exterior Helmholtz problem: −Δ𝑢−𝑘²𝑢=𝑓,inΩ^𝑐,𝜕𝑢𝜕𝑛=𝑔,onΓ_𝑅, (1) where Ω is a bounded domain in ℝ³, Ω^𝑐=ℝ³⧵Ω, and Γ_𝑅=𝜕Ω. 𝑘 denotes the wave number related to the wavelength 𝜆 of the incident wave through 𝑘=2𝜋/𝜆. And assume that 𝑔 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: lim_𝑟→∞𝑟𝜕𝑢𝜕𝑟−i𝑘𝑢=0, (2) where i=√−1, 𝑟=|𝐱|=𝑥²+𝑦²+𝑧². Take a sphere Γ₁={(𝑟,𝜃,𝜑)∣𝑟=𝑅₁,0⩽𝜃⩽2𝜋,0⩽𝜑⩽𝜋} which surrounds Γ_𝑅, such that dist(Γ₁,Γ_𝑅)>0. Then the artificial boundary Γ₁ 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 Γ_𝑅 and Γ₁, Ω_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 𝜆⁰∈𝐻^1/2(Γ₁), 𝑛=0.(ii)Step 2. Solve the Dirichlet problem in Ω_ext: −Δ𝑢^𝑛₂−𝑘²𝑢^𝑛₂=0,inΩ_ext,𝑢^𝑛₂=𝜆^𝑛,onΓ₁,lim_𝑟→∞𝑟𝜕𝑢^𝑛₂𝜕𝑟−i𝑘𝑢^𝑛₂=0. (3)(iii)Step 3. Solve the Neumann problem in Ω_int: −Δ𝑢^𝑛₁−𝑘²𝑢^𝑛₁=𝑓,inΩ_int,𝜕𝑢^𝑛₁𝜕𝑛₁=𝑔,onΓ_𝑅,𝜕𝑢^𝑛₁𝜕𝑛₁=−𝜕𝑢^𝑛₂𝜕𝑛₂,onΓ₁. (4)(iv)Step 4. Input 𝜃_𝑛, and let 𝜆^𝑛+1=𝜃_𝑛𝑢^𝑛_1Γ1+1−𝜃_𝑛𝜆^𝑛. (5)(v)Step 5. Let 𝑛=𝑛+1, and go to Step 2,where 𝑢^𝑛₁ and 𝑢^𝑛₂ are the 𝑛th approximate solutions in Ω_int and Ω_ext, respectively. 𝑛₁ and 𝑛₂ denote the unit outward normals of Γ₁ with respect to two neighboring subdomains; 𝜃_𝑛 denotes the 𝑛th relaxation factor and 𝜆⁰ is an arbitrary function in 𝐻^1/2(Γ₁). Note that, on interface Γ₁, 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 𝜕𝑢^𝑛₂/𝜕𝑛₂ directly from 𝜆^𝑛 by making use of the following natural integral equation [6]: 𝜕𝑢^𝑛₂𝜕𝑛₂=𝒦𝜆^𝑛, (6) where 𝒦 is the natural integral operator of Helmholtz equation in Ω_ext: 𝜕𝑢^𝑛₂𝜕𝑛₂=−𝜕𝑢^𝑛₂=𝜕𝑟𝜃,𝜑^∞𝑙_{𝑙=0𝑚=−𝑙}𝑘ℎ^′_𝑙𝑘𝑅₁ℎ_𝑙𝑘𝑅₁Λ^𝑚_𝑙𝑌^𝑚_𝑙𝜃,𝜑, (7)ℎ_𝑙(⋅) is the first-kind Hankel function of order 𝑙 and Λ^𝑚_𝑙=₀^2𝜋𝜋₀𝜆^𝑛𝜃^′,𝜑^′𝑌^𝑚_𝑙𝜃^′,𝜑^′sin𝜃^′d𝜃^′d𝜑^′,𝑌^𝑚_𝑙𝜃,𝜑=(𝑙−𝑚)!(𝑙+𝑚)!2𝑙+14𝜋⋅𝑃^𝑚_𝑙(cos𝜃)𝑒^𝑖𝑚𝜑,_{∞𝑙𝑙=0𝑚=−𝑙}1+𝑙^21/2Λ^𝑚_𝑙²<+∞, (8) and 𝑃^𝑚_𝑙(⋅) denotes the associated Legendre function of the first kind, 𝑌^𝑚_𝑙(⋅) denotes the complex conjugate of 𝑌^𝑚_𝑙(⋅).

3. Equivalent Iterative Method

To discuss the convergence of the D-N alternating algorithm in Section 2, we first establish a law of deciding. Let 𝑢 be the exact solution of (1). Denote 𝜆=𝛾₀𝑢=𝑢|_Γ1, 𝑢₁=𝑢|_Ωint, 𝑢₂=𝑢|_Ωext; it is easy to know 𝜆∈𝐻^1/2(Γ₁), and 𝑢₁ and 𝑢₂ are solutions of the following equations, respectively: −Δ𝑢₁−𝑘²𝑢₁=𝑓,inΩ_int,𝜕𝑢₁𝜕𝑛₁=𝑔,onΓ_𝑅,𝑢₁=𝜆,onΓ₁,−Δ𝑢₂−𝑘²𝑢₂=0,inΩ_ext,𝑢₂=𝜆,onΓ₁,lim_𝑟→∞𝑟𝜕𝑢₂𝜕𝑟−i𝑘𝑢₂=0. (9)(10) By Green formula, we have 𝜕𝑢₁𝜕𝑛_1Γ1=−𝜕𝑢₂𝜕𝑛_2Γ1. (11) According to superimposition principle, the solutions of (9) and (10) can be denoted 𝑢_𝑠=ℛ_𝑠𝜆+𝒯_𝑠𝑓,𝑠=1,2, (12) where 𝒯_𝑠𝑓(𝑠=1,2) satisfy, respectively,−Δ𝒯₁𝑓−𝑘²𝒯₁𝑓=𝑓,inΩ_int,𝜕𝒯₁𝑓𝜕𝑛₁=𝑔,onΓ_𝑅,𝒯₁𝑓=0,onΓ₁,−Δ𝒯₂𝑓−𝑘²𝒯₂𝑓=0,inΩ_ext,𝒯₂𝑓=0,onΓ₁,lim_𝑟→∞𝑟𝜕𝒯₂𝑓𝜕𝑟−i𝑘𝒯₂𝑓=0. (13) It is easy to verify that 𝜕𝒯₂𝑓𝜕𝑛₁=0. (14) Extend elements of Ω₁ and Ω₂ by zeros to ℛ₁𝜆 and ℛ₂𝜆, respectively; therefore ℛ₁𝜆 satisfies −Δℛ₁𝜆−𝑘²ℛ₁𝜆=0,inΩ_int,𝜕ℛ₁𝜆𝜕𝑛₁=0,onΓ_𝑅,ℛ₁𝜆=𝜆,onΓ₁, (15) while ℛ₂𝜆 satisfies −Δℛ₂𝜆−𝑘²ℛ₂𝜆=0,inΩ_ext,ℛ₂𝜆=𝜆,onΓ₁,lim_𝑟→∞𝑟𝜕ℛ₂𝜆𝜕𝑟−i𝑘ℛ₂𝜆=0, (16) where ℛ₂ is the Possion integral operator of Helmholtz equation in the exterior spherical domain Ω_ext. Define 𝒮_𝑖=𝜕ℛ_𝑖𝜕𝑛_𝑖,𝑖=1,2, (17)𝒮₂ is 𝒦, which is the natural integral operator of Helmholtz equation in the exterior spherical domain Ω_ext. Substituting (12) in (11), we have 𝜕ℛ₁−ℛ₂𝜆𝜕𝑛₁=𝜕𝒯₂−𝒯₁𝑓𝜕𝑛₁=−𝜕𝒯₁𝑓𝜕𝑛₁. (18) Namely, 𝜕ℛ₁−ℛ₂𝜆𝜕𝑛₁=−𝜕𝒯₁𝑓𝜕𝑛₁. (19) Define the right-hand term in (19) 𝜒=−𝜕𝒯₁𝑓/𝜕𝑛₁; obviously, 𝜒, a function independent of 𝜆, can be solved independently by the problem in subdomains in advance. Then, (19) is equivalent to equation 𝒮𝜆=𝜒, (20) where 𝒮=𝒮₁+𝒮₂=𝒮₁+𝒦 is the Poincaré-Steklov operator [8] on Γ₁. Note that 𝜆∈𝐻^1/2(Γ₁), 𝒮=𝒮₁+𝒮₂=𝒮₁+𝒦; it comes that 𝒮₁∶𝐻^1/2Γ₁⟶𝐻^−1/2Γ₁. (21) Then, 𝒮∶𝐻^1/2Γ₁⟶𝐻^−1/2Γ₁. (22)

Theorem 1. 𝒮_𝑖∶𝐻^1/2(Γ₁)→𝐻^−1/2(Γ₁),𝑖=1,2, is symmetric positive operator.

Proof. We note that 𝐻^1/2(Γ₁)⊂𝐿²(Γ₁)⊂𝐻^−1/2(Γ₁); therefore, for any 𝜆,𝜇∈𝐻^1/2(Γ₁)⊂𝐿²(Γ₁). By Green formula, we know 𝒮_𝑖𝜆,𝜇=_Γ1𝒮_𝑖𝜆𝜇d𝑠=_Γ1𝜕ℛ_𝑖𝜆𝜕𝑛_𝑖𝜇d𝑠=𝑎_𝑖ℛ_𝑖𝜆,ℛ_𝑖𝜇, (23) where 𝑎_𝑖(𝑢,𝑣)=³_𝑠=1Ωext𝜕𝑢𝜕𝑥_𝑠⋅𝜕𝑣𝜕𝑥_𝑠d𝐱+_Ωext𝑘²𝑢𝑣d𝐱. (24) This proves 𝒮_𝑖 is symmetric positive operator.

Corollary 2. There exists 𝒮_𝑖⁻¹, which is the inverse operator of 𝒮_𝑖.

Corollary 3. The bilinear form 𝑎_𝑖(⋅,⋅) on a normed linear space 𝑉_𝑘 is said to be coercive on 𝐻^1/2(Γ₁); that is, there exists a positive constant 𝐶 such that 𝒮_𝑘𝜆,𝜆⩾𝐶‖𝜆‖²_𝐻1/2(Γ1), (25) where 𝑉_𝑘={𝑣∣𝑣∈𝐻¹₀(Ω_𝑘),𝑣|_{𝜕Ω∩𝜕Ω𝑘}=0,𝑘=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 𝒮₁⁻¹ than to solve the inverse of 𝒮. In fact, 𝒮₁𝜆=ℎ,ℎ∈𝐻^−1/2Γ₁ (26) is equivalent to finding 𝑢₁∈𝑉₁, such that 𝑎₁𝑢₁,𝑣=_Γ1ℎ𝑣d𝑠,∀𝑣∈𝑉₁. (27) Let 𝜆=𝑢₁|_Γ1. (28) It is easy to verify that 𝜆=𝒮₁⁻¹ℎ∈𝐻^1/2Γ₁. (29) Now we introduce 𝒮₁ as preprocessor and preconditioned Richardson iterative method. Given 𝜆⁰∈𝐻^1/2(Γ₁), the procedure of iterative 𝒮₁𝜆^𝑛+1−𝜆^𝑛=𝜃_𝑛𝜒−𝒮𝜆^𝑛, (30) where 𝜕𝜒=−𝜕𝑛₁𝒯₁𝑓. (31) Following [9], to analyze the convergence of algorithm (30), what we need to do is just to estimate the eigenvalues of operator 𝒮₁⁻¹𝒮, namely, to estimate the upper bound and the lower bound of ⟨𝒮𝜆,𝜆⟩/⟨𝒮₁𝜆,𝜆⟩.

Lemma 4. There exist two positive constants 𝜎 and 𝜏, such that 𝜎=inf_{𝜆∈𝐻1/2(Γ1)}⟨𝒮𝜆,𝜆⟩𝒮₁𝜆,𝜆,𝜏=sup_{𝜆∈𝐻1/2(Γ1)}⟨𝒮𝜆,𝜆⟩𝒮₁.𝜆,𝜆 (32)

Lemma 5. For any 𝜆∈𝐻^1/2(Γ₁), one has 𝜎⩽⟨𝒮𝜆,𝜆⟩/⟨𝒮₁𝜆,𝜆⟩⩽𝜏.

Theorem 6. If one puts 𝜀=2/𝜏 and 0<𝜃_𝑛<𝜀, preconditioned Richardson iterative algorithm (30) is convergent. Particularly, optimal relaxation factor is 𝜃_𝑜𝑝𝑡=2/(𝜎+𝜏). Correspondingly, optimal compression ratio of the iteration is 𝛿_𝑜𝑝𝑡=(𝜏−𝜎)/(𝜎+𝜏).

Theorem 7. The D-N alternating algorithm (3)–(5) is equivalent to the preconditioned Richardson iterative method (30).

Proof. Let 𝑒^𝑛₁=𝑢₁−𝑢^𝑛₁, let 𝑒^𝑛₂=𝑢₂−𝑢^𝑛₂, and let 𝑣^𝑛=𝜆−𝜆^𝑛. Following (3)–(10), we have −Δ𝑒^𝑛₂−𝑘²𝑒^𝑛₂=0,inΩ_ext,𝑒^𝑛₂=𝑣^𝑛,onΓ₁,lim_𝑟→∞𝑟𝜕𝑒^𝑛₂𝜕𝑟−i𝑘𝑒^𝑛₂=0,−Δ𝑒^𝑛₁−𝑘²𝑒^𝑛₁=0,onΩ_int,𝑒^𝑛₁=0,onΓ_𝑅,𝜕𝑒^𝑛₁𝜕𝑛₁=−𝜕𝑒^𝑛₂𝜕𝑛₂,onΓ₁,𝑣^𝑛+1=𝜃_𝑛𝑒^𝑛_1Γ1+1−𝜃_𝑛𝑣^𝑛. (33)(34)(35) Meanwhile, 𝑒^𝑛₂=ℛ₂𝑣^𝑛,𝑒^𝑛₁=ℛ₁𝑒^𝑛₁. (36) Substituting the above equation to the third equation of (34), we can obtain 𝒮₁𝑒^𝑛₁=−𝒮₂𝑣^𝑛=−𝒦𝑣^𝑛. (37) Therefore, 𝒮₁𝑣^𝑛+1−𝑣^𝑛=𝒮₁𝜃_𝑛𝑒^𝑛_1Γ1−𝜃_𝑛𝑣^𝑛=−𝜃_𝑛𝒦𝑣^𝑛−𝜃_𝑛𝒮₁𝑣^𝑛=−𝜃_𝑛𝒮𝑣^𝑛; (38) thus, 𝒮₁𝜆^𝑛+1−𝜆^𝑛=𝜃_𝑛𝒮𝜆−𝜆^𝑛. (39) It comes that 𝜕𝒮𝜆=𝜕𝑛₁ℛ₁𝜆+𝜕𝜕𝑛₂ℛ₂𝜆=𝜕𝜕𝑛₁𝑢₁−𝒯₁𝑓+𝜕𝑢₂𝜕𝑛₂=−𝜕𝒯₁𝑓𝜕𝑛₁=𝜒. (40) Substituting the above equation in (39), we have 𝒮₁𝜆^𝑛+1−𝜆^𝑛=𝜃_𝑛𝜒−𝒮𝜆^𝑛. (41)

4. Convergence of the Algorithm

Considering a special case in which Γ₀={(𝑟,𝜃,𝜑)∣𝑟=𝑅₀,0⩽𝜃⩽2𝜋,0⩽𝜑⩽𝜋} and the interface Γ₁={(𝑟,𝜃,𝜑)∣𝑟=𝑅₁,0⩽𝜃⩽2𝜋,0⩽𝜑⩽𝜋}, following [10], we have the following.

Lemma 8. Suppose 𝑓_𝑛(𝑥)=𝑛+2𝑛+1𝑛+1𝑥2𝑛+1^𝑛+1,𝑥∈[0,1],𝑛∈𝑍⁺={0,1,2,…}, (42) and set 𝑈(𝑥)=sup_𝑛∈𝑍+𝑓_𝑛(𝑥),𝐿(𝑥)=inf_𝑛∈𝑍+𝑓_𝑛(𝑥); (43) then, we have 𝑈(𝑥)=121,0⩽𝑥⩽2,1𝑥,2𝐿(𝑥)=𝑓<𝑥⩽1,₀(𝑥),0⩽𝑥⩽𝑥₀,𝑓_𝑘+1(𝑥),𝑥_𝑘<𝑥⩽𝑥_𝑘+1, (44) where 𝑥_𝑘 is the only root of equation 𝑓_𝑘+1(𝑥)−𝑓_𝑘(𝑥)=0 in interval (0,1) and 𝑥_𝑘<𝑥_𝑘+1, 𝑘 is a nonnegative integer.

For 𝜆(𝜃,𝜑)∈𝐻^1/2(Γ₁), in the 𝐿²-sense, the following expansion holds 𝜆𝜃,𝜑=^∞𝑙_{𝑙=0𝑚=−𝑙}Λ^𝑚_𝑙𝑌^𝑚_𝑙𝜃,𝜑, (45) where 𝑃^𝑚_𝑙(𝑥) denotes the associated Legendre function of the first kind; 𝑌^𝑚_𝑙(𝜃,𝜑) denotes the complex conjugate of 𝑌^𝑚_𝑙. On the one hand, 𝒦^𝑐𝜆,𝜆=𝑅₁^∞𝑙_{𝑙=0𝑚=−𝑙}(1+𝑙)Λ^𝑚_𝑙² (46) satisfies 𝑅₁^∞𝑙_{𝑙=0𝑚=−𝑙}1+𝑙^21/2Λ^𝑚_𝑙^21/2<+∞;Λ^𝑚_𝑙=(−1)^|𝑚|Λ_𝑙^−𝑚, (47) on the other hand, set 𝐻^𝑁₁𝜆=^𝑁𝑙_{𝑙=0𝑚=−𝑙}𝑅₁^𝑙+1Λ^𝑚_𝑙𝑅₁^2𝑙+1−𝑅₀^2𝑙+1𝑟^𝑙−𝑅₀^2𝑙+1𝑟^𝑙+1𝑌^𝑚_𝑙,𝑅𝜃,𝜑₀<𝑟<𝑅₁. (48) Following [11] and (48), we have lim_𝑁→∞𝐻₁𝜆−𝐻^𝑁₁𝜆_𝐻1(Ω1)=0. (49) According to Green formula and orthogonality of 𝑌^𝑚_𝑙(𝜃,𝜑), it comes that 𝒮₁𝜆,𝜆=_Ω1∇𝐻₁𝜆⋅∇𝐻₁𝜆d𝐱=lim_{𝑀→∞,𝑁→∞Ω1}∇𝐻^𝑀₁𝜆⋅∇𝐻^𝑀₁𝜆d𝐱=lim_{𝑀→∞,𝑁→∞Γ1}𝜕𝐻^𝑀₁𝜆𝐻𝜕𝑟^𝑁1𝜆d𝑠=lim_{𝑀→∞,𝑁→∞}^{min{𝑀,𝑁}𝑙}_{𝑙=0𝑚=−𝑙}𝑙𝑅₁^2𝑙+1+(𝑙+1)𝑅₀^2𝑙+1𝑅₁^2𝑙+1−𝑅₀^2𝑙+1×𝑅₁Λ^𝑚_𝑙⋅Λ^𝑚_𝑙. (50)

Lemma 9. Suppose 𝒮₁𝜆,𝜆=𝑅₁^∞𝑙_{𝑙=0𝑚=−𝑙}𝑙𝑅₁^2𝑙+1+(𝑙+1)𝑅₀^2𝑙+1𝑅₁^2𝑙+1−𝑅₀^2𝑙+1Λ^𝑚_𝑙²;⟨𝒦𝜆,𝜆⟩=𝑅₁^∞𝑙_{𝑙=0𝑚=−𝑙}(𝑙+1)Λ^𝑚_𝑙², (51) then, 1𝑈𝑅₀/𝑅₁⩽⟨𝒮𝜆,𝜆⟩𝒮₁⩽1𝜆,𝜆𝐿𝑅₀/𝑅₁. (52)

Theorem 10. If relaxation factor 𝜃_𝑛 satisfies 0<𝜃_𝑛𝑅<2𝐿⁰𝑅₁, (53) then the D-N alternating algorithm (3)–(5) converges. Particularly, the optimal relaxation factor is 𝜃_𝑜𝑝𝑡=2𝑈𝑅₀/𝑅₁𝐿𝑅₀/𝑅₁𝑈𝑅₀/𝑅₁+𝐿𝑅₀/𝑅₁, (54) while the optimal compression ratio of the iteration is 𝛿_𝑜𝑝𝑡=𝑈𝑅₀/𝑅₁−𝐿𝑅₀/𝑅₁𝑈𝑅₀/𝑅₁+𝐿𝑅₀/𝑅₁. (55)

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 𝑅₁ and the geometry of Ω^𝑐 in some way, which is confirmed by our numerical results.

5. The Weak Form and Discretization

The weak form of problem (4) is Find 𝑢^𝑛₁∈𝐻¹_𝑔(Ω_int), such that 𝐷_Ωint𝑢^𝑛₁,𝑣=_Ωint𝑓𝑣d𝐱−_Γ1𝑣𝜕𝑢^𝑛₂𝜕𝑛₂d𝑠,∀𝑣∈𝐻¹₀Ω_int, (56) where 𝐷_Ωint𝑢^𝑛₁,𝑣=_Ωint∇𝑢^𝑛₁⋅∇𝑣+𝑘²𝑢^𝑛₁𝑣𝐻d𝐱,_1𝑔Ω_int=𝑣∣𝑣∈𝐻¹Ω_int,𝑣|_Γ𝑅,𝐻=𝑔₁₀Ω_int=𝑣∣𝑣∈𝐻¹Ω_int,𝑣|_Γ0=0. (57)

Obviously, 𝑣|_Γ1∈𝐻^1/2(Γ₁). Following (6), the weak form (56) can be rewritten as Find 𝑢^𝑛₁∈𝐻¹_𝑔(Ω_int), such that 𝐷_Ωint𝑢^𝑛₁,𝑣=_Ωint𝑓𝑣d𝐱−𝐷^𝑐𝜆^𝑛,𝑣,∀𝑣∈𝐻¹₀Ω_int,𝑖=1,2, (58) where 𝐷^𝑐𝜆^𝑛,𝑣=𝒦𝜆^𝑛,𝑣=−^∞𝑙_{𝑙=0𝑚=−𝑙}𝑘𝑅²₁ℎ^′_𝑙(𝑘𝑅)ℎ_𝑙()𝜆𝑘𝑅^𝑛𝑚𝑙𝑉^𝑚_𝑙,𝑉^𝑚_𝑙=₀^2𝜋𝜋₀𝑣𝜃^′,𝜑^′𝑌^𝑚_𝑙𝜃^′,𝜑^′sin𝜃^′d𝜃^′d𝜑^′,𝜆^𝑛𝑚_𝑙=₀^2𝜋𝜋₀𝜆^𝑛𝜃^′,𝜑^′𝑌^𝑚_𝑙𝜃^′,𝜑^′sin𝜃^′d𝜃^′d𝜑^′. (59) (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 𝑆_ℎ(Ω_int) denote the linear subspace of 𝐻¹₀(Ω_int) corresponding to this partition. Define 𝑀₀, 𝑀₁, and 𝑀_𝑖 to be the sets of all nodes belonging to Γ_𝑅, Γ₁, and Ω_int, respectively. Denote 𝑈^𝑛_𝐴 the approximate value and 𝐿_𝐴(𝐱) a basis function at node 𝐴. Obviously, 𝑆_ℎΩ_int=span𝐿_𝐴^𝑖(𝐱),𝐿_𝐴¹(𝐱)∣𝐴^𝑖∈𝑀_𝑖,𝐴¹∈𝑀₁. (60) The finite element approximation of 𝑢^𝑛₁ can be expressed as 𝑢^𝑛_1ℎ=_{𝐴𝑖∈𝑀𝑖}𝑈^𝑛_𝐴𝑖𝐿_𝐴^𝑖(𝐱)+_{𝐴1∈𝑀1}𝑈^𝑛_𝐴1𝐿_𝐴¹(𝐱)+_{𝐴0∈𝑀0}𝑔𝐱_𝐴⁰𝐿_𝐴⁰(𝐱). (61) By (61) and (5), we can obtain the following discrete D-N alternating algorithm: 𝑄¹¹𝑄^1𝑖𝑄¹⁰𝑄^𝑖1𝑄^𝑖𝑖𝑄^𝑖0𝑄⁰¹𝑄^0𝑖𝑄⁰⁰𝑈^𝑛₁𝑈^𝑛_𝑖𝑈^𝑛₀=−𝐵Λ^𝑛0𝐹⁰,Λ^𝑛+1=𝜃_𝑛𝑈^𝑛_𝐴1+1−𝜃_𝑛Λ^𝑛, (62)(63) where Λ^𝑛=𝜆^𝑛𝐴^1𝑇,𝑈^𝑛_𝑠=𝑈^𝑛_𝐴^𝑠𝑇,𝑄^𝑠𝑡=𝐷_Ωint𝐿_𝐴^𝑠,𝐿_𝐵^𝑡,𝑄𝑠,𝑡=0,1,𝑖,₀₁=𝑄¹⁰=0,𝐵=𝑏𝛾₀𝑁_𝐴¹,𝛾₀𝑁_𝐵¹,𝐹⁰=_Γ𝑅𝑔𝑁_𝐵⁰d𝑠^𝑇,𝐴^𝑠,𝐵^𝑠∈𝑀_𝑠,𝑠=0,1,𝑖. (64)

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: 𝑆_ℎ⁽¹⁾Λ^𝑛+1−Λ^𝑛=𝜃_𝑛𝐹¹−𝒮_ℎΛ^𝑛, (65) where 𝑆_ℎ⁽¹⁾=𝑄¹¹−𝑄^1𝑖𝑄^𝑖𝑖−𝑄^𝑖0𝑄⁰⁰⁻¹𝑄^0𝑖−1𝑄^𝑖1,𝑆_ℎ=𝑆_ℎ⁽¹⁾𝐹+𝐵,₁=𝑄^1𝑖𝑄^𝑖𝑖−𝑄^𝑖0𝑄⁰⁰⁻¹𝑄^0𝑖−1𝑄^𝑖0𝑄⁰⁰⁻¹𝐹⁰. (66)

Proof. Let [𝑈₁,𝑈_𝑖,𝑈₀]^𝑇 be the solution of the equation 𝑄¹¹+𝐵𝑄^1𝑖𝑄¹⁰𝑄^𝑖1𝑄^𝑖𝑖𝑄^𝑖0𝑄⁰¹𝑄^0𝑖𝑄⁰⁰𝑈₁𝑈_𝑖𝑈₀=00𝐹⁰. (67) Simplified to 𝑄¹¹−𝑄^1𝑖𝑄^𝑖𝑖−𝑄^𝑖0𝑄⁰⁰⁻¹𝑄^0𝑖−1𝑄^𝑖1+𝐵𝑈₁=𝑄^1𝑖𝑄^𝑖𝑖−𝑄^𝑖0𝑄⁰⁰⁻¹𝑄^0𝑖−1𝑄^𝑖0𝑄⁰⁰⁻¹𝐹⁰, (68) namely, 𝑆_ℎ𝑈₁=𝐹¹. (69) Then, (67) can be rewritten as 𝑄¹¹𝑄^1𝑖𝑄¹⁰𝑄^𝑖1𝑄^𝑖𝑖𝑄^𝑖0𝑄⁰¹𝑄^0𝑖𝑄⁰⁰𝑈₁𝑈_𝑖𝑈₀=−𝐵𝑈₁0𝐹⁰. (70) By (62) and (70), 𝑄¹¹𝑄^1𝑖𝑄¹⁰𝑄^𝑖1𝑄^𝑖𝑖𝑄^𝑖0𝑄⁰¹𝑄^0𝑖𝑄⁰⁰𝑈^𝑛₁−𝑈₁𝑈^𝑛_𝑖−𝑈_𝑖𝑈^𝑛₀−𝑈₀=−𝐵Λ^𝑛−𝑈₁00. (71) We have 𝑆_ℎ⁽¹⁾𝑈^𝑛₁−𝑈₁=𝐵𝑈₁−Λ^𝑛. (72) Therefore, 𝑆_ℎ⁽¹⁾Λ^𝑛+1−Λ^𝑛=𝜃_𝑛𝑆_ℎ⁽¹⁾𝑈^𝑛₁−Λ^𝑛=𝜃_𝑛𝑆_ℎ⁽¹⁾𝑈^𝑛₁−𝑈₁+𝑆_ℎ⁽¹⁾𝑈₁−Λ^𝑛=𝜃_𝑛𝐵+𝑆_ℎ⁽¹⁾𝑈₁−Λ^𝑛=𝜃_𝑛𝑆_ℎ𝑈₁−Λ^𝑛=𝜃_𝑛𝐹¹−𝑆_ℎΛ^𝑛. (73)

Theorem 13. Let 𝛼 be spectral radius of (𝑆_ℎ⁽¹⁾)⁻¹𝑆_ℎ, which is iterative matrix of preconditioned Richardson iteration. Then, there is a positive constant 𝛽, which is independent of finite element mesh parameter ℎ of subdomain Ω_int, such that 𝛼⩽𝛽.

Theorem 14. Eigenvalues of iterative matrix (𝑆_ℎ⁽¹⁾)⁻¹𝑆_ℎ are real and greater than 1.

Theorem 15. Put 𝜃_𝑛=𝜂(𝑛=0,1,…); then, there exists a constant 𝛿(0<𝛿<1), which is independent of finite element mesh parameter ℎ 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 ℎ of subdomain Ω_int.

Proof. By (73), we have 𝑈₁−Λ^𝑛+1=𝐼−𝜂𝑆_ℎ⁽¹⁾⁻¹𝑆_ℎ𝑈₁−Λ^𝑛=𝐼−𝜂𝑆_ℎ⁽¹⁾⁻¹𝑆_ℎ^𝑛+1𝑈₁−Λ⁰; (74) it comes that 𝑈₁−Λ^𝑛+1₂⩽𝛿^𝑛+1𝑈₁−Λ⁰₂, (75) where 𝛿=𝐼−𝜂(𝑆_ℎ⁽¹⁾)⁻¹𝑆_ℎ2. (76) Following Theorem 13 and Theorem 14, there exists a constant 𝛿(0<𝛿<1), which is independent of ℎ. For 0<𝜂<𝜎, spectral radius of 𝐼−𝜂(𝑆_ℎ⁽¹⁾)⁻¹𝑆_ℎ is less than 1, and spectral norm 𝛿<1; therefore, lim_𝑛→+∞𝑈₁−Λ^𝑛+1₂=0. (77) 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 ℎ 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 Section 2. 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 ℎ.

Suppose that Ω^𝑐={(𝑥,𝑦,𝑧)∣|𝑥|>1,|𝑦|>1,|𝑧|>1} is the exterior domain of cube [−1,1]×[−1,1]×[−1,1]. The artificial boundary is Γ₁={(𝑟,𝜃,𝜑)∣𝑟=𝑅₁,𝑅₁>√3}, 𝑟=𝑥²+𝑦²+𝑧². Split the bounded domain Ω_int into finite element mesh as follows. Firstly, make every inner boundary edge into 𝑁₁ equivalent parts and connect nodes by using segments which are parallel to the corresponding coordinate axis; then, the grid located in the inner boundary Γ₁ 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 𝑁₂ equivalent parts. At last, we get an eight-node trilinear isoparametric finite elements of Ω_int.

Denote 𝑒 the maximum node-error on Ω_𝑖: 𝑒(𝑛)=sup_{𝑃𝑖∈Ω𝑖}𝑢𝑃_𝑖−𝑢^𝑛_1ℎ𝑃_𝑖; (78)𝑒_ℎ denotes the maximum node-error of the abjacent two-steps on nodes 𝑒_ℎ(𝑛)=sup_{𝑃𝑖∈Ω𝑖}𝑢^𝑛+1_1ℎ𝑃_𝑖−𝑢^𝑛_1ℎ𝑃_𝑖; (79)𝑞_ℎ is the convergence rate 𝑞_ℎ(𝑛)=𝑒_ℎ(𝑛−1)𝑒_ℎ(𝑛). (80)

We substitute ∑^𝑀_𝑙=0 for ∑^∞_𝑙=0 in the computing of the entries of 𝐾_ℎ. Denote by 𝑚 the total number of nodes on Ω_int. By computing, the results are as follows.

Example 16. In the problem (1), we take 𝑔=3𝑘4𝜋²i−2i−2𝑘𝑘²𝑒^i𝑘cos𝜃; (81) the exact solution is 𝑢(𝑟,𝜃)=−34𝜋𝑘𝑟+i𝑘²𝑟²𝑒^i𝑘𝑟cos𝜃. (82) Let 𝑘=2.0, 𝑁₁=𝑁₂=2, 𝑚=78, 𝑀=4, the numerical results are shown in Tables 1 and 2.
Tables 1 and 2 show that the D-N alternating algorithm (62)-(63) is convergent. Table 1 shows that the speed of convergence is faster while 𝜃_𝑛=0.5. Conversely, Table 2 shows the different phenomenon that the speed of convergence is faster while 𝜃_𝑛=0.3. It proves Theorem 10 and shows relaxation factor 𝜃_𝑛 is dependent on 𝑅₁. But it is valuable to know that 𝜃_𝑛=0.5 unsatisfied the requirement of Theorem 10  (2𝐿(𝑅₀/𝑅₁)=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 Ω_𝑖 for 𝑅₁=2.0.

Table 2 

Maximum node-error on Ω_𝑖 for 𝑅₁=6.0.

Example 17. In the problem (1), we take 𝑔=15𝑘32𝜋³+4𝑘²i−9𝑘−9i𝑘³𝑒^i(𝑘+𝜙)sin2𝜃; (83) the exact solution is 𝑢(𝑟,𝜃)=−1532𝜋3i𝑘³𝑟³+3𝑘²𝑟²−i𝑒𝑘𝑟^{i(𝑘𝑟+𝜑)}sin2𝜃. (84) Let 𝑘=2.0, 𝑅₁=6.0, 𝑁₁=𝑁₂=2, 𝑚=78, 𝑀=4, the numerical results are shown in Table 3.
Table 3 shows that in the scope of relaxation factor, the bigger the relaxation factor 𝜃_𝑛 is, the faster the convergence rate is.
While we put 𝜃_𝑛=0.4 and 𝑅₁=6.0 and denote Mesh I: 𝑁₁=𝑁₂=1, 𝑚=16, 𝑀=4; Mesh II: 𝑁₁=𝑁₂=2, 𝑚=78, 𝑀=4; Mesh III: 𝑁₁=𝑁₂=4, 𝑚=490, 𝑀=4. The numerical results are shown in Table 4.
Table 4 shows that the convergence rate of the discrete D-N alternating algorithm is independent of mesh parameter ℎ. 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 𝑂(ℎ) order.

Table 3 

The relation between convergence rate and relaxation factor.

Table 4 

The relation between convergence rate and mesh.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The research is supported by the National Natural Science Foundantion of China, contact/Grant no. 11371198, Jiangsu Provincial Natural Science Foundation of China, contact/Grant no. BK20141008, and the Jiangsu Provincial Key Laboratory for Numerical Simulation of Large Scale Complex Systems contract/Grant no. 201305.