Meshless Galerkin method using radial basis functions based on domain decomposition

Applied Mathematics and Computation 179 (2006) 750–762 www.elsevier.com/locate/amc

Meshless Galerkin method using radial basis functions based on domain decomposition Yong Duan School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu 610054, PR China

Abstract In this present paper, we will give an overview about the combination of RBF meshless method with domain decomposition, especially Galerkin method. Based on the analysis for this method before, we will consider the combination of projection domain decomposition method with radial basis function. Both of Galerkin approximation and collocation method are discussed. Some examples are given to show the efficiency of this method.  2005 Elsevier Inc. All rights reserved. Keywords: Radial basis function; Projection domain decomposition; Collocation; Meshless; Galerkin; Parallel

1. Introduction In the last decade, the so-called meshless method for partial differential equations is focused on because that this method can deal with those problem with complex domain or high dimensional domain. Moreover, it can dramatically save time in mesh-generation such as finite element method (FEM). Such complication has been described recently by Griebel and Schweitzer [11] as ‘‘a very time-consuming portion of overall computation is the mesh-generation from CAD input data. Typically, more than 70 percent of overall computation is spent by mesh generators’’. The popular meshless methods include moving least square method (see [4,5]), generalized finite element method (GFEM) (see [3]) and radial basis function (RBF) method (see [6,7]). We will restrict us to RBF meshless method in this paper. RBFs are a powerful tool in multi-variable approximation. Schaback and his group have done many researches about the interpolation theory using RBF and its applications for PDEs, see, for example, [19– 22,29,30]. Compared with the classical numerical method, such as finite difference method (FDM), finite volume method (FVM), finite element method (FEM), the RBF does not need the construction of mesh, which is a complicated work for FEM. Actually, it is a true meshless method (see [13]). RBF meshless methods , including collocation and Galerkin method, are discussed in a series of papers, for example [8,9,23,27]. Compared with collocation method, the Galerkin approximation cannot carry out for both Dirichlet and mixed problem because that the finite dimensional space spanned by RBF does not belong to H 1C ðXÞ, although it has a better E-mail address: [email protected] 0096-3003/$ - see front matter  2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.11.153

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

751

