- Email: [email protected]
A regularized multi-level technique for solving potential problems by the method of fundamental solutions
Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A regularized multi-level technique for solving potential problems by the method of fundamental solutions Csaba Gáspár Széchenyi István University, Egyetem tér 1, H-9026 Győr, Hungary
art ic l e i nf o
a b s t r a c t
Article history: Received 6 October 2013 Accepted 7 May 2014
The method of fundamental solutions is investigated in the case when the source points are located along the boundary of the domain of the original problem and coincide with the collocation points. The appearing singularities are eliminated by several techniques: by using approximate but continuous fundamental solutions (regularization) and via auxiliary subproblems to avoid the stronger singularities that appear in the normal derivatives of the fundamental solution (desingularization). Both monopole and dipole formulations are investigated. A special iterative solution algorithm is presented, which converts the original (mixed) problem to a sequence of pure Dirichlet and pure Neumann subproblems. The pure subproblems can be handled efficiently by using conjugate gradients. The efficiency is significantly increased by embedding the resulting method in a natural multi-level context. At the same time, the problem of the use of highly ill-conditioned matrices is also avoided. & 2014 Published by Elsevier Ltd.
Keywords: Meshless method Method of fundamental solutions Regularization Desingularization Multi-level method
1. Introduction Elliptic partial differential equations play an essential role in a lot of fields of application. Modelling stationary phenomena such as diffusion in fluids or in gases, heat transfer in machines e.g. in parts of traditional or hybrid cars, seepage through porous media lead to solving elliptic problems. The usual implicit time discretization techniques of time-dependent problems result in elliptic problems as well, at every time step. To handle elliptic problems in a meshless way, the popular Boundary Element Method is not suitable, since it requires a boundary mesh structure. The strength of the meshless methods is to circumvent the generation of both domain and boundary mesh. So far, a number of boundary meshless methods have been developed e. g. the boundary knot method [3,4] which uses nonsingular general solutions, or the method of fundamental solutions (MFS, see e.g. [1]), which is based on the fundamental solution of the applied partial differential operator, i.e. on solutions with singularities. In this paper, we restrict ourselves to the MFS applied to homogeneous problem. (For non-homogeneous problems, the approach can be combined with the well-known principle of the Method of Particular Solutions). Consider a second-order elliptic homogeneous linear partial differential equation: Lu ¼ 0
in Ω
E-mail address: [email protected]
ð1Þ
defined in a sufficiently smooth domain Ω supplied with mixed boundary conditions: ∂u ujΓD ¼ u0 ; ¼ v0 ; ð2Þ ∂nΓN where Γ≔∂Ω denotes the boundary of Ω, which has a disjoint decomposition into a Dirichlet part Γ D and a Neumann part Γ N . In its traditional form, the MFS produces an approximate solution of the problem (1) and (2) in the following form: N
uN ðxÞ ¼ ∑ αj Φðx x~ j Þ;
ð3Þ
j¼1
where Φ is a fundamental solution of the operator L, i.e. ΔΦ ¼ δ (here δ denotes the Dirac distribution concentrated at the origin). The predefined points x~ 1 ; …; x~ N (the source points) are located outside of the domain Ω. The a priori unknown coefficients α1 ; …; αN can be computed by enforcing the boundary conditions: N
∑ αj Φðxk x~ j Þ ¼ u0 ðxk Þ
j¼1 N
∑ αj
j¼1
∂Φ ðx x~ j Þ ¼ v0 ðxk Þ ∂nk k
ðxk A Γ D Þ; ð4Þ ðxk A Γ N Þ;
where x1 ; …; xN A Γ are predefined boundary collocation points. Then the function uN exactly satisfies Eq. (1), and exhibits singularities at the source points. Instead of the formulation (3) (called monopole formulation hereafter), it is often more advantageous to use the dipole
http://dx.doi.org/10.1016/j.enganabound.2014.05.002 0955-7997/& 2014 Published by Elsevier Ltd.
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
formulation, where the approximate solution of (1) and (2) is sought in the following form: ð5Þ
j¼1
∂Φ ðx x~ j Þ ¼ u0 ðxk Þ ∂nj k
∂2 Φ ∑ αj ðxk x~ j Þ ¼ v0 ðxk Þ j ¼ 1 ∂nk ∂nj
ð6Þ ðxk A Γ N Þ:
Note that the monopole and the dipole formulations can be considered as meshless discretizations of the indirect BEM based on single layer and double layer potentials, respectively. A common disadvantage of the traditional forms (3) and (5) is the use of external source points, the location of which can be hardly automatized. Though the MFS has excellent accuracy properties (see e.g. [10] and references therein), the systems (4) and (6) are highly illconditioned in general, which is a severe drawback of the method. This is the case especially when the source points are located far from the boundary. On the other hand, if they are close to the boundary, the systems (4) and (6) become much better conditioned, however, the accuracy goes wrong due to the appearance of numerical singularities at the boundary collocation points. A usual technique is to define the source and the boundary collocation points to coincide. Special techniques are required to avoid the problem of singularity (regularization and desingularization, see e.g. [14,7,12,8]). In this paper, we investigate some regularized versions of both the monopole and the dipole formulations. It turns out that, from a computational point of view, the dipole formulation is much more advantageous for handling pure Dirichlet problems, while for pure Neumann problems, the monopole formulation overperforms the dipole formulation. For mixed boundary conditions, a special iterative technique is proposed which converts the original mixed problem to a convergent sequence of pure Dirichlet and pure Neumann subproblems. This results in a computationally efficient method, which avoids also the problem of highly ill-conditioned linear systems. The efficiency can be increased further by embedding the method in a natural multi-level context.
In the rest of the paper, suppose that the source points and the boundary collocation points x1 ; …; xN A Γ coincide. Then the monopole formulation has the form uN ðxÞ ¼ ∑ αj Φðx xj Þ;
ð7Þ
j¼1
N
∑ αj Akj ¼ u0 ðxk Þ
j¼1
ðxk A Γ D Þ; ð8Þ
∑ αj Bkj ¼ v0 ðxk Þ
ðxk A Γ N Þ:
Similarly, the dipole formulation has the form N
uN ðxÞ ¼ ∑ αj j¼1
Here the entries of the matrices A, B, C, Q are defined as follows:
C kj ¼
∂Φ ðx xj Þ; ∂nj
∂Φ ðx xj Þ; ∂nj k
Bkj ¼
∂Φ ðx xj Þ; ∂nk k
Q kj ¼
∂2 Φ ðx xj Þ ∂nk ∂nj k
ð11Þ
for j ak. Due to the singularity of the fundamental solution at the origin, the diagonal entries cannot be computed by the above definition. For the proper definition of Akk, one should replace the fundamental solution Φ with an approximate fundamental solution Φ, which has no singularity at the origin. Such an approximate fundamental solution can be defined e.g. by truncation. In polar coordinates: 8 1 > > > if r Z < ΦðrÞ c ð12Þ ΦðrÞ≔ 1 1 > > if r o ; Φ > : c c provided that Φ is a radial function i.e. it depends only on r ¼ jjxjj, which is often the case. Here c denotes a carefully chosen scaling constant which should remain inversely proportional to the characteristic distance of the boundary collocation points, when N varies. Another regularization technique is to replace Φ with the fundamental solution of the singularly perturbed fourth-order operator LðI ð1=c2 ÞLÞ, where I denotes the identity operator and c is again a scaling constant, see [8] for details. Thus, the diagonal terms Akk can be computed without difficulty. Using the simplest truncation, Akk ¼ Φð0Þ, while for j a k, Akj ¼ Φðxk xj Þ. The proper definition of the diagonal terms Bkk is somewhat more difficult since the derivatives of Φ have stronger singularities at the origin than Φ itself. Let w be a smooth, easily computable particular solution of Eq. (1), then w can be approximated by the monopole formulation: N
wN ðxÞ≔ ∑ βj Φðx xj Þ; j¼1
N ∂wN ∂Φ ðxÞ ¼ ∑ βj ðx xj Þ: ∂n ∂n j¼1
Hence ∑ βj Bkj þ βk Bkk ¼
jak
where the coefficients α1 ; …; αN can be computed by solving the algebraic system:
j¼1
j¼1
ðxk A Γ N Þ:
from where the coefficients β1 ; …; βN can be computed by solving (1) supplied with a pure Dirichlet condition. Computing the normal derivative of the particular solution w, we have the following:
2. Regularization and desingularization
N
ðxk A Γ D Þ; ð10Þ
∑ αj Q kj ¼ v0 ðxk Þ
Akj ¼ Φðxk xj Þ;
ðxk A Γ D Þ;
N
N
j¼1 N
Again, Φ is the fundamental solution of the operator L, and the source points x~ 1 ; …; x~ N are located outside of Ω. Now the coefficients α1 ; …; αN can be computed by solving the linear system: N
N
∑ αj C kj ¼ u0 ðxk Þ
∂Φ uN ðxÞ ¼ ∑ αj ðx x~ j Þ: ∂n j j¼1 N
∑ αj
where the coefficients α1 ; …; αN solve the algebraic system:
ð9Þ
∂w ðx Þ ∂nk k
ðk ¼ 1; 2; …; NÞ:
Thus, the diagonal terms Bkk can be defined as ! 1 ∂w ∑ βj Bkj þ ðxk Þ : Bkk ≔ βk ∂nk jak This is the desingularization idea, see e.g. [14,11,5] for details. Note that in the simplest case of the Laplace equation w can often be chosen e.g. by w : 1; in the case of the modified 2D Helmholtz equation Δu λ2 u ¼ 0, a particular solution wðxÞ≔I 0 ðλ jjxjjÞ can be used, where I0 denotes the familiar modified Bessel function of the first kind, and so on. By a proper choice of the particular solution w, one can ensure that the coefficients βk do not vanish, so that the
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
definition of Bkk is correct (using possibly several particular solutions). The desingularization can be improved by redefining not only the diagonal terms but also the neighbouring entries [9]. In the dipole formulation, a straightforward desingularization is C kk ≔Bkk (while for ja k, C kj ¼ Bjk ). Finally, the diagonal terms Qkk can be defined by a desingularization procedure similar to the monopole formulation, by solving an auxiliary pure Dirichlet problem. See [8] for details. Remark. The definition of the diagonal entries C kk ≔Bkk is straightforward but not trivial, so that it needs some explanation. For simplicity, assume that Ω is a polygon with sides Γj (j ¼ 1; 2; …; N) and the boundary collocation point xj is the midpoint of Γj (more precisely let us approximate Ω in this way provided that N is large enough). Let us express the solution of the original problem in the form of a single layer potential: Z 1 uðξÞ ¼ ðlog ξ yÞ wðyÞ dΓ y ðξ A ΩÞ; 2π Γ with some continuous density function w. It is well known that the normal derivative of a single layer potential has a jump at the boundary i.e. Z ∂u 1 〈xk y; nk 〉 1 lim ðξÞ ¼ wðyÞ dΓ y wðxk Þ: ξ-xk ∂nk 2π Γ jjxk yjj2 2 (Here 〈; 〉 denotes the Euclidean scalar product.) The integral in the right-hand side can be expressed as a sum of integrals taken along the sides Γj, whence Z ∂u 1 〈xk y; nk 〉 ðξÞ ¼ ∑ wðyÞ dΓ y lim 2 ξ-xk ∂nk j a k2π Γ j jjxk yjj Z 1 〈xk y; nk 〉 1 þ wðyÞ dΓ y wðxk Þ: 2π Γ k jjxk yjj2 2 The last integral (taken along Γk) vanishes since xk y is orthogonal to the normal vector nk. Approximating the preceding integrals by the midpoint rule, we obtain the following discrete form of the equations: 1 ∑ Bkj hj wj wk ¼ vk ; 2 jak where wj ≔wðxj Þ, hj ≔jΓ j j (the length of Γj) and Bkj ¼
1 〈xk xj ; nk 〉 ∂Φ ¼ ðx xj Þ: 2π jjxk xj jj2 ∂nk k
Denoting by αj ≔hj wj , we have ∑ Bkj αj
jak
1 α ¼ vk ; 2hk k
whence Bkk ¼ 1=ð2hk Þ. On the other hand, expressing the solution of the original problem in the form of a double layer potential: Z 1 〈ξ y; ny 〉 uðξÞ ¼ wðyÞ dΓ y ðξ A ΩÞ; 2π Γ jjξ yjj2 (with a density function w which is different from the previous case). Using the well-known theorem of the jump of the double
3
layer potential at the boundary, we have Z 1 〈xk y; ny 〉 1 lim uðξÞ ¼ wðyÞ dΓ y wðxk Þ ¼ uk : ξ-xk 2π Γ jjxk yjj2 2 A procedure similar to the case of the single layer potential yields that these equations can be approximated by the discrete forms: ∑ C kj αj
jak
1 α ¼ uk ; 2hk k
where 1 〈xk xj ; nj 〉 C kj ≔ 2π jjxk xj jj2
ðj a kÞ
and the diagonal entries are C kk ¼ 1=ð2hk Þ. Thus, C kk ¼ Bkk , and for ja k: Bjk ¼
1 〈xj xk ; nj 〉 ¼ C kj : 2π jjxj xk jj2
Note that the above definitions of Bkk and Ckk mimic the jump of the normal derivative of the single layer potential at the boundary and the jump of the double layer potential, respectively. Both the monopole and the dipole formulations are suitable for solving mixed boundary problems using a regularized fundamental solution; the singularities of the diagonal entries of the systems (8) and (10) are avoided by regularization and desingularization. To illustrate the above techniques, consider the simplest 2D Laplace equation Δu ¼ 0
ð13Þ
in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution uðxÞ≔log ððxð1Þ þ2Þ2 þ ðxð2Þ þ 3Þ2 Þ; ð1Þ
ð14Þ
ð2Þ
where x , x denote the components of the vector variable x. The boundary was discretized by N collocation points in an equidistant way. The scaling parameter was set as follows: c≔N, which is a quasioptimal choice (more precisely, define c≔const:=h, where h denotes the characteristic distance of the boundary collocation points, and const: means a dimensionless constant which is independent of the discretization of the boundary). The problem was supplied with mixed boundary conditions; along the half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. Table 1 shows the relative L2-errors of the approximate solution computed on the boundary Γ and also the condition numbers of the corresponding algebraic systems. The method applied was the desingularized MFS based on monopoles (7) and dipoles (9). This seems an acceptable compromise: the accuracy is still acceptable and the condition numbers are still moderate; they are in fact much less than in the case of a traditional MFS with external source points located from the boundary at a distance which is comparable with the size of the domain. Remark. For Neumann and mixed problems, it is sufficient to compute the error of the approximate solution along the boundary only. Once the boundary values of the solution have been (approximately) computed, the solution inside the domain can be reconstructed by solving an additional pure Dirichlet problem.
Table 1 MFS with truncation. Mixed boundary conditions. Relative L2-errors on the boundary and condition numbers. N
8
16
32
64
128
256
512
Monopoles: relative L2-error (%) Condition number Dipoles: relative L2-error (%) Condition number
0.4541 8 0.4982 3
0.2009 21 0.1516 6
0.1058 53 0.0458 13
0.0552 127 0.0177 29
0.0282 294 0.0081 63
0.0142 670 0.0040 135
0.0071 1.5E þ3 0.0020 287
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
Table 2 MFS with truncation and desingularization. Pure Neumann boundary condition (using monopoles) and pure Dirichlet boundary condition (using dipoles). Relative L2-errors on the boundary and condition numbers. N
8
16
32
64
128
256
512
Monopoles: relative L2-error (%) Condition number Dipoles: relative L2-error (%) Condition number
0.5580 15 4.8602 3
0.3566 17 1.3709 3
0.2057 18 0.4715 3
0.1111 19 0.1869 4
0.0568 21 0.0821 4
0.0286 30 0.0385 4
0.0143 44 0.0186 4
This latter problem can be treated by a regularized MFS and needs no desingularization. The condition numbers dramatically decrease if the above methods are applied to pure problems. More precisely, the regularizeddesingularized MFS based on monopoles (7) is quite efficient for solving pure Neumann problems, while the dipole formulation (9) overperforms the regularized MFS when solving pure Dirichlet problems. This is in accordance with the indirect boundary element methods: pure Dirichlet (resp. Neumann) problems should be treated by double layer (resp. single layer) potentials. To illustrate this, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as earlier. Table 2 shows the relative L2-errors of the approximate solution computed on the boundary Γ and the condition numbers of the corresponding algebraic system (8), when a pure Neumann condition was prescribed. Here the monopole formulation was applied. Table 2 shows also the relative L2-errors of the normal derivative of the approximate solution computed on the boundary Γ and also the condition numbers of the corresponding algebraic system (10), when a pure Dirichlet condition was prescribed. Here the dipole formulation was applied. Since the algebraic system is always singular (due to the fact that the solution of the pure Neumann problem is unique up to an additive constant only), the first row of the matrix was substituted by ones, which make the matrix nonsingular (in fact, this prescription defines well the above additive constant). Comparing with the previous examples presented in Table 1, the accuracy of the approximate solutions remain more or less the same, but the condition numbers are much lower, so that the numerical properties are much more advantageous. This is not surprising; in both cases, the discretization mimics a boundary integral equation, which is a Fredholm equation of the second kind; the operators are compact perturbations of the identity, which implies e.g. that the conjugate gradient method converges quite rapidly: the speed of convergence is superlinear in fact (cf. [2]). This shows that in the case of a pure Dirichlet and a pure Neumann condition, the computational cost of the solution of the discrete systems (10) and (8) can be significantly reduced by applying e.g. the simple conjugate gradient method (or another Krylov subspace method). In addition to this, the problems of dense and ill-conditioned matrices are also avoided. In the case of mixed problems, the above advantages can be preserved if the original mixed problem is converted to ‘pure’ problems, as shown in the next section.
3. Converting mixed problems to pure problems The following technique was originally proposed in connection with the boundary element method, see [6]. Here we briefly recall the main idea but it should be pointed out that the appearing subproblems can be treated in principle by any other method; here the regularized and desingularized MFS are used. For the sake of simplicity, we restrict ourselves to the 2D Laplace equation; note, however, that the technique essentially remains the same for more general and/or 3D problems as well.
Consider the mixed problem: ∂u Δu ¼ 0 in Ω; ujΓD ¼ u0 ; ¼ v0 ∂n ΓN in the Sobolev space H 1 ðΩÞ. Let P be a (not necessarily orthogonal) projector of the closed subspace of the functions of H 1=2 ðΓÞ vanishing along Γ D . Then the operators P 1 ≔I P, P 2 ≔P n (the adjoint of P) are also projectors in the spaces H 1=2 ðΓÞ and H 1=2 ðΓÞ, respectively, and can be interpreted as certain extensions from Γ D to Γ and from Γ N to Γ, respectively. Define the following iteration: ΔU n þ 1=2 ¼ 0; ΔU n þ 1 ¼ 0;
U n þ 1=2 jΓ ¼ un þ P 1 ðu0 un Þ; ∂U n þ 1 ¼ vn þ 1=2 þP 2 ðv0 vn þ 1=2 Þ; ∂n Γ
ð15Þ ð16Þ
(n ¼0, 1, 2,…), where un ≔U n jΓ ;
∂U n þ 1=2 vn þ 1=2 ≔ : ∂n Γ
The first equation is a pure Dirichlet, while the second one is a pure Neumann problem. The iteration can be interpreted as follows: in the first half-step, the Dirichlet boundary condition along Γ D is exactly satisfied, while in the second half-step, it is the Neumann condition along Γ N that is exactly fulfilled. Introducing the Dirichlet-to-Neumann operator (also referred to as Poincaré–Steklov operator) J by Ju≔∂U=∂n, where ΔU ¼ 0 in Ω and UjΓ ≔u, (note that J is an invertible mapping which maps H 1=2 ðΓÞ onto a one-codimensional subspace of H 1=2 ðΓÞ) the above iteration can be written in the form un þ 1 ≔J 1 ðI P 2 ÞJðI P 1 Þun þ b; If the operator where b≔J 1 ðP 2 v0 þ ðI P 2 ÞJP 1 u0 Þ. 1 A≔J ðI P 2 ÞJðI P 1 Þ is a contraction in the space H 1=2 ðΓÞ, the iteration is convergent. The ‘ideal’ choice of the extension operators is as follows. Let P 1 u≔wjΓ , where w is the (unique) solution of the special mixed problem: ∂w Δw ¼ 0 in Ω; wjΓD ¼ ujΓ D ; ¼ 0: ∂n ΓN Similarly, let P 2 u≔∂w=∂n, where now w is the (unique) solution of the mixed problem: ∂w Δw ¼ 0 in Ω; wjΓD ¼ 0; ¼ vjΓ N : ∂n ΓN Then, as can be easily seen, the operator A ¼ J 1 ðI P 2 ÞJðI P 1 Þ equals to the zero operator. That is, the iteration (15) and (16) results in the exact solution after a single iteration step. Unfortunately, to numerically realize these extension operators, mixed problems should be solved, i.e. from a computational point of view, the implementation is at least as costly as that of the solution of the original problem itself. However, in practice, P1, P2 can be defined by a coarse-grid approximation of the above mixed problems, so that the computation of the projections can be embedded in the multi-level context without difficulty. In this case, the operator A is not zero any more, but might remain contractive, so the iteration is still convergent.
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Once the extension operators have been defined, the steps of the iteration can be realized by solving pure Dirichlet and pure Neumann subproblems using the dipole and the monopole formulations as discussed earlier. The extension operators P1, P2 can also be defined in many other ways and also independently of each other. Our experience was that even the extension operators P1, P2 defined by a simple Shepard interpolation technique from Γ D to Γ and a constant extension from Γ N to Γ, respectively, result in excellent convergence properties. To illustrate the method, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as earlier. The problem (13) was supplied with mixed boundary conditions. Along half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. We have implemented the above conversion to a sequence of ‘pure’ subproblems (15) and (16). Table 3 shows the relative L2-errors of the approximate solution computed on the boundary Γ with different numbers of boundary collocation points (N). The parameter iteration indicates the applied number of iterations (15) and (16). The extension operator P1 was defined by Shepard interpolation: ∑x A ΓD uj wkj uk ≔ j ∑xj A ΓD wkj
ðxk A Γ N Þ;
5
ðk þ 1Þ þ 1Þ point xðkÞ , … , xðk of j A Sk the ‘unification’ of some points xj j 1
r
the finer level Sk þ 1 (e.g. xðkÞ j can be defined as the barycentre of the
þ 1Þ þ 1Þ points xðk , … , xðk which are located in a small neighbourj j 1
r
hood of the point xðkÞ j ). Moreover, between the consecutive levels, ‘inter-grid’ transfer operators (restrictions and prolongations) should be defined. This task is also independent of the concrete partial differential equation itself. Suppose that on the level Sk, the components of the vector ðα1 ; …; αNk Þ denote the coefficients of the ðkÞ ðkÞ monopoles (or dipoles) associated with the points xðkÞ 1 ; x2 ; …; xN k .
The simplest projection is defined by averaging r 1 ðk þ 1Þ αðkÞ j ≔ r ∑ αjp p¼1
ðj ¼ 1; 2; …; N k Þ;
while the simplest prolongation is the uniform (piecewise constant) prolongation: þ 1Þ 1 αðk ≔ αðkÞ jp r j
ðp ¼ 1; 2; …; rÞ;
As a next step, the problem has to be discretized at each level, i. e. the matrices A, B, C, Q as well as (8) and (10) should be defined for each collocation point set S0, …, SL. In the simplest cascade technique, only prolongations are needed. The algorithm is as follows:
where wkj ≔1=jjxk xj jj2 . The extension operator P2 was defined to be a constant extension:
On the coarsest level (S0), solve the discrete problem ((8) or
vk ≔const:
Prolongate the solution to the next finer level. Improve the solution by applying e.g. several conjugate gradi-
ðxk A Γ D Þ;
where the const: is defined in such a way that the integral of the extended function v over the boundary vanishes (which is the case when v is the normal derivative of a harmonic function). Comparing with the results of Table 1, it can be seen that the relative errors of the latter methods are achieved after 4 and 5 iteration steps of the solution of ‘pure’ problems (15) and (16). The above ‘pure’ problems can be treated efficiently from a computational point of view, as pointed out in the previous section. Thus, a remarkably efficient method has been obtained. The efficiency can be increased by embedding this iteration into a multi-level technique. 3.1. Multi-level solution To handle the problem in a multi-level way, first of all, a sequence of finite sets of collocation points ðkÞ ðkÞ Sk ≔fxðkÞ 1 ; x2 ; …; xNk g
ðk ¼ 0; 1; …; LÞ;
is needed. Here S0 is considered the ‘coarsest grid’ and SL is the ‘finest grid’ (note that, in the meshless context they need not have any grid or mesh structure). The inclusions Sk Sk þ 1 are not required; however, it is often more comfortable to consider a Table 3 MFS with truncation. Mixed boundary condition. Iteration by converting to pure subproblems. Relative L2-errors on the boundary. N Iteration
8
16
32
64
128
256
512
1 2 3 4 5 6 7 8
2.9570 0.4598 0.5038 0.5018 0.5019 0.5019 0.5019 0.5019
6.0692 0.2344 0.1933 0.1849 0.1855 0.1856 0.1856 0.1856
8.6811 0.3843 0.1042 0.0862 0.0874 0.0874 0.0874 0.0874
11.069 0.5148 0.0739 0.0421 0.0435 0.0435 0.0435 0.0435
13.374 0.6287 0.0668 0.0201 0.0217 0.0216 0.0216 0.0216
15.635 0.7378 0.0681 0.0095 0.0108 0.0107 0.0107 0.0107
17.904 0.8021 0.0691 0.0053 0.0054 0.0053 0.0053 0.0053
(10)) exactly.
ent iterations to the normal equation.
Repeat the procedure until the finest level (SL ) is achieved. The number of conjugate gradient iterations at each level can be kept moderate, which may significantly reduce the computational complexity. Note that, once the ‘inter-grid’ transfer operators have been constructed, the usual items of the multi-level methods (coarse grid correction, multigrid cycle, etc., see e.g. [13] for details) can be defined in a straightforward way. As an example, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as before. The problem (13) was supplied with mixed boundary conditions. Along the half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. We have applied the above cascade type multi-level technique based on monopoles (see (8)). The next finer ‘grid’ was always constructed by doubling the number of boundary collocation points in an equidistant manner. At each level, the approximate solution was defined by applying 4 iteration steps of ‘pure’ problems (15) and (16), where the initial guess was transferred from the previous coarser level. At each ‘pure’ problem, 5 conjugate gradient iteration steps were applied to the corresponding normal equation. Table 4 shows the relative L2-errors of the approximate solution computed on the boundary Γ at the different levels (l) and the corresponding numbers of boundary collocation points (N). The results indicate that the same accuracy can be achieved as in the iterative solution technique (cf. Table 3), but the number of necessary arithmetic operations is significantly reduced.
4. Summary and conclusions Some versions of the method of fundamental solutions have been developed. The source and boundary collocation points were assumed to coincide. This causes singularities, which have been
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
Table 4 MFS with truncation using monopoles and multi-level method. Mixed boundary condition. Relative L2-errors on the boundary. l N Relative L2-error (%)
0 8 0.5018
1 16 0.1856
2 32 0.0875
eliminated by several techniques. The singularities caused by the singularity of the original fundamental solution have been avoided by using a truncated (or other approximate) fundamental solution controlled by a scaling constant. The stronger singularities caused by the Neumann boundary condition i.e. the normal derivatives of the fundamental solution have been eliminated by a special desingularization technique. The same procedure was performed also for the dipole formulation. From a computational point of view, the monopole formulation has proven more advantageous for pure Neumann problems, while the dipole formulation overperforms the traditional monopole formulation for pure Dirichlet problems. For a mixed boundary condition, a special iterative solution algorithm has been presented, converting the original (mixed) problem to a sequence of pure Dirichlet and pure Neumann subproblems, the solutions of which converge rapidly to the solution of the original mixed problem. These pure subproblems can be handled efficiently by using conjugate gradients. The efficiency can be increased further by embedding the resulting method in a multi-level context. In addition to it, the problem of highly ill-conditioned matrices is also avoided. The use of the fast multipole summation technique may speed up the method further.
3 64 0.0436
[2]
[3] [4] [5] [6]
[7]
[8] [9]
[10] [11]
Acknowledgments [12]
The research was partly supported by the European Union (co-financed by the European Social Fund) under the project TÁMOP-4.2.2.A-11/1/KONV-2012-0012. References
[13]
[14]
4 128 0.0216
5 256 0.0109
6 512 0.0057
boundary elements. Proceedings of the 24th international conference on the boundary element method incorporating meshless solution seminar, vol. 13. Southampton, Boston: WitPress; 2002. p. 67–76. Axelsson O, Karátson J. Superlinearly convergent CG methods via equivalent preconditioning for nonsymmetric elliptic operators. Numer Math 2004;99: 197–223. Chen W. Symmetric boundary knot method. Eng Anal Bound Elem 2002;26: 489–94. Chen W, Shen LJ, Shen ZJ, Yuan GW. Boundary knot method for Poisson equations. Eng Anal Bound Elem 2005;29:756–60. Chen W, Wang FZ. A method of fundamental solutions without fictitious boundary. Eng Anal Bound Elem 2010;34:530–2. Gáspár C. Multigrid and multipole techniques in the boundary integral equation method. In: Hackbusch W, Wittum G, editors. Notes on numerical fluid mechanics, vol. 54. Braunschweig, Wiesbaden: Vieweg-Verlag; 1996. p. 102–14. Gáspár C. A multi-level regularized version of the method of fundamental solutions. In: Chen CS, Karageorghis A, Smyrlis YS, editors. The method of fundamental solutions – a meshless method. Atlanta, USA: Dynamic Publishers, Inc.; 2008. p. 145–64. Gáspár C. Regularization techniques for the method of fundamental solutions. Int J Comput Methods 2013;10(2) p. 1341004-1–21. Gáspár C. Combining regularization and desingularization techniques in the method of fundamental solutions. In: Sellier A, Aliabadi MH, editors. Advances in boundary element and meshless techniques XIV. UK: EC, Ltd.; 2013. p. 219–24. Li X. On convergence of the method of fundamental solutions for solving the Dirichlet problem of Poisson's equation. Adv Comput Math 2005;23:265–77. Sarler B. Desingularised method of double layer fundamental solutions for potential flow problems. In: Skerget L, Brebbia CA, editors. Boundary elements and other mesh reduction methods XXX. Southampton, Boston: WitPress; 2008. p. 159–68. Sarler B. Solution of potential flow problems by the modified method of fundamental solutions: formulations with the single layer and the double layer fundamental solutions. Eng Anal Bound Elem 2009;33:1374–82. Stüben K, Trottenberg U. Multigrid methods: fundamental algorithms, model problem analysis and applications. GDM-Studien, No. 96, Birlinghoven, Germany; 1984. Young DL, Chen KH, Lee CW. Novel meshless method for solving the potential problems with arbitrary domain. J Comput Phys 2005;209:290–321.
[1] Alves CJS, Chen CS, Sarler B. The method of fundamental solutions for solving Poisson problems. In: Brebbia CA, Tadeu A, Popov V. Int. series on advances in
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A regularized multi-level technique for solving potential problems by the method of fundamental solutions Csaba Gáspár Széchenyi István University, Egyetem tér 1, H-9026 Győr, Hungary
art ic l e i nf o
a b s t r a c t
Article history: Received 6 October 2013 Accepted 7 May 2014
The method of fundamental solutions is investigated in the case when the source points are located along the boundary of the domain of the original problem and coincide with the collocation points. The appearing singularities are eliminated by several techniques: by using approximate but continuous fundamental solutions (regularization) and via auxiliary subproblems to avoid the stronger singularities that appear in the normal derivatives of the fundamental solution (desingularization). Both monopole and dipole formulations are investigated. A special iterative solution algorithm is presented, which converts the original (mixed) problem to a sequence of pure Dirichlet and pure Neumann subproblems. The pure subproblems can be handled efficiently by using conjugate gradients. The efficiency is significantly increased by embedding the resulting method in a natural multi-level context. At the same time, the problem of the use of highly ill-conditioned matrices is also avoided. & 2014 Published by Elsevier Ltd.
Keywords: Meshless method Method of fundamental solutions Regularization Desingularization Multi-level method
1. Introduction Elliptic partial differential equations play an essential role in a lot of fields of application. Modelling stationary phenomena such as diffusion in fluids or in gases, heat transfer in machines e.g. in parts of traditional or hybrid cars, seepage through porous media lead to solving elliptic problems. The usual implicit time discretization techniques of time-dependent problems result in elliptic problems as well, at every time step. To handle elliptic problems in a meshless way, the popular Boundary Element Method is not suitable, since it requires a boundary mesh structure. The strength of the meshless methods is to circumvent the generation of both domain and boundary mesh. So far, a number of boundary meshless methods have been developed e. g. the boundary knot method [3,4] which uses nonsingular general solutions, or the method of fundamental solutions (MFS, see e.g. [1]), which is based on the fundamental solution of the applied partial differential operator, i.e. on solutions with singularities. In this paper, we restrict ourselves to the MFS applied to homogeneous problem. (For non-homogeneous problems, the approach can be combined with the well-known principle of the Method of Particular Solutions). Consider a second-order elliptic homogeneous linear partial differential equation: Lu ¼ 0
in Ω
E-mail address: [email protected]
ð1Þ
defined in a sufficiently smooth domain Ω supplied with mixed boundary conditions: ∂u ujΓD ¼ u0 ; ¼ v0 ; ð2Þ ∂nΓN where Γ≔∂Ω denotes the boundary of Ω, which has a disjoint decomposition into a Dirichlet part Γ D and a Neumann part Γ N . In its traditional form, the MFS produces an approximate solution of the problem (1) and (2) in the following form: N
uN ðxÞ ¼ ∑ αj Φðx x~ j Þ;
ð3Þ
j¼1
where Φ is a fundamental solution of the operator L, i.e. ΔΦ ¼ δ (here δ denotes the Dirac distribution concentrated at the origin). The predefined points x~ 1 ; …; x~ N (the source points) are located outside of the domain Ω. The a priori unknown coefficients α1 ; …; αN can be computed by enforcing the boundary conditions: N
∑ αj Φðxk x~ j Þ ¼ u0 ðxk Þ
j¼1 N
∑ αj
j¼1
∂Φ ðx x~ j Þ ¼ v0 ðxk Þ ∂nk k
ðxk A Γ D Þ; ð4Þ ðxk A Γ N Þ;
where x1 ; …; xN A Γ are predefined boundary collocation points. Then the function uN exactly satisfies Eq. (1), and exhibits singularities at the source points. Instead of the formulation (3) (called monopole formulation hereafter), it is often more advantageous to use the dipole
http://dx.doi.org/10.1016/j.enganabound.2014.05.002 0955-7997/& 2014 Published by Elsevier Ltd.
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
formulation, where the approximate solution of (1) and (2) is sought in the following form: ð5Þ
j¼1
∂Φ ðx x~ j Þ ¼ u0 ðxk Þ ∂nj k
∂2 Φ ∑ αj ðxk x~ j Þ ¼ v0 ðxk Þ j ¼ 1 ∂nk ∂nj
ð6Þ ðxk A Γ N Þ:
Note that the monopole and the dipole formulations can be considered as meshless discretizations of the indirect BEM based on single layer and double layer potentials, respectively. A common disadvantage of the traditional forms (3) and (5) is the use of external source points, the location of which can be hardly automatized. Though the MFS has excellent accuracy properties (see e.g. [10] and references therein), the systems (4) and (6) are highly illconditioned in general, which is a severe drawback of the method. This is the case especially when the source points are located far from the boundary. On the other hand, if they are close to the boundary, the systems (4) and (6) become much better conditioned, however, the accuracy goes wrong due to the appearance of numerical singularities at the boundary collocation points. A usual technique is to define the source and the boundary collocation points to coincide. Special techniques are required to avoid the problem of singularity (regularization and desingularization, see e.g. [14,7,12,8]). In this paper, we investigate some regularized versions of both the monopole and the dipole formulations. It turns out that, from a computational point of view, the dipole formulation is much more advantageous for handling pure Dirichlet problems, while for pure Neumann problems, the monopole formulation overperforms the dipole formulation. For mixed boundary conditions, a special iterative technique is proposed which converts the original mixed problem to a convergent sequence of pure Dirichlet and pure Neumann subproblems. This results in a computationally efficient method, which avoids also the problem of highly ill-conditioned linear systems. The efficiency can be increased further by embedding the method in a natural multi-level context.
In the rest of the paper, suppose that the source points and the boundary collocation points x1 ; …; xN A Γ coincide. Then the monopole formulation has the form uN ðxÞ ¼ ∑ αj Φðx xj Þ;
ð7Þ
j¼1
N
∑ αj Akj ¼ u0 ðxk Þ
j¼1
ðxk A Γ D Þ; ð8Þ
∑ αj Bkj ¼ v0 ðxk Þ
ðxk A Γ N Þ:
Similarly, the dipole formulation has the form N
uN ðxÞ ¼ ∑ αj j¼1
Here the entries of the matrices A, B, C, Q are defined as follows:
C kj ¼
∂Φ ðx xj Þ; ∂nj
∂Φ ðx xj Þ; ∂nj k
Bkj ¼
∂Φ ðx xj Þ; ∂nk k
Q kj ¼
∂2 Φ ðx xj Þ ∂nk ∂nj k
ð11Þ
for j ak. Due to the singularity of the fundamental solution at the origin, the diagonal entries cannot be computed by the above definition. For the proper definition of Akk, one should replace the fundamental solution Φ with an approximate fundamental solution Φ, which has no singularity at the origin. Such an approximate fundamental solution can be defined e.g. by truncation. In polar coordinates: 8 1 > > > if r Z < ΦðrÞ c ð12Þ ΦðrÞ≔ 1 1 > > if r o ; Φ > : c c provided that Φ is a radial function i.e. it depends only on r ¼ jjxjj, which is often the case. Here c denotes a carefully chosen scaling constant which should remain inversely proportional to the characteristic distance of the boundary collocation points, when N varies. Another regularization technique is to replace Φ with the fundamental solution of the singularly perturbed fourth-order operator LðI ð1=c2 ÞLÞ, where I denotes the identity operator and c is again a scaling constant, see [8] for details. Thus, the diagonal terms Akk can be computed without difficulty. Using the simplest truncation, Akk ¼ Φð0Þ, while for j a k, Akj ¼ Φðxk xj Þ. The proper definition of the diagonal terms Bkk is somewhat more difficult since the derivatives of Φ have stronger singularities at the origin than Φ itself. Let w be a smooth, easily computable particular solution of Eq. (1), then w can be approximated by the monopole formulation: N
wN ðxÞ≔ ∑ βj Φðx xj Þ; j¼1
N ∂wN ∂Φ ðxÞ ¼ ∑ βj ðx xj Þ: ∂n ∂n j¼1
Hence ∑ βj Bkj þ βk Bkk ¼
jak
where the coefficients α1 ; …; αN can be computed by solving the algebraic system:
j¼1
j¼1
ðxk A Γ N Þ:
from where the coefficients β1 ; …; βN can be computed by solving (1) supplied with a pure Dirichlet condition. Computing the normal derivative of the particular solution w, we have the following:
2. Regularization and desingularization
N
ðxk A Γ D Þ; ð10Þ
∑ αj Q kj ¼ v0 ðxk Þ
Akj ¼ Φðxk xj Þ;
ðxk A Γ D Þ;
N
N
j¼1 N
Again, Φ is the fundamental solution of the operator L, and the source points x~ 1 ; …; x~ N are located outside of Ω. Now the coefficients α1 ; …; αN can be computed by solving the linear system: N
N
∑ αj C kj ¼ u0 ðxk Þ
∂Φ uN ðxÞ ¼ ∑ αj ðx x~ j Þ: ∂n j j¼1 N
∑ αj
where the coefficients α1 ; …; αN solve the algebraic system:
ð9Þ
∂w ðx Þ ∂nk k
ðk ¼ 1; 2; …; NÞ:
Thus, the diagonal terms Bkk can be defined as ! 1 ∂w ∑ βj Bkj þ ðxk Þ : Bkk ≔ βk ∂nk jak This is the desingularization idea, see e.g. [14,11,5] for details. Note that in the simplest case of the Laplace equation w can often be chosen e.g. by w : 1; in the case of the modified 2D Helmholtz equation Δu λ2 u ¼ 0, a particular solution wðxÞ≔I 0 ðλ jjxjjÞ can be used, where I0 denotes the familiar modified Bessel function of the first kind, and so on. By a proper choice of the particular solution w, one can ensure that the coefficients βk do not vanish, so that the
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
definition of Bkk is correct (using possibly several particular solutions). The desingularization can be improved by redefining not only the diagonal terms but also the neighbouring entries [9]. In the dipole formulation, a straightforward desingularization is C kk ≔Bkk (while for ja k, C kj ¼ Bjk ). Finally, the diagonal terms Qkk can be defined by a desingularization procedure similar to the monopole formulation, by solving an auxiliary pure Dirichlet problem. See [8] for details. Remark. The definition of the diagonal entries C kk ≔Bkk is straightforward but not trivial, so that it needs some explanation. For simplicity, assume that Ω is a polygon with sides Γj (j ¼ 1; 2; …; N) and the boundary collocation point xj is the midpoint of Γj (more precisely let us approximate Ω in this way provided that N is large enough). Let us express the solution of the original problem in the form of a single layer potential: Z 1 uðξÞ ¼ ðlog ξ yÞ wðyÞ dΓ y ðξ A ΩÞ; 2π Γ with some continuous density function w. It is well known that the normal derivative of a single layer potential has a jump at the boundary i.e. Z ∂u 1 〈xk y; nk 〉 1 lim ðξÞ ¼ wðyÞ dΓ y wðxk Þ: ξ-xk ∂nk 2π Γ jjxk yjj2 2 (Here 〈; 〉 denotes the Euclidean scalar product.) The integral in the right-hand side can be expressed as a sum of integrals taken along the sides Γj, whence Z ∂u 1 〈xk y; nk 〉 ðξÞ ¼ ∑ wðyÞ dΓ y lim 2 ξ-xk ∂nk j a k2π Γ j jjxk yjj Z 1 〈xk y; nk 〉 1 þ wðyÞ dΓ y wðxk Þ: 2π Γ k jjxk yjj2 2 The last integral (taken along Γk) vanishes since xk y is orthogonal to the normal vector nk. Approximating the preceding integrals by the midpoint rule, we obtain the following discrete form of the equations: 1 ∑ Bkj hj wj wk ¼ vk ; 2 jak where wj ≔wðxj Þ, hj ≔jΓ j j (the length of Γj) and Bkj ¼
1 〈xk xj ; nk 〉 ∂Φ ¼ ðx xj Þ: 2π jjxk xj jj2 ∂nk k
Denoting by αj ≔hj wj , we have ∑ Bkj αj
jak
1 α ¼ vk ; 2hk k
whence Bkk ¼ 1=ð2hk Þ. On the other hand, expressing the solution of the original problem in the form of a double layer potential: Z 1 〈ξ y; ny 〉 uðξÞ ¼ wðyÞ dΓ y ðξ A ΩÞ; 2π Γ jjξ yjj2 (with a density function w which is different from the previous case). Using the well-known theorem of the jump of the double
3
layer potential at the boundary, we have Z 1 〈xk y; ny 〉 1 lim uðξÞ ¼ wðyÞ dΓ y wðxk Þ ¼ uk : ξ-xk 2π Γ jjxk yjj2 2 A procedure similar to the case of the single layer potential yields that these equations can be approximated by the discrete forms: ∑ C kj αj
jak
1 α ¼ uk ; 2hk k
where 1 〈xk xj ; nj 〉 C kj ≔ 2π jjxk xj jj2
ðj a kÞ
and the diagonal entries are C kk ¼ 1=ð2hk Þ. Thus, C kk ¼ Bkk , and for ja k: Bjk ¼
1 〈xj xk ; nj 〉 ¼ C kj : 2π jjxj xk jj2
Note that the above definitions of Bkk and Ckk mimic the jump of the normal derivative of the single layer potential at the boundary and the jump of the double layer potential, respectively. Both the monopole and the dipole formulations are suitable for solving mixed boundary problems using a regularized fundamental solution; the singularities of the diagonal entries of the systems (8) and (10) are avoided by regularization and desingularization. To illustrate the above techniques, consider the simplest 2D Laplace equation Δu ¼ 0
ð13Þ
in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution uðxÞ≔log ððxð1Þ þ2Þ2 þ ðxð2Þ þ 3Þ2 Þ; ð1Þ
ð14Þ
ð2Þ
where x , x denote the components of the vector variable x. The boundary was discretized by N collocation points in an equidistant way. The scaling parameter was set as follows: c≔N, which is a quasioptimal choice (more precisely, define c≔const:=h, where h denotes the characteristic distance of the boundary collocation points, and const: means a dimensionless constant which is independent of the discretization of the boundary). The problem was supplied with mixed boundary conditions; along the half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. Table 1 shows the relative L2-errors of the approximate solution computed on the boundary Γ and also the condition numbers of the corresponding algebraic systems. The method applied was the desingularized MFS based on monopoles (7) and dipoles (9). This seems an acceptable compromise: the accuracy is still acceptable and the condition numbers are still moderate; they are in fact much less than in the case of a traditional MFS with external source points located from the boundary at a distance which is comparable with the size of the domain. Remark. For Neumann and mixed problems, it is sufficient to compute the error of the approximate solution along the boundary only. Once the boundary values of the solution have been (approximately) computed, the solution inside the domain can be reconstructed by solving an additional pure Dirichlet problem.
Table 1 MFS with truncation. Mixed boundary conditions. Relative L2-errors on the boundary and condition numbers. N
8
16
32
64
128
256
512
Monopoles: relative L2-error (%) Condition number Dipoles: relative L2-error (%) Condition number
0.4541 8 0.4982 3
0.2009 21 0.1516 6
0.1058 53 0.0458 13
0.0552 127 0.0177 29
0.0282 294 0.0081 63
0.0142 670 0.0040 135
0.0071 1.5E þ3 0.0020 287
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
Table 2 MFS with truncation and desingularization. Pure Neumann boundary condition (using monopoles) and pure Dirichlet boundary condition (using dipoles). Relative L2-errors on the boundary and condition numbers. N
8
16
32
64
128
256
512
Monopoles: relative L2-error (%) Condition number Dipoles: relative L2-error (%) Condition number
0.5580 15 4.8602 3
0.3566 17 1.3709 3
0.2057 18 0.4715 3
0.1111 19 0.1869 4
0.0568 21 0.0821 4
0.0286 30 0.0385 4
0.0143 44 0.0186 4
This latter problem can be treated by a regularized MFS and needs no desingularization. The condition numbers dramatically decrease if the above methods are applied to pure problems. More precisely, the regularizeddesingularized MFS based on monopoles (7) is quite efficient for solving pure Neumann problems, while the dipole formulation (9) overperforms the regularized MFS when solving pure Dirichlet problems. This is in accordance with the indirect boundary element methods: pure Dirichlet (resp. Neumann) problems should be treated by double layer (resp. single layer) potentials. To illustrate this, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as earlier. Table 2 shows the relative L2-errors of the approximate solution computed on the boundary Γ and the condition numbers of the corresponding algebraic system (8), when a pure Neumann condition was prescribed. Here the monopole formulation was applied. Table 2 shows also the relative L2-errors of the normal derivative of the approximate solution computed on the boundary Γ and also the condition numbers of the corresponding algebraic system (10), when a pure Dirichlet condition was prescribed. Here the dipole formulation was applied. Since the algebraic system is always singular (due to the fact that the solution of the pure Neumann problem is unique up to an additive constant only), the first row of the matrix was substituted by ones, which make the matrix nonsingular (in fact, this prescription defines well the above additive constant). Comparing with the previous examples presented in Table 1, the accuracy of the approximate solutions remain more or less the same, but the condition numbers are much lower, so that the numerical properties are much more advantageous. This is not surprising; in both cases, the discretization mimics a boundary integral equation, which is a Fredholm equation of the second kind; the operators are compact perturbations of the identity, which implies e.g. that the conjugate gradient method converges quite rapidly: the speed of convergence is superlinear in fact (cf. [2]). This shows that in the case of a pure Dirichlet and a pure Neumann condition, the computational cost of the solution of the discrete systems (10) and (8) can be significantly reduced by applying e.g. the simple conjugate gradient method (or another Krylov subspace method). In addition to this, the problems of dense and ill-conditioned matrices are also avoided. In the case of mixed problems, the above advantages can be preserved if the original mixed problem is converted to ‘pure’ problems, as shown in the next section.
3. Converting mixed problems to pure problems The following technique was originally proposed in connection with the boundary element method, see [6]. Here we briefly recall the main idea but it should be pointed out that the appearing subproblems can be treated in principle by any other method; here the regularized and desingularized MFS are used. For the sake of simplicity, we restrict ourselves to the 2D Laplace equation; note, however, that the technique essentially remains the same for more general and/or 3D problems as well.
Consider the mixed problem: ∂u Δu ¼ 0 in Ω; ujΓD ¼ u0 ; ¼ v0 ∂n ΓN in the Sobolev space H 1 ðΩÞ. Let P be a (not necessarily orthogonal) projector of the closed subspace of the functions of H 1=2 ðΓÞ vanishing along Γ D . Then the operators P 1 ≔I P, P 2 ≔P n (the adjoint of P) are also projectors in the spaces H 1=2 ðΓÞ and H 1=2 ðΓÞ, respectively, and can be interpreted as certain extensions from Γ D to Γ and from Γ N to Γ, respectively. Define the following iteration: ΔU n þ 1=2 ¼ 0; ΔU n þ 1 ¼ 0;
U n þ 1=2 jΓ ¼ un þ P 1 ðu0 un Þ; ∂U n þ 1 ¼ vn þ 1=2 þP 2 ðv0 vn þ 1=2 Þ; ∂n Γ
ð15Þ ð16Þ
(n ¼0, 1, 2,…), where un ≔U n jΓ ;
∂U n þ 1=2 vn þ 1=2 ≔ : ∂n Γ
The first equation is a pure Dirichlet, while the second one is a pure Neumann problem. The iteration can be interpreted as follows: in the first half-step, the Dirichlet boundary condition along Γ D is exactly satisfied, while in the second half-step, it is the Neumann condition along Γ N that is exactly fulfilled. Introducing the Dirichlet-to-Neumann operator (also referred to as Poincaré–Steklov operator) J by Ju≔∂U=∂n, where ΔU ¼ 0 in Ω and UjΓ ≔u, (note that J is an invertible mapping which maps H 1=2 ðΓÞ onto a one-codimensional subspace of H 1=2 ðΓÞ) the above iteration can be written in the form un þ 1 ≔J 1 ðI P 2 ÞJðI P 1 Þun þ b; If the operator where b≔J 1 ðP 2 v0 þ ðI P 2 ÞJP 1 u0 Þ. 1 A≔J ðI P 2 ÞJðI P 1 Þ is a contraction in the space H 1=2 ðΓÞ, the iteration is convergent. The ‘ideal’ choice of the extension operators is as follows. Let P 1 u≔wjΓ , where w is the (unique) solution of the special mixed problem: ∂w Δw ¼ 0 in Ω; wjΓD ¼ ujΓ D ; ¼ 0: ∂n ΓN Similarly, let P 2 u≔∂w=∂n, where now w is the (unique) solution of the mixed problem: ∂w Δw ¼ 0 in Ω; wjΓD ¼ 0; ¼ vjΓ N : ∂n ΓN Then, as can be easily seen, the operator A ¼ J 1 ðI P 2 ÞJðI P 1 Þ equals to the zero operator. That is, the iteration (15) and (16) results in the exact solution after a single iteration step. Unfortunately, to numerically realize these extension operators, mixed problems should be solved, i.e. from a computational point of view, the implementation is at least as costly as that of the solution of the original problem itself. However, in practice, P1, P2 can be defined by a coarse-grid approximation of the above mixed problems, so that the computation of the projections can be embedded in the multi-level context without difficulty. In this case, the operator A is not zero any more, but might remain contractive, so the iteration is still convergent.
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Once the extension operators have been defined, the steps of the iteration can be realized by solving pure Dirichlet and pure Neumann subproblems using the dipole and the monopole formulations as discussed earlier. The extension operators P1, P2 can also be defined in many other ways and also independently of each other. Our experience was that even the extension operators P1, P2 defined by a simple Shepard interpolation technique from Γ D to Γ and a constant extension from Γ N to Γ, respectively, result in excellent convergence properties. To illustrate the method, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as earlier. The problem (13) was supplied with mixed boundary conditions. Along half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. We have implemented the above conversion to a sequence of ‘pure’ subproblems (15) and (16). Table 3 shows the relative L2-errors of the approximate solution computed on the boundary Γ with different numbers of boundary collocation points (N). The parameter iteration indicates the applied number of iterations (15) and (16). The extension operator P1 was defined by Shepard interpolation: ∑x A ΓD uj wkj uk ≔ j ∑xj A ΓD wkj
ðxk A Γ N Þ;
5
ðk þ 1Þ þ 1Þ point xðkÞ , … , xðk of j A Sk the ‘unification’ of some points xj j 1
r
the finer level Sk þ 1 (e.g. xðkÞ j can be defined as the barycentre of the
þ 1Þ þ 1Þ points xðk , … , xðk which are located in a small neighbourj j 1
r
hood of the point xðkÞ j ). Moreover, between the consecutive levels, ‘inter-grid’ transfer operators (restrictions and prolongations) should be defined. This task is also independent of the concrete partial differential equation itself. Suppose that on the level Sk, the components of the vector ðα1 ; …; αNk Þ denote the coefficients of the ðkÞ ðkÞ monopoles (or dipoles) associated with the points xðkÞ 1 ; x2 ; …; xN k .
The simplest projection is defined by averaging r 1 ðk þ 1Þ αðkÞ j ≔ r ∑ αjp p¼1
ðj ¼ 1; 2; …; N k Þ;
while the simplest prolongation is the uniform (piecewise constant) prolongation: þ 1Þ 1 αðk ≔ αðkÞ jp r j
ðp ¼ 1; 2; …; rÞ;
As a next step, the problem has to be discretized at each level, i. e. the matrices A, B, C, Q as well as (8) and (10) should be defined for each collocation point set S0, …, SL. In the simplest cascade technique, only prolongations are needed. The algorithm is as follows:
where wkj ≔1=jjxk xj jj2 . The extension operator P2 was defined to be a constant extension:
On the coarsest level (S0), solve the discrete problem ((8) or
vk ≔const:
Prolongate the solution to the next finer level. Improve the solution by applying e.g. several conjugate gradi-
ðxk A Γ D Þ;
where the const: is defined in such a way that the integral of the extended function v over the boundary vanishes (which is the case when v is the normal derivative of a harmonic function). Comparing with the results of Table 1, it can be seen that the relative errors of the latter methods are achieved after 4 and 5 iteration steps of the solution of ‘pure’ problems (15) and (16). The above ‘pure’ problems can be treated efficiently from a computational point of view, as pointed out in the previous section. Thus, a remarkably efficient method has been obtained. The efficiency can be increased by embedding this iteration into a multi-level technique. 3.1. Multi-level solution To handle the problem in a multi-level way, first of all, a sequence of finite sets of collocation points ðkÞ ðkÞ Sk ≔fxðkÞ 1 ; x2 ; …; xNk g
ðk ¼ 0; 1; …; LÞ;
is needed. Here S0 is considered the ‘coarsest grid’ and SL is the ‘finest grid’ (note that, in the meshless context they need not have any grid or mesh structure). The inclusions Sk Sk þ 1 are not required; however, it is often more comfortable to consider a Table 3 MFS with truncation. Mixed boundary condition. Iteration by converting to pure subproblems. Relative L2-errors on the boundary. N Iteration
8
16
32
64
128
256
512
1 2 3 4 5 6 7 8
2.9570 0.4598 0.5038 0.5018 0.5019 0.5019 0.5019 0.5019
6.0692 0.2344 0.1933 0.1849 0.1855 0.1856 0.1856 0.1856
8.6811 0.3843 0.1042 0.0862 0.0874 0.0874 0.0874 0.0874
11.069 0.5148 0.0739 0.0421 0.0435 0.0435 0.0435 0.0435
13.374 0.6287 0.0668 0.0201 0.0217 0.0216 0.0216 0.0216
15.635 0.7378 0.0681 0.0095 0.0108 0.0107 0.0107 0.0107
17.904 0.8021 0.0691 0.0053 0.0054 0.0053 0.0053 0.0053
(10)) exactly.
ent iterations to the normal equation.
Repeat the procedure until the finest level (SL ) is achieved. The number of conjugate gradient iterations at each level can be kept moderate, which may significantly reduce the computational complexity. Note that, once the ‘inter-grid’ transfer operators have been constructed, the usual items of the multi-level methods (coarse grid correction, multigrid cycle, etc., see e.g. [13] for details) can be defined in a straightforward way. As an example, consider again the problem (13) in the square Ω≔½ 1; 1 ½ 1; 1 with the test solution (14). Set c≔N, as before. The problem (13) was supplied with mixed boundary conditions. Along the half of the boundary, a Dirichlet condition was prescribed, while along the remaining part of the boundary, a Neumann condition was imposed. We have applied the above cascade type multi-level technique based on monopoles (see (8)). The next finer ‘grid’ was always constructed by doubling the number of boundary collocation points in an equidistant manner. At each level, the approximate solution was defined by applying 4 iteration steps of ‘pure’ problems (15) and (16), where the initial guess was transferred from the previous coarser level. At each ‘pure’ problem, 5 conjugate gradient iteration steps were applied to the corresponding normal equation. Table 4 shows the relative L2-errors of the approximate solution computed on the boundary Γ at the different levels (l) and the corresponding numbers of boundary collocation points (N). The results indicate that the same accuracy can be achieved as in the iterative solution technique (cf. Table 3), but the number of necessary arithmetic operations is significantly reduced.
4. Summary and conclusions Some versions of the method of fundamental solutions have been developed. The source and boundary collocation points were assumed to coincide. This causes singularities, which have been
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i
C. Gáspár / Engineering Analysis with Boundary Elements ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
Table 4 MFS with truncation using monopoles and multi-level method. Mixed boundary condition. Relative L2-errors on the boundary. l N Relative L2-error (%)
0 8 0.5018
1 16 0.1856
2 32 0.0875
eliminated by several techniques. The singularities caused by the singularity of the original fundamental solution have been avoided by using a truncated (or other approximate) fundamental solution controlled by a scaling constant. The stronger singularities caused by the Neumann boundary condition i.e. the normal derivatives of the fundamental solution have been eliminated by a special desingularization technique. The same procedure was performed also for the dipole formulation. From a computational point of view, the monopole formulation has proven more advantageous for pure Neumann problems, while the dipole formulation overperforms the traditional monopole formulation for pure Dirichlet problems. For a mixed boundary condition, a special iterative solution algorithm has been presented, converting the original (mixed) problem to a sequence of pure Dirichlet and pure Neumann subproblems, the solutions of which converge rapidly to the solution of the original mixed problem. These pure subproblems can be handled efficiently by using conjugate gradients. The efficiency can be increased further by embedding the resulting method in a multi-level context. In addition to it, the problem of highly ill-conditioned matrices is also avoided. The use of the fast multipole summation technique may speed up the method further.
3 64 0.0436
[2]
[3] [4] [5] [6]
[7]
[8] [9]
[10] [11]
Acknowledgments [12]
The research was partly supported by the European Union (co-financed by the European Social Fund) under the project TÁMOP-4.2.2.A-11/1/KONV-2012-0012. References
[13]
[14]
4 128 0.0216
5 256 0.0109
6 512 0.0057
boundary elements. Proceedings of the 24th international conference on the boundary element method incorporating meshless solution seminar, vol. 13. Southampton, Boston: WitPress; 2002. p. 67–76. Axelsson O, Karátson J. Superlinearly convergent CG methods via equivalent preconditioning for nonsymmetric elliptic operators. Numer Math 2004;99: 197–223. Chen W. Symmetric boundary knot method. Eng Anal Bound Elem 2002;26: 489–94. Chen W, Shen LJ, Shen ZJ, Yuan GW. Boundary knot method for Poisson equations. Eng Anal Bound Elem 2005;29:756–60. Chen W, Wang FZ. A method of fundamental solutions without fictitious boundary. Eng Anal Bound Elem 2010;34:530–2. Gáspár C. Multigrid and multipole techniques in the boundary integral equation method. In: Hackbusch W, Wittum G, editors. Notes on numerical fluid mechanics, vol. 54. Braunschweig, Wiesbaden: Vieweg-Verlag; 1996. p. 102–14. Gáspár C. A multi-level regularized version of the method of fundamental solutions. In: Chen CS, Karageorghis A, Smyrlis YS, editors. The method of fundamental solutions – a meshless method. Atlanta, USA: Dynamic Publishers, Inc.; 2008. p. 145–64. Gáspár C. Regularization techniques for the method of fundamental solutions. Int J Comput Methods 2013;10(2) p. 1341004-1–21. Gáspár C. Combining regularization and desingularization techniques in the method of fundamental solutions. In: Sellier A, Aliabadi MH, editors. Advances in boundary element and meshless techniques XIV. UK: EC, Ltd.; 2013. p. 219–24. Li X. On convergence of the method of fundamental solutions for solving the Dirichlet problem of Poisson's equation. Adv Comput Math 2005;23:265–77. Sarler B. Desingularised method of double layer fundamental solutions for potential flow problems. In: Skerget L, Brebbia CA, editors. Boundary elements and other mesh reduction methods XXX. Southampton, Boston: WitPress; 2008. p. 159–68. Sarler B. Solution of potential flow problems by the modified method of fundamental solutions: formulations with the single layer and the double layer fundamental solutions. Eng Anal Bound Elem 2009;33:1374–82. Stüben K, Trottenberg U. Multigrid methods: fundamental algorithms, model problem analysis and applications. GDM-Studien, No. 96, Birlinghoven, Germany; 1984. Young DL, Chen KH, Lee CW. Novel meshless method for solving the potential problems with arbitrary domain. J Comput Phys 2005;209:290–321.
[1] Alves CJS, Chen CS, Sarler B. The method of fundamental solutions for solving Poisson problems. In: Brebbia CA, Tadeu A, Popov V. Int. series on advances in
Please cite this article as: Gáspár C. A regularized multi-level technique for solving potential problems by the method of fundamental solutions. Eng. Anal. Boundary Elem. (2014), http://dx.doi.org/10.1016/j.enganabound.2014.05.002i