A D-N alternating algorithm for exterior 3-D problem with ellipsoidal artificial boundary

1.   Introduction

In this paper, we consider the following exterior Dirichlet problem:

where $\Omega_0\subset \mathbb{R}^3$ be a bounded Lipschitz domain, and $\Gamma_0 = \partial\Omega_0$ be its boundary, $\Omega = R^3\backslash\overline{\Omega}_0$. $\Delta$ is the Laplace operator. Assume that the given function $g$ satisfies $g\in H^{\frac{1}{2}}(\Gamma_0)$ and $f(\textbf{x})\in L^2(\Omega)$ has the compact supported set.

It is well known that partial differential equation boundary value problems with bounded or unbounded domains widely arise in many fields of scientific and engineering. For partial differential equation boundary value problems with bounded domain, the commonly finite element method (FEM) ^[1] and finite difference method (FDM) ^[2] were used to solve these problems. Meanwhile, for unbounded domains partial differential equation, boundary element method (BEM) ^[3,4,5], natural boundary element method (NBEM) ^[6,7], artificial boundary method ^{[8,9,10,11,12,13]}, and the coupling of NBEM (or BEM) and FEM ^[14,15] were all efficient numerical methods. Both the methods on bounded domains and the methods on unbounded domains, they all have their own advantages and disadvantages.

It should be pointed out that the coupling method of FEM and BEM proposed by Feng and Yu ^[6,7] is fully compatible with the FEM. It is easy to be implemented on the computer and can be coupled with FEM naturally and directly by natural boundary reduction (NBR). The coupling of FEM and NBEM was first applied to solve two dimensional (2-D) unbounded domain of the elliptic problems ^[16,17]. A circle ^[11] or an ellipse ^[12,16] was usually selected as the artificial boundary for exterior 2-D problems. For exterior 3-D problems, a sphere was usually selected as the artificial boundary ^[12,18,19]. For some special shapes, such as cigar-shaped, rotating ellipsoid boundary, some special artificial boundary methods ^{[20,21,22,23]} are given. Since it leads to a smaller computational domain and does not increase the computational cost of the stiffness matrix from boundary reduction, it shows that these methods are very efficient. In the last decade, mesh-less methods which involve the radial basis functions (RBF) are also used to solve higher dimensions exterior problems. The motivations is that such functions or kernels are easily applicable and implementable in higher order partial differential equations, independent from geometry, etc^[24,25,26].

In this paper, based on FEM and NBEM, a Dirichlet-Neumann (D-N) alternating algorithm has been devised to solve exterior 3-D Poisson problem in an infinite region with an ellipsoidal boundary. The standard procedure of the method is described as follows: Firstly, let $\Omega$ be an unbounded domain with boundary $\Gamma_0$. Then we introduce an artificial boundary $\Gamma_1$, where $\Gamma_0$ is surrounded by $\Gamma_1$. We denote the bounded domain between $\Gamma_0$ and $\Gamma_1$ as $\Omega_1$, and $\Omega_2$ is the unbounded domain with boundary $\Gamma_1$. Secondly, in $\Omega_2$, we use the NBEM to solve this problem. Furthermore, we use FEM to solve this problem in $\Omega_1$. Finally, on the interface between $\Omega_1$ and $\Omega_2$, a sequence of boundary conditions is generated iteratively until convergence to the solution of the original problem.

The outline of the paper is as follows: In Section 2, the original domain $\Omega$ is divide two non-overlapping regions by an ellipsoidal artificial boundary $\Gamma_1$ and the corresponding D-N alternating algorithm of the problem is constructed. In Section 3 and Section 4, the convergence analysis of the algorithm and the error estimates were given, respectively. In Section 5, two numerical examples are carried out to demonstrate the effectiveness and accuracy of this method. Finally, we state the conclusions in Section 6.

2.   A D-N alternating algorithm based on NBR

Let $\Gamma_1 = \{(x, y, z)\; |\; \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2} = 1, a > b > c > 0\}$ denote an ellipsoid. Then an unbounded domain outside the boundary $\Gamma_1$ is $\Omega_2$.

Let $\Gamma_0$ be a surface inside the $\Gamma_1$, and dist$(\Gamma_1, \Gamma_0) > 0$. Then for problem (1.1), $\Omega$ is divided into two non-overlapping subdomains $\Omega_1$ and $\Omega_2$.

Let $\Omega_1$ be the bounded domain among $\Gamma_0$ and $\Gamma_1$, and $\Omega_2$ be the unbounded domain outside $\Gamma_1$. Then the original problem (1.1) is decomposed into two subproblems over subdomains $\Omega_1$ and $\Omega_2$.

Firstly, based on ellipsoidal system of coordinates $(\lambda_1, \lambda_2, \lambda_3)$ which is defined in ^[27] and $\Gamma_1$ coincides with the ellipsoid $\lambda_1 = a$, the Cartesian coordinates $(x, y, z)$ related to the ellipsoidal coordinates $(\lambda_1, \lambda_2, \lambda_3)$ can be expressed as follows:

where $h^2 = a^2-b^2$ and $k^2 = a^2-c^2$.

Secondly, the Laplace's equation $\Delta u = 0$ in the following problem:

can be expressed in terms by the ellipsoidal coordinates $(\lambda_1, \lambda_2, \lambda_3)$. In the orthogonal set of coordinates $(\lambda_1, \lambda_2, \lambda_3)$, the Laplace's equation is expressed as in ^[27]

In exterior $\Omega_2$, the normal solution of the Laplace's equation ^[27] is

As we know, the functions $E_1$–$E_3$ must satisfy Lam$\acute{\rm e}$'s equation

where $p$ and $q$ are called the separation parameters.

Then, the Laplace equation's analytic solution over the domain $\Omega_2$ is

where

and ${\rm d}s$ is the surface element about the Cartesian coordinates on $\Gamma_1$. Let $E_n^p(\lambda)$ denote $n$ order Lam$\acute{\rm e}$ functions of the first kind with the eigenvalue $p$. As suggested by the Lam$\acute{\rm e}$ function of the first kind, there exist the Lam$\acute{\rm e}$ functions of the second kind. The Lam$\acute{\rm e}$'s functions of the second kind are then defined by

such that $F_n^p(\lambda_1)\rightarrow 0$ when $\lambda_1\rightarrow \infty$ and they have the same eigenvalue $p$.

The expression (2.3) is called the Poisson integral formula. From (2.3) and the Lemma 1 of Huang ^[21], an expression of the normal derivative $\frac{\partial u}{\partial\textbf{n}}$ can be obtained. If we set $\lambda_1 = a$, the exact artificial boundary condition is

Here, the expression (2.4) is also called the natural integral equation. So we can obtain the solution of problem (2.2) directly from (2.3).

Based on (2.3) and (2.4), we construct a D-N alternating algorithm as follows:

Step 1. Put $k = 0$, and pick an initial value $\lambda^0\in H^{\frac{1}{2}}(\Gamma_1)$.

Step 2. In exterior domain $\Omega_2$, solve a Dirichlet boundary value problem:

Step 3. In interior domain $\Omega_1$, solve a mixed boundary value problem:

Step 4. On $\Gamma_1$, update the boundary value by

Step 5. Put $k = k+1$, go to Step 2.

For the D-N alternating algorithm, we need know the relaxation factor $\theta_k$, which is called a suitable real number. The algorithm's main work is on the Step 3, where a mixed boundary value problem in bounded domain $\Omega_1$ is solved by standard finite element method. In order to solve the Dirichlet problem with exterior domain $\Omega_2$ in the Step 2, we can use the normal derivative of the solution of the problem (2.5) on $\Gamma_1$. Similarly, $\frac{\partial u_2^{k}}{\partial\textbf{n}}$ can be found from $\lambda^k$ by using the natural integral equation (2.4).

3.   Analysis of convergence

Let

Then the following coupled variational problem is equivalent to the problem (1.1).

where $D_1(u, v) = \int_{\Omega_1}\nabla u\cdot\nabla v{\rm d}x{\rm d}y{\rm d}z$. From ^[20], we know

where

Let $S_h(\Omega_1)\subset V_0(\Omega_1)$ be an linear finite element space of $V_0(\Omega_1)$, and subdivide $\Omega_1$ into hexahedrons or tetrahedrons. Then, the approximation variational problem of (3.1) can be written as:

From the problem (3.3), we can obtain the algebraic equations as follow:

where $U$ and $V$ are vectors whose components are function values at interior nodes of $\Omega_1$ and nodes on $\Gamma_1$, respectively.

From FEM in $\Omega_1$, we can get the matrix

which is a stiffness matrix. From the NBEM on $\Gamma_1$, we can obtain the $K_h$. The expression (3.4) can also be rewritten as follow:

Then, the iterative algorithm can be given as follow:

with

Since $A$ is a symmetric and positive definite matrix, we know that $A_{11}^{-1}$ is nonsingular. From the expression (3.4), we know

Let

then the expression (3.8) can also be rewritten $S_h\cdot V = \overline{b}_2$, where $S_h$ is the discrete analogue of the Steklov-Poincar$\acute{\rm e}$ operator on $\Gamma_1$. Therefore, preconditioned Richardson iteration is constructed as follow:

Theorem 3.1. ^[23] The discrete D-N alternating algorithms (3.6), (3.7) and the iteration (3.9) are equivalent.

Theorem 3.2. ^[23] The condition number of the iterative matrix $(S_h^{(1)})^{-1}S_h$ of the discrete D-N alternating algorithm is independent of the finite element mesh size $h$.

Theorem 3.3. ^[23] If $0 < \min\theta_k\leq\max\theta_k < 1$, then the D-N alternating methods (3.6) and (3.7) are convergent, and the convergence rate is independent of the finite element mesh size $h$.

4.   Error estimates

Let $V_h(\Omega_1)$ be a piecewise linear finite element space and $\mathcal{N}_h$ be the node set in $\Omega_1\cup\Omega_0\cup\Gamma_1$. We denote

The discrete variational problems of the variational problems (3.1) and (3.3) are the followings respectively:

and

Similarly, the following error estimation can be obtained from ^[21].

Theorem 4.1. Suppose that $u\in H^2(\Omega_1)\cap V_g(\Omega_1)$ and $u^{Nh}\in V_g(\Omega_1)$ are the solutions of problems (3.1) and (4.2) respectively, and $u|_{\Gamma_1}\in H^{\frac{3}{2}}(\Gamma_1)$. Then there exists a positive constant $C$ independent of $h$ and $N$ such that

Theorem 4.2. Suppose that $u\in H^2(\Omega_1)\cap V_g(\Omega_1)$ and $u^{Nh}\in V_g(\Omega_1)$ are the solutions of problems (3.1) and (4.2) respectively, and $u|_{\Gamma_1}\in H^{\frac{3}{2}}(\Gamma_1)$. Then there exists a positive constant $C$ independent of $h$ and $N$ such that

where $\alpha$ is related to the smoothness convexity of the domain $\Omega_1$ and the differentiability of the solution of boundary value problem.

5.   Numerical experiments

In this section, two numerical examples are used to show the effectiveness of the D-N alternating algorithm, where the exact solutions of all two examples are known. Let

be the maximal error of all node functions on $\overline{\Omega}_1$,

be the maximum node-error of the adjacent two-steps on nodes and

be the approximation of the convergence rate. We use three meshes: Mesh I, Mesh II and Mesh III to express our computations. Mesh I is a coarse mesh by a lot of small tetrahedrons. Mesh II is refined from Mesh I in such a way what every tetrahedron of Mesh I is divided into eight "equal" tetrahedron. Mesh III is the refined mesh of Mesh II in the same way.

Example 5.1. Let $f(\textbf{x}) = 0$, $\Gamma_0 = \{(x, y, z): \frac{x^2}{6.25}+\frac{y^2}{4.0}+\frac{z^2}{2.25} = 1\}$, the ellipsoidal artificial boundary $\Gamma_1 = \{(x, y, z): \frac{x^2}{6.25}+\frac{y^2}{4.0}+\frac{z^2}{2.25} = a^2, a > 1\}$ and the exact solution of problem (1.1) is

Taking $g(\textbf{x}) = u(\textbf{x})|_{\Gamma_0}$. The correspondent results are shown in Table 1, Table 2 and Figure 1.

Example 5.2. Let $f(\textbf{x}) = 0$, $\Omega = \{(x, y, z):|x|\geq 2.5, |y|\geq 2.0, |z|\geq 1.5\}$ and the ellipsoidal artificial boundary $\Gamma_1 = \{(x, y, z): \frac{x^2}{6.25}+\frac{y^2}{4.0}+\frac{z^2}{2.25} = a^2, a > \sqrt{3}\}$. The exact solution of problem (1.1) be

Taking $g(\textbf{x}) = u(\textbf{x})|_{\Gamma_0}$. The correspondent results is shown by Table 3, Table 4 and Figure 2.

In the numerical experiments, we use the $\Gamma_1 = \{(x, y, z): \frac{x^2}{6.25}+\frac{y^2}{4.0}+\frac{z^2}{2.25} = 4, \}$ as the numerical example 1 and example 2. The radius $a = 2$ just represents the shape of the $\Gamma_1$. The larger $a$ is, the larger domain of $\Gamma_1$ and $\Gamma_0$ enclosed will be. So the more computation time will be spent. On the contrary, The smaller $a$ is, the smaller area of $\Gamma_1$ and $\Gamma_0$ enclosed will be. So The less computation time is spent for the computer.

6.   Conclusions

In this paper, we have proposed a D-N alternating algorithm based on the NBR for exterior 3-D Poisson problem with ellipsoidal artificial boundary. The maximal error of the approximate solution decreases with the mesh refining from the numerical results Table 1 and Table 3. When we fixed convergence factor $\theta_k$, the convergence rate can be basically considered to maintain invariable in a same mesh. Table 2 and Table 4 have demonstrated that the convergence of the algorithm is the best as $\theta_k$ approaches to 0.5. Figure 1 and Figure 2 have shown that maximal errors in relative maximum norm will quickly dwindle until approach to stable states. The numerical results agree with the theoretic results.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 11401054, 61972055), the Hunan Provincial Natural Science Foundation of China (No. 2018JJ3520), and the general Project of Hunan provincial education department (No. 18C0220).

Conflict of interest

We declare that we have no conflict of interest.