performance in theoretical analysis [23] for Neumann problem of Helmholtz equation. Some skills, like penalty method, Lagrange method, are advised to deal with Dirichlet problem and mixed problem to get better numerical solution at the cost of larger condition number than before [27]. It is well known that the stiffness matrix for RBF meshless method will be highly ill-conditioned and cannot be expected well-conditioned whenever higher accuracy is expected because of the uncertainty [19] property of RBF. Although RBF has these defects, its perspective is mainly charming. Considering the ill-condition property of RBF, which may hinder the applications for large scale problem, DDM should be a instinctive plea for the meshless method using RBF. Of course, parallel computation is another important application. Hon [14,24] consider the overlapping and non-overlapping domain decomposition using RBF collocation. The numerical result is quite excited. However, the theoretical analysis of this method is absent. Actually, we cannot prove that the linear system arising from RBF collocation is solvable. Wu [28] consider the additive Schwarz domain decomposition and obtain an estimate like which in FEM, with a constant dependent on density h. Duan [25,26] had tried combine RBF meshless Galerkin method with nonoverlapping DDM. For Dirichlet–Neumann iteration, the convergent rate is dramatically dependent on density h. Also, the numerical result is not good, especially for the point near corner. According to the property of RBF meshless Galerkin method, we hope that the DDM should be: 1. Without iteration, which is due to the bad condition number of RBF meshless methods. 2. We can obtain the error estimate. For Robin iteration, we cannot. 3. Whenever we consider Helmholtz equations with Neumann boundary condition, every sub-problem should be Neumann type. This can avoid the approximation for essential BVP. Based on these request, we will consider an improved projection domain decomposition (PDM), opposed to the original PDM. Projection domain decomposition method (PDM) is introduced firstly by Agoshkov et al. in [1]. The most important aspect of this method is how to obtain the well-conditioned, piecewise polynomial bases of the corresponding space arising in non-overlapping domain decomposition. In a series of papers [10,16,17], the possibility to construct these bases in a fast and accurate way has been proven for the Laplace and Stokes operators. Compared with another domain decomposition method (DDM), the PDM need neither the suitable choice of the acceleration parameters to ensure the convergent rate of the iteration arising in DDM nor the iteration at each subdomains. This property is meaningful for meshless method (using radial basis function (RBF)) coupled with DDM. In this paper, we will combine RBF meshless method with PDM. All analysis will focus on RBF meshless Galerkin method. Considering that collocation method is more successful in numerical computation than Galerkin approximation, we will introduce the collocation method and give some examples to show the efficiency. This paper is organized as follows. In Section 2.1, we will introduce the improved PDM at continuous case. In Section 2.2, we will introduce the RBF finite space and give the error estimate of this method. In Section 2.3, we will give some examples for this method. In Section 3, we will introduce RBF collocation method without theoretical analysis. Also, some examples will be given. 2. PDM using RBF Galerkin methods 2.1. Non-overlapping domain decomposition Consider the model problem: ( Du þ u ¼ f in X; ou ¼ 0 on oX; on

ð1Þ

where X is a d-dimensional domain, with a Lipschitz boundary oX, whose outer unit normal direction is denoted by n, f is a given function of L2(X). Assuming that X is conformally partitioned into two

752

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

non-overlapping convex subdomains X1 and X2, and denoting by C ¼ X1 \ X2 . We indicate by ui the restriction to Xi, i = 1, 2, of the solution u of Eq. (1), and ni by the normal direction on oXi \ C, oriented outward X2. Set n = n1. A standard domain decomposition skill [18] leads to the following equation: 8 Du1 þ u1 ¼ f in X1 ; > > T > ou1 > ¼0 on oX1 oX; > > on > > < u1 ¼ u2 on C; ð2Þ ou2 ou1 > on C; > on ¼ on > > T > > ou2 ¼ 0 on oX2 oX; > on > : Du2 þ u2 ¼ f in X2 . Denoted by k ¼ ou , the unknown trace of the normal derivatives of the solution u on C, ui ¼ ui þ Hi k, where onjC 8 Du þ ui ¼ f in Xi ; > > <  i T oui ¼0 on oXi oX; ð3Þ on > > : oui ¼ 0 on C. on

Hi k is the harmonic extension of k in Xi 8 > < DHi k þ Hi k ¼ 0 in Xi ; T oHi k ¼0 on oXi oX; on > : oHi k ¼k on C. on

ð4Þ

The multi-domain form is equivalent to the original equation if and only if, the following Poincare´–Steklov equation holds, Tk ¼ v on C; where v ¼

ou2 on



ou1 . on

ð5Þ

T is entitled by Poincare´–Steklov operator, opposed to the Steklov–Poincare´ operator.

12

1

T : H ðCÞ 7! H 2 ðCÞ : Tk ¼ H1 k  H2 k.

Remark 2.1. Operator T had been extensively discussed in FEM, which maps Neumann data to Dirichlet data. For the discussion in FEM, see [15]. Introduce the bilinear form, "u, v 2 H1(X) Z Aðu; vÞ ¼ ru  rv þ uv; ZX Ai ðu; vÞ ¼ ru  rv þ uv Xi

and bounded linear functional Z Z FðvÞ ¼ fv; Fi ðvÞ ¼ fv. X

Xi 1

Theorem 2.1. Operator T is coercive and continuous i.e., there exist constants C1, C2, such that 8k; g 2 H 2 ðCÞ hTk; gi 6 C 1 kgk C 2 kkk

2 H

1 2 ðCÞ

H

1 2 ðCÞ

kkk

6 hTk; ki.

H

1 2 ðCÞ

;

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

753

Proof Z

Z Z oH1 g oH2 g ðH1 k  H2 kÞg ¼ H1 k þ H2 k on on2 C C C C Z Z oH1 g oH2 g þ ¼ A1 ðH1 k; H1 gÞ þ A2 ðH2 k; H2 gÞ. ¼ H1 k H2 k on on oX1 oX2

hTk; gi ¼

Tkg ¼

Z

Because that bilinear form Ai is coercive and continuous, we have hTk; gi 6 C 1 ðkH1 gkH 1 ðX1 Þ kH1 kkH 1 ðX1 Þ þ kH2 gkH 1 ðX2 Þ kH2 kkH 1 ðX2 Þ Þ. By using the priori estimates, we have hTk; gi 6 C 1 kkk

H

1 2 ðCÞ

kgk

H

1 2 ðCÞ

.

By using the trace theorem, we have hTg; gi P C 2 kgk This finish the proof.

2 H

1 2 ðCÞ

.

h

Remark 2.2. For Dirichlet BVP  Lu ¼ f in X; u¼0

ð6Þ

on oX;

where L is any elliptic operator, and oX 2 C1 . By using the penalty methods [2], we have the following approximation: ( Lua ¼ f in X; ð7Þ oua þ aua ¼ 0 on oX. onL Also, the following error estimate holds C ku  ua;h kH 1 ðXÞ 6 pffiffiffi kf kH 0 ðXÞ þ Cahkþ1 kf kH k ðXÞ ; a

ð8Þ

where ua,h is the RBF-solution for Eq. (7). Remark 2.3. We can also consider Lagrange multipliers to deal with Eq. (1) with Dirichlet boundary condition u = g on oX. Let f 2 Hk(X), k P 0, g 2 Hs(oX), s P 12 and u is the solution of the original problem. Further, let uh and kh be the approximate solution of the RBF method with Lagrangian multiplier. Assuming that the centers of both X and oX are quasi-uniform with density h1, h2 respectively. The RBF on oX is w 2 Cl . h1 h2l3dþ3 is 2 sufficiently small. Then      ou kþ1 sþ32  ku  uh kH 1 ðXÞ þ kh   6 C h kf k kgk k ðXÞ þ h sþ1 ðoXÞ ; H H on 1 H

2 ðoXÞ

where h = max{h1, h2}, C = min{h1, h2}2l5d+6. 2.2. Discreted by RBF For a given discrete centers X  X and some RBF /(r), we can define the finite dimensional space o n    V ih ¼ span / x  xij  ; j ¼ 1; 2; . . . ; N i ;

754

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

the space spanned by the RBFs with centers in Xi , Ni is the number of centers in Xi , XC the centers on C. Assuming that the density in Xi are equivalent to h1 ¼ supx2X minxi 2X kx  xi k, h2 ¼ supx2C minxi 2X C kx  xi k are the density of the distinct centers X = {x1, x2, . . . , xN} and XC, respectively. kÆk is the Euclidean norm. Kh ¼ spanfð/ðkx  xj kÞÞ; xj 2 X C g. In sake of completeness, we introduce some radial basis function below: 2

expðck  k Þ

Gaussians Multiquadrics

c > 0;

2 b

2

ðc þ k  k Þ

b 2 N ; c 6¼ 0;

b

Thin-plate splines k  k log k  k b Thin-plate splines k  k Sobolev splines

b 2 2N ; b is odd; 2

ekk ð3 þ 3k  k þ 3k  k Þ.

The RBF introduced above are global RBF. One of the most popular classes of CS-RBFs is the one introduced by Wendland [22]. Another is given by Wu [30]. These functions are strictly positive definite in Rd for all d less than or equal   to some fixed value d0, and can be constructed to have any desired amount of smoothness 2j. For l ¼ d2 þ j þ 1 and j = 0, 1, 2, 3, Wendland’s function reads: : l /l;0 ðrÞ¼ð1  rÞþ ; : lþ1 /l;1 ðrÞ¼ð1  rÞþ ½ðl þ 1Þr þ 1; : 2 2 /l;2 ðrÞ¼ð1  rÞlþ2 þ ½ðl þ 4l þ 3Þr þ ð3l þ 6Þr þ 3; : lþ3 /l;3 ðrÞ¼ð1  rÞþ ½ðl3 þ 9l2 þ 32l þ 15Þr3 þ ð6l2 þ 36l þ 45Þr2 þ ð15l þ 45Þr þ 15. Here and throughout, G denotes equality up to a constant factor. For example, the functions : 2 /2;0 ¼ð1  rÞþ ; : /3;1 ðrÞ¼ð1  rÞ4þ ½4r þ 1; : 6 /4;2 ðrÞ¼ð1  rÞþ ½35r2 þ 18r þ 3; : /5;3 ðrÞ¼ð1  rÞ8þ ½32r3 þ 25r2 þ 8r þ 1; which are in C0 ; C2 ; C4 ; C6 , respectively, and strictly positive definite in R3. Wu’s CS-RBF reads: 6

/ðrÞ ¼ ð1  rÞþ ð5r5 þ 30r4 þ 72r3 þ 82r2 þ 36r þ 6Þ 2 C6 . ^ satisfying Assuming that the RBF satisfying / 2 Cl with Fourier transform / ^  ð1 þ ktkÞ2l . /ðtÞ The following interpolation error holds. Lemma 2.1 (RBF-interpolant error estimate theorem [21,23]). Let X  Rd be an open and bounded domain, having Lipschitz boundary. Denote by Su the interpolant on X = {x1, x2, . . . , xN}  X to a function u 2 Hk(X) with k > d/2. Then there exists a constant h0 > 0 such that for all X with h < h0 where h is the density of X, the estimate ku  su kH j ðXÞ 6 Chkj kukH k ðXÞ is valid for 0 6 j 6 k. ^  ð1 þ ktkÞld , then for "g 2 Kh, the following inequality holds Theorem 2.2. Assuming that wðtÞ kgkH 1=2 ðoXÞ 6 Ch22l3dþ3 kgkH 1=2 ðoXÞ ; where C is a constant independent of the density h2.

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

755

Proof. We will denote by C, a generic constant independent of h2. According to the definition of kgkH 1=2 ðoXÞ , we have kgkH 1=2 ðoXÞ ¼

kgk2H 0 ðoXÞ hg; li hg; gi P ¼ . klkH 1=2 ðoXÞ kgkH 1=2 ðoXÞ kgkH 1=2 ðoXÞ

sup l2H 1=2 ðoXÞnf0g

Assuming that g ¼

P

i wðkx

 xi kÞki , by using the skill in [26,28], we have

2 2 kgkH 0 ðoXÞ P Ch22lþ2d2 k~ kk ;

furthermore kgk2H 1 ðoXÞ

2   X   ¼ wðkx  xi kÞki   1  i

6

X

Pn

i¼1 jai jÞ

2

6n

Pn

i¼1 jai j

2

kwðkx  xi kÞki kH 1 ðoXÞ

6

i

H ðoXÞ

by using ð

!2 max kwðkx  xi kÞk2H 1 ðoXÞ i

X

!2 jki j

i

, we have

2 2 2 2 2 kgkH 1 ðoXÞ 6 N max kwðkx  xi kÞkH 1 ðoXÞ k~ kk 6 Chdþ1 max kwðkx  xi kÞkH 1 ðoXÞ k~ kk . 2 i

i

This yield kgk2H 0 ðoXÞ P Ch22lþ3d3 kgk2H 1=2 ðoXÞ . Then kgkH 1=2 ðoXÞ 6 Ch22l3dþ3 kgkH 1=2 ðoXÞ .



Let /i,C be the radial basis function on C, the algorithm reads Step 1. Solve Eq. (3) to obtain ui . Step 2. Assuming that kh ¼

P

j C j /j;C ,

solve Eq. (4) for k = /j,C.

Step 3. SolvePinterface equation (15) to obtain the coefficients Cj, the original solution is u~i ¼ ui þ j C j Hi /j;C . Theorem 2.3. Assuming that the original solution u satisfy that density of discrete centers in X and C, respectively.  1  kþ kui  u~i kH 1 ðXi Þ 6 C h1 2 þ hk1 þ h1 h22l3dþ3 kuk

H

kþ3 2 ðXÞ

ou on

2 H k ðCÞ; k P 1; /C 2 Cl , h1 and h2 are the

.

Proof. Following the analysis of Cea´ theorem, we have kHi kh  Hi kkH 1 ðXi Þ 6 inf kgh  kkH 0 ðCÞ 6 ChkC kkkH k ðCÞ 6 ChkC kuk gh 2Kh

H

kþ3 2 ðXÞ

.

Denoted by Hi kh , the discrete solution, By using the priori estimate [12] and interpolation error estimates, we have      H kh  Hi k 1 6 Hi kh  Hi kh H 1 ðXi Þ þ kHi kh  Hi kkH 1 ðXi Þ 6 Chk2 kuk kþ32 þ Ch1 kHi kh kH 2 ðXi Þ ; i H ðXi Þ H

ðXÞ

756

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

according to the inverse inequality, we have kkh kH 1=2 ðCÞ 6 Ch22l3dþ3 kkh k

H

1 2 ðCÞ

this yield    H kh  Hi k 1 6 Chk2 kuk i H ðXi Þ

H

;

kþ3 2 ðXÞ

þ Ch1 h22l3dþ3 ð1 þ Chk2 Þkuk

H

kþ3 2 ðXÞ

;

then kþ12

kui  u~i kH 1 ðXi Þ 6 Cðh1 This end the proof.

þ hk2 þ h1 h22l3dþ3 Þkuk

H

kþ3 2 ðXÞ

.

h

Remark 2.4. The condition for the density of X and interface C is quite rigorous. This is induced by the inverse inequality (Theorem 2.3). In FEM, the corresponding form reads kgkH 1=2 ðCÞ 6 Ch1 kgkH 1=2 ðCÞ

8g 2 V h ;

where Vh is a suitable finite element space. Nevertheless, these do not effect the numerical computation. 2.3. Numerical examples Throughout this paper, we will denote by R 2 ðuh  ui Þ 2 2 ðL error iÞ ¼ Xi R 2 ; u Xi the relative error of L2-norm in Xi. uh is the numerical solution using Lagrangian multiplier, u is the exact solution of corresponding problem. Example 2.1. Consider Helmholtz equation: ( Du þ u ¼ f in X; ou ¼0 on oX. on

ð9Þ

Let 2

p  p  p þ 1 sin x sin y ; f ðx; yÞ ¼ 2 2 2 the exact solution is p  p  uðx; yÞ ¼ sin x sin y . 2 2 Domain X = [1, 1; 1, 1]. Assuming that X is partitioned into two subdomains X1 = [1, 1; 1, 0], X2 = [1, 1; 0, 1], interface C = [1, 1; 0, 0], Taken 45 points xi = 1 + (i  1) · 0.25, yj = 1 + (j  1) · 0.25, i = 1, . . . , 9, j = 1, . . . , 5 on X1, 45 points xi = 1 + (i  1) · 0.25, yj = (j  1) · 0.25, i = 1, . . . , 9, j = 1, . . . , 5, on X2 and 11 points xi =  1 + (i  1) · 0.2, yj = 0, i = 1, . . . , 9 on C/(r) = r5 as the basis function, Tables 2.1 and 2.2 are the pointwise relative error at discrete centers. Example 2.2. Domain in Example 2.1 is changed into X = [1, 1; 1, 0] [ [0, 1; 0, 1] = X1 [ X2, the equation reads:

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

757

Table 2.1 L2 error 1 = 0.0060 0.0016 0.0002 0.0064 0.0076 0 0.0074 0.0063 0.0002 0.0015

0.0006 0.0000 0.0061 0.0143 0 0.0138 0.0061 0.0000 0.0006

0.0026 0.0025 0.0057 0.0149 0 0.0148 0.0056 0.0025 0.0025

0.0007 0.0052 0.0025 0.0052 0 0.0045 0.0023 0.0051 0.0006

0 0 0 0 0 0 0 0 0

0.0001 0.0002 0.0097 0.0200 0 0.0199 0.0097 0.0001 0.0001

0.0039 0.0019 0.0090 0.0212 0 0.0211 0.0090 0.0019 0.0039

0.0037 0.0047 0.0064 0.0090 0 0.0089 0.0063 0.0047 0.0038

0 0 0 0 0 0 0 0 0

Table 2.2 L2 error 2 = 0.0091 0.0020 0.0011 0.0081 0.0162 0 0.0161 0.0081 0.0011 0.0020

8 Du þ u ¼ f > > > ou < ¼ 0; on  ou > ¼  p2 sin p2 y ; > on >  : ou p ¼ 2 sin p2 x ; on

in X; x ¼ 1; 1

or

y ¼ 1; 1;

ð10Þ

x ¼ 0; y ¼ 0.

Centers in X1 are the same as Example 2.1, centers in X2 are xi = (i  1) · 0.25, yj = (j  1) · 0.25, i = 1, . . . , 5, j = 1, . . . , 5 (Tables 2.3 and 2.4). Example 2.3. Consider 8 Du þ 2u ¼ 0 > > > > ou 1þy > > < on ¼ e ; ou ¼ e1þy ; on > > ou 1þx > ; > > on ¼ e > : ou 1þx ¼ e ; on

the BVP in X; x ¼ 1; x ¼ 1; y ¼ 1;

ð11Þ

y ¼ 1.

Table 2.3 L2 error 1 = 0.0191 0.0130 0.0106 0.0264 0.0558 0 0.0327 0.0125 0.0013 0.0020

0.0070 0.0015 0.0277 0.0452 0 0.0079 0.0043 0.0000 0.0032

0.0151 0.0034 0.0320 0.0441 0 0.0294 0.0019 0.0051 0.0073

0.0025 0.0035 0.0491 0.0488 0 0.0424 0.0037 0.0007 0.0112

0 0 0 0 0 0 0 0 0

758

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

Table 2.4 L2 error 2 = 0.0047 0 0 0 0 0

0 0.0176 0.0023 0.0064 0.0081

0 0.0007 0.0024 0.0018 0.0039

0 0.0080 0.0020 0.0023 0.0006

0 0.0027 0.0044 0.0005 0.0014

Table 2.5 L2 error 1 = 0.0129 0.0017 0.0210 0.0218 0.0023 0.0076 0.0007 0.0161 0.0046 0.0006

0.0068 0.0299 0.0181 0.0109 0.0041 0.0003 0.0154 0.0044 0.0003

0.0173 0.0199 0.0116 0.0296 0.0079 0.0228 0.0096 0.0044 0.0026

0.0171 0.0249 0.0385 0.0183 0.0036 0.0014 0.0223 0.0264 0.0158

0.0106 0.0401 0.0004 0.0550 0.0095 0.0064 0.0048 0.0081 0.0064

0.0407 0.0532 0.0298 0.0463 0.0082 0.0103 0.0170 0.0029 0.0205

0.0110 0.0163 0.0094 0.0095 0.0032 0.0067 0.0114 0.0004 0.0010

0.0125 0.0152 0.0100 0.0179 0.0063 0.0063 0.0026 0.0018 0.0017

0.0019 0.0101 0.0137 0.0080 0.0000 0.0030 0.0075 0.0077 0.0055

0.0002 0.0033 0.0007 0.0101 0.0004 0.0035 0.0013 0.0011 0.0021

0.0082 0.0038 0.0045 0.0048 0.0019 0.0015 0.0029 0.0001 0.0031

Table 2.6 L2 error 2 = 0.0051 0.0006 0.0101 0.0227 0.0003 0.0033 0.0047 0.0139 0.0085 0.0009

The exact solution is u = ex+y. Centers in X1 are xi =  1 + (i  1) · 0.25, yj =  1 + (j  1) · 0.25, i = 1, . . . , 9, j = 1,    , 6. Centers in X2 are xi = 1 + (i  1) · 0.25, yj = (j  1) · 0.25, i = 1, . . . , 9, j = 1,    , 6, centers on C are xi = 1 + (i  1) · 0.2, yj = 0, i = 1, . . . , 11, Taken /(r) = r5 as the basis function, the numerical result are given in Tables 2.5 and 2.6. There are some difficulties if we expect higher accuracy. Which is due to the condition number in solving the interface equation. One way to avoid it is choose the function introduced by Agoshkov [1]. This is not our target because that we need a true meshless method. All of the linear system are solved by GMRES. 3. PDM using RBF collocation methods As that we have pointed out, RBF collocation methods are lack of theoretical analysis. Nevertheless, it perform better in numerical computation. In this section, we will just introduce the algorithm of PDM using RBF collocation and give some examples. 3.1. Algorithm For the general elliptic BVP:  Lu ¼ f in X; u¼0 on oX. Domain X and Xi are the same as which in Section 2.

ð12Þ

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

759

Denoting by k = ujC, the unknown trace of the solution of Eq. (12) on the interface C, ui ¼ ui þ Hi k, with 8  > < Lui ¼ f in Xi ; T ui ¼ 0 on oXi oX; ð13Þ > :  on C. ui ¼ 0 Hi k is the harmonic extension of the trace k 8 > < LHi k ¼ 0 in Xi ; T Hi k ¼ 0 on oXi oX; > : on C. Hi k ¼ k

ð14Þ

An additional equation, named Steklov–Poincare´ equation, should be satisfied to ensure ui þ Hi k is the solution of Eq. (12). Sk ¼ v on C; where v ¼ Sk ¼

ou2 on



ou1 on

ð15Þ

. S is the Steklov–Poincare´ operator,

oH1 k oH2 k  . on on

Algorithm Step 1. Solving Eq. (14) by taking k = /(x  xk), we can obtain Hi /ðx  xk Þ. where xk 2 C. i.e. Solve linear system 8 LHi /ðxj  xk Þ ¼ 0; xj 2 Xi ; > > < T Hi /ðxj  xk Þ ¼ 0; xj 2 oXi oX; ð16Þ > > : Hi /ðxj  xk Þ ¼ /ðxj  xk Þ; xj 2 C. Step 2. Solving Eq. (13) 8  Lu ðxj Þ ¼ f ðxj Þ; > > < i ui ðxj Þ ¼ 0; > > :  ui ðxj Þ ¼ 0;

to get ui ¼ xj 2 X i ; xj 2 oXi

k C k /ðx

 xk Þ, with xk 2 Xi , i.e,

ð17Þ

oX;

xj 2 C.

Step 3. Assuming that Hi k ¼ (15), i.e. Skðxj Þ ¼ vðxj Þ;

T

P

P

j cj Hi /ðx

 xj Þ with xj 2 C, we can get cj by solving the Steklov–Poincare´ Eq.

xj 2 C.

Step 4. The approximated solution is obtained by ui ¼ ui þ 3.2. Numerical examples Throughout this section, we will denote by R ð~ ui  uexact Þ2 2 2 ; ðL ðXi Þerror Þ ¼ XiR 2 ðuexact Þ Xi

ð18Þ P

j cj Hi /ðx

 xj Þ.

760

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

the relative error of L2-norm in Xi. e u i is the numerical solution in Xi, denote by L1 ðXi Þerror ¼ sup Xi

je u i  uexact j ; juexact j

the maximum pointwise relative error at centers. Example 3.1. Consider the Poisson equation  Du ¼ f in X; u¼0

ð19Þ

on oX;

with X = X1 [ X2 = [1, 1; 1, 0] [ [1, 1; 0, 1], C = {(x, y)jx 2 [1, 1], y = 0}, the right-hand side f1 ðx; yÞ ¼ ex ðy 2  1Þðx2 þ 4x þ 1Þ þ 2ex ðx2  1Þ; f2 ðx; yÞ ¼ 4ðx2 þ y 2 Þ sinðx2  1Þ sinðy 2  1Þ þ 2 sinðy 2 þ x2  2Þ; respectively and the corresponding exact solutions are u1 ðx; yÞ ¼ ðx2  1Þðy 2  1Þex ; u2 ðx; yÞ ¼ sinðx2  1Þ sinðy 2  1Þ. Take TPS /(r) = r6 log (r) and Sobolev spline /(r) = e r(3 + 3r + 3r2) as the RBF. The discrete centers is uniformly distributed in X with density h. Tables 3.1 and 3.2 are the numerical results. Example 3.2. Consider Poisson equation  Du ¼ f in X; u¼0

ð20Þ

on oX;

with X = X1 [ X2 = [1, 1; 1, 0] [ {(x, y)jy P 0, x2 + y2 6 1}, the discrete centers are taken as Fig. 1. Assuming that f ðx; yÞ ¼ cosðx2 þ y 2  1Þ½ð10x2  2Þðy þ 1Þ þ ð6y þ 2Þðx2  1Þ  sinðx2 þ y 2  1Þ½ð4y 3 þ 4y 2 Þðx2  1Þ þ ð4x4  4x2  2Þðy þ 1Þ; the exact solution of this problem is (see Table 3.3) u ¼ ðx2  1Þðy þ 1Þ sinðx2 þ y 2  1Þ. Example 3.3. Consider the groundwater flow problem in a heterogeneous and anisotropic medium: (

Þ þ oyo ðbðx; yÞ ouðx;yÞ Þ ¼ f ðx; yÞ in X; Lu ¼ oxo ðaðx; yÞ ouðx;yÞ ox oy u¼0 on oX;

ð21Þ

Table 3.1 Numerical results for Poisson equation (TPS) f1(x, y) 2

L (X1)error = 0.0082 L1(X1)error = 0.0463

f2(x, y) 2

L (X2)error = 0.0082 L1(X2)error = 0.0466

L2(X1)error = 0.0053 L1(X2)error = 0.0141

L2(X2)error = 0.0053 L1(X1)error = 0.0142

Table 3.2 Numerical results for Poisson equation (Sobolev spline) f1(x, y) 2

L (X1)error = 0.0071 L1(X1)error = 0.0320

f2(x, y) 2

L (X2)error = 0.0071 L1(X2)error = 0.0332

L2(X1)error = 0.0027 L1(X2)error = 0.0088

L2(X2)error = 0.0027 L1(X1)error = 0.0093

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

761

1

0.8

0.6

0.4

0.2

0

–0.2

–0.4

–0.6

–0.8

–1 –1

–0.8

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

0.8

1

Fig. 1. The distribution of distinct centers for Example 3.2. Table 3.3 Numerical result for Example 3.2 (TPS) L2(X1)error = 0.0038

L2(X2)error = 0.0121

Table 3.4 Numerical results for groundwater flow equation (TPS) h = 0.2 2

L (X1)error = 0.0412 L1(X1)error = 0.0714

h = 0.1 2

L (X2)error = 0.0356 L1(X2)error = 0.0658

L2(X1)error = 0.0044 L1(X2)error = 0.0226

L2(X2)error = 0.0037 L1(X1)error = 0.0217

with X = X1 [ X2 = [1, 1; 1, 0] [ [1, 1; 0, 1]. Assuming that aðx; yÞ ¼ 2  x2  y 2 ; bðx; yÞ ¼ exy ; the exact solution uðx; yÞ ¼ xyð1  x2 Þð1  y 2 Þ; the corresponding right side function f ðx; yÞ ¼ exy xð1  x2 Þð3y 2  6y  1Þ þ yð1  y 2 Þð12x3 þ 6xy 2  14xÞ. Table 3.4 is the numerical result of Example 3.3 for uniformly discrete centers with density h. 4. Conclusions We can also use others RBFs, for example, multi-quadric function, Gaussian basis and compactly supported RBFs (CS-RBF). We can obtain a sparse stiff matrix if we use CS-RBF as the basis. The condition number of the linear system will also be better than global RBF. But numerical integration is complex. For

762

Y. Duan / Applied Mathematics and Computation 179 (2006) 750–762

multi-quadric function and Gaussian basis, the numerical result is very sensitive to the control parameter c. Compared with the PDM using FEM or finite difference, the PDM using RBF does not need the pre-construction of the basis on the interface C and the mesh-generation for FEM at the cost that we will solve a linear system with large condition number, which will be larger when the density h of the discrete centers become smaller. As that we have pointed out, it may be something inevitable for the meshless method using RBF for its uncertainty property. References [1] V.I. Agoshkov, E. Ovtchinnikov, Projection decomposition method, Math. Models Methods Appl. Sci. 4 (1994) 773–794. [2] I. Babusˇka, The finite element method with Penalty, Math. Comput. 27 (1973) 221–228. [3] I. Babusˇka, U. Banerjee, J.E. Osborn, Survey of meshless and generalized finite element methods: a unified approach, Acta Numer. (2003) 1–125. [4] T. Belytschko, Y.Y. Lu, L. Gu, P. Krysl, Element-free Galerkin Methods, Int. J. Numer. Methods Eng. 37 (1994) 229–256. [5] T. Belytschko, Y. Krongauze, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Methods Appl. M. 139 (1996) 3–47 (special issue on Meshless Methods). [6] M.D. Buhmann, Radial basis functions, Acta. Numer. (2000) 1–38. [7] M.D. Buhmann, Radial Basis Functions: Theory and Implementations, Cambridge University Press, 2003. [8] G.E. Fasshauer, Solving differential equations with radial basis functions: multilevel methods and smoothing, Adv. Comput. Math. 11 (1999) 139–159. [9] C. Franke, R. Schaback, Convergence order estimates of meshless collocation methods using radial basis functions, Adv. Comput. Math. 8 (1998) 381–399. [10] P. Gervasio, E. Ovtchinnikov, A. Quarteroni, The spectral projection decomposition method for elliptic equations in two dimensions, SIAM J. Numer. Anal. 34 (1997) 1616–1639. [11] M. Griebel, M.A. Schweitzer, A particle-partition of unity method for the solution of elliptic, parabolic and hyperbolic PDEs, SIAM J. Sci. Comput. 22 (2000) 853–890. [12] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, 1985. [13] G.R. Liu, Mesh Free Methods: Moving beyond the Finite Element Method, CRC Press, 2003. [14] J. Li, Y.C. Hon, Domain decomposition for radial basis meshless methods, Numer. Methods PDEs 20 (2004) 450–462. [15] B.N. Khoromskij, J.M. Melenk, Boundary concentrated finite element methods, SIAM J. Numer. Anal. 41 (2003) 1–36. [16] E. Ovtchinnikov, The construction of a well-conditioned basis for the projection decomposition method, Calcolo 30 (1993) 255–271. [17] E. Ovtchinnikov, Projection decomposition method with a well-conditioned basis for the generalised Stokes problem, in: V.K.W. Morton, M.J. Baines (Eds.), Numerical Methods for Fluid Dynamics, Clarendon Press, Oxford, 1995, pp. 515–521. [18] A. Quarteroni, A. Valli, Domain decomposition method for partial differential equation, Oxford University Press, England, 1999. [19] R. Schaback, Error estimates and condition numbers for radial basis function interpolation, Adv. Comput. Math 3 (1995) 251–264. [20] R. Schaback, Improved error bounds for scattered data interpolation by radial basis functions, Math. Comput. 68 (1999) 201–216. [21] H. Wendland, Sobolev-type error estimates for interpolation by radial basis functions, in: A. LeMe´haute´, C. Rabut, L.L. Schumaker (Eds.), Surface Fitting and Multiresolution Methods, Vanderbilt University Press, Nashville, 1997, pp. 337–344. [22] H. Wendland, Error estimates for interpolation by compactly supported radial basis functions of minimal degree, J. Approx. Theory 93 (1998) 258–272. [23] H. Wendland, Meshless Galerkin methods using radial basis functions, Math. Comput. 68 (1999) 1521–1531. [24] X. Zhou, Y.C. Hon, Jichun Li, Overlapping domain decomposition method by radial basis functions, Appl. Numer. Math. 44 (2003) 241–255. [25] Y. Duan, Y.J. Tan, Meshless collocation method based on Dirichlet–Neumann substructure iteration, Appl. Math. Comput. 162 (2005) 317–327. [26] Y. Duan, Y.J. Tan, Meshless Galerkin method based on regions partitioned into subdomains using global RBF, in: Proceedings of ICNAAM 2004, Wiley-VCH, 2004, pp. 128–133. [27] Y. Duan, Y.J. Tan, A Meshless Galerkin methods for Dirichlet problems using radial basis functions, J. Comput. Appl. Math., in press. [28] Y.C. Hon, Z. Wu, Additive Schwarz domain decomposition with radial basis approximation, Int. J. Appl. Math. 4 (1) (2000) 81–98. [29] Z. Wu, R. Schaback, Local error estimates for radial basis function interpolation of scattered data, IMA J. Numer. Anal. 13 (1993) 13–27. [30] Z. Wu, Compactly supported positive definite radial functions, Adv. Comput. Math. 4 (1995) 283–292.