Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 418426, 8 pages

http://dx.doi.org/10.1155/2014/418426

Research Article

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

Qing Chen,¹ Baoqing Liu,² and Qikui Du¹

¹School of Mathematical Sciences, Jiangsu Provincial Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing 210023, China

²School of Applied Mathematics, Jiangsu Provincial Key Laboratory for NSLSCS, Nanjing University of Finance and Economics, Nanjing 210023, China

Received 6 May 2014; Revised 3 July 2014; Accepted 3 July 2014; Published 20 July 2014

Academic Editor: Xiao-Qiao He

Copyright © 2014 Qing Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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:

where is a bounded domain in , , and . denotes the wave number related to the wavelength of the incident wave through . 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:

where , . Take a sphere which surrounds , such that . Then the artificial boundary divides the original unbounded region into two subregions and , where is the outside subregion and is the bounded subregion surrounded by and , and are nonoverlaping domains. Then the original problem is turned into two subproblems over subdomains and . For exterior boundary value problem (1), we construct the following D-N alternating algorithm:

(i)Step 1. Select the initial , .(ii)Step 2. Solve the Dirichlet problem in :

(iii)Step 3. Solve the Neumann problem in :

(iv)Step 4. Input , and let

(v)Step 5. Let , and go to Step 2,

where and are the th approximate solutions in and , 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 . 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]:

where is the natural integral operator of Helmholtz equation in :

is the first-kind Hankel function of order and

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 , , ; it is easy to know , and and are solutions of the following equations, respectively:

By Green formula, we have

According to superimposition principle, the solutions of (9) and (10) can be denoted

where satisfy, respectively,

It is easy to verify that

Extend elements of and by zeros to and , respectively; therefore satisfies

while satisfies

where is the Possion integral operator of Helmholtz equation in the exterior spherical domain . Define

is , which is the natural integral operator of Helmholtz equation in the exterior spherical domain . Substituting (12) in (11), we have

Namely,

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

where is the Poincaré-Steklov operator [8] on . Note that , ; it comes that

Then,

Theorem 1. , is symmetric positive operator.

Proof. We note that ; therefore, for any . By Green formula, we know

where

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 ; that is, there exists a positive constant such that

where .

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,

is equivalent to finding , such that

Let

It is easy to verify that

Now we introduce as preprocessor and preconditioned Richardson iterative method. Given , the procedure of iterative

where

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

Lemma 5. For any , one has .

Theorem 6. If one puts and , preconditioned Richardson iterative algorithm (30) is convergent. Particularly, optimal relaxation factor is . 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

Meanwhile,

Substituting the above equation to the third equation of (34), we can obtain

Therefore,

thus,

It comes that

Substituting the above equation in (39), we have

4. Convergence of the Algorithm

Considering a special case in which and the interface , following [10], we have the following.

Lemma 8. Suppose

and set

then, we have

where is the only root of equation in interval and , is a nonnegative integer.

For , in the -sense, the following expansion holds

where denotes the associated Legendre function of the first kind; denotes the complex conjugate of . On the one hand,

satisfies

on the other hand, set

Following [11] and (48), we have

According to Green formula and orthogonality of , it comes that

Lemma 9. Suppose

then,

Theorem 10. If relaxation factor satisfies

then the D-N alternating algorithm (3)–(5) converges. Particularly, the optimal relaxation factor is

while the optimal compression ratio of the iteration is

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 , such that

where

Obviously, . Following (6), the weak form (56) can be rewritten as Find , such that

where

(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 into a finite number of element domains. Let denote the linear subspace of corresponding to this partition. Define , , and to be the sets of all nodes belonging to , , and , respectively. Denote the approximate value and a basis function at node . Obviously,

The finite element approximation of can be expressed as

By (61) and (5), we can obtain the following discrete D-N alternating algorithm:

where

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:

where

Proof. Let be the solution of the equation

Simplified to

namely,

Then, (67) can be rewritten as

By (62) and (70),

We have

Therefore,

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 , such that .

Theorem 14. Eigenvalues of iterative matrix are real and greater than .

Theorem 15. Put ; then, there exists a constant , which is independent of finite element mesh parameter of subdomain . For , 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 .

Proof. By (73), we have

it comes that

where

Following Theorem 13 and Theorem 14, there exists a constant , which is independent of . For , spectral radius of is less than , and spectral norm ; therefore,

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 .

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 is the exterior domain of cube . The artificial boundary is , . Split the bounded domain 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 .

Denote the maximum node-error on :

denotes the maximum node-error of the abjacent two-steps on nodes

is the convergence rate

We substitute for in the computing of the entries of . Denote by the total number of nodes on . By computing, the results are as follows.

Example 16. In the problem (1), we take

the exact solution is

Let , , , , 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 . Conversely, Table 2 shows the different phenomenon that the speed of convergence is faster while . It proves Theorem 10 and shows relaxation factor is dependent on . But it is valuable to know that unsatisfied the requirement of Theorem 10  , 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 .

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.06529

Table 2: Maximum node-error on for .

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.63435

Example 17. In the problem (1), we take

the exact solution is

Let , , , , , 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 and and denote Mesh I: , , ; Mesh II: , , ; Mesh III: , , . 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.

Iterative times

0
1
2
3
4
12

0.2

0.13381
0.07449
0.05900
0.05302
0.05072
0.04928

0.02358
0.05223
0.01157
0.00256
0.00057

4.51432
4.51329
4.51206
4.51101

0.3

0.13381
0.54416
0.04968
0.04330
0.04927
0.04927

0
0.50323
0.0806
0.01294
0.00209
0.00031

0
0
6.2432
6.2301
6.1952
6.1809

0.4

0.13381
0.03434
0.05268
0.04849
0.04945
0.04927

0
0.01512
0.00176
0.00021
0.00002
0

0
0
8.5961
8.5532
8.3201
7.9326

Table 4: The relation between convergence rate and mesh.

Mesh
Iterative times

0
1
2
3
4
12

I

0.16651
0.12025
0.11342
0.11354
0.11339
0.11327

0.11547
0.10300
0.08651
0.07381
0.06525

1.12105
1.19058
1.17205
1.13114

II

0.13381
0.03434
0.05268
0.04849
0.04945
0.04927

0.02834
0.01814
0.01162
0.00745
0.00478

1.56223
1.56134
1.56001
1.56018

III

0.08874
0.01239
0.01753
0.01623
0.01648
0.01644

0.01056
0.00449
0.00191
0.00083
0.00036

2.35073
2.35327
2.30124
2.29801

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.

References

D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, The Netherlands, 1992.

I. Harari and T. J. R. Hughes, “Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains,” Computer Methods in Applied Mechanics and Engineering, vol. 97, no. 1, pp. 103–124, 1992.

D. Givoli and J. B. Keller, “A finite element method for large domains,” Computer Methods in Applied Mechanics and Engineering, vol. 76, no. 1, pp. 41–66, 1989.

D. H. Yu, “Discretization of non-overlapping domain decomposition method for un-bounded domains and its covergence,” Mathematica Numerica Sinica, vol. 18, no. 3, pp. 328–336, 1996.

D. H. Yu and Z. P. Jia, “The nonoverlapping domain decomposition method based on natural boundary reduction for a 2-dimensional exterior Helmholtz problem,” Mathematica Numerica Sinica. Jisuan Shuxue, vol. 22, no. 2, pp. 227–240, 2000.

J. M. Wu and D. H. Yu, “The natural integral equations of 3-D exterior Helmholtz problem and its numerical solution,” Chinese Journal of Computational Physics, vol. 16, no. 5, pp. 449–456, 1999.

Z. P. Jia, J. M. Wu, and D. H. Yu, “A coupled natural boundary element and finite element method for solving a 3-dimensional exterior Helmholtz problem,” Mathematica Numerica Sinica, vol. 23, no. 3, pp. 357–368, 2001.

J. C. Xu and S. Zhang, “Preconditioning the Poincaré-Steklov operator by using Green's function,” Mathematics of Computation, vol. 66, no. 217, pp. 125–138, 1997.

T. Läu, J. M. Shi, and Z. B. Lin, Domain Decomposition Method. A New Technique for Numereical Solution of Partial Differential Equations, Science Press, Beijing, China, 1992.

D. Yu and J. M. Wu, “A nonoverlapping domain decomposition method for exterior 3-D problem,” Journal of Computational Mathematics, vol. 19, no. 1, pp. 77–86, 2001.

J. M. Wu and D. H. Yu, “The natural integral equations of 3-D harmonic problems and their numerical solutions,” Mathematica Numerica Sinica, vol. 20, no. 4, pp. 419–430, 1998.