A D-N Alternating Algorithm for Solving 3D Exterior Helmholtz Problems
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.
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  and Helmholtz equation . Wu and Yu  studied the natural integral equations of 3D Helmholtz problems. Jia et al.  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 : 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  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 , 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 .
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 , 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  and (48), we have According to Green formula and orthogonality of , it comes that
Lemma 9. Suppose then,
Remark 11. We have analyzed the convergence for exterior spherical domain. Following , 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 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.
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
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.
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.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar