Neumann problems of 2D Laplace's equation by method of fundamental solutions

Applied Numerical Mathematics 119 (2017) 126–145

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Neumann problems of 2D Laplace’s equation by method of fundamental solutions Zi-Cai Li a , Ming-Gong Lee b , Hung-Tsai Huang c,∗,1 , John Y. Chiang d,e,∗ a

National Sun Yat-sen University, Kaohsiung, 80424, Taiwan Department of Leisure and Recreation Management, Ph.D. Program in Engineering Science, Chung Hua University, Hsin-Chu, 30012, Taiwan Department of Financial and Computational Mathematics, I-Shou University, Kaohsiung, 84001, Taiwan d Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung, 80424, Taiwan e Department of Healthcare Administration and Medical Informatics, Kaohsiung Medical University, Kaohsiung, 80708, Taiwan b c

a r t i c l e

i n f o

Article history: Received 4 August 2009 Received in revised form 8 February 2017 Accepted 13 April 2017 Available online 20 April 2017 Keywords: Stability analysis Error analysis Condition number Method of fundamental solutions Neumann problems Truncated singular value decomposition

a b s t r a c t The method of fundamental solutions (MFS) was first used by Kupradze in 1963 [21]. Since then, there have appeared numerous reports of the MFS. Most of the existing analysis for the MFS are confined to Dirichlet problems on disk domains. It seems to exist no analysis for Neumann problems. This paper is devoted to Neumann problems in non-disk domains, and the new stability analysis and the error analysis are made. The bounds for both condition numbers and errors are derived in detail. The optimal convergence rates in L 2 and H 1 norms in S are achieved, and the condition number grows exponentially as the number of fundamental functions increases. To reduce the huge condition numbers, the truncated singular value decomposition (TSVD) may be solicited. Numerical experiments are provided to support the analysis made. The analysis for Neumann problems in this paper is intriguing due to its distinct features. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction To solve elliptic boundary value problems, the fundamental solutions (FS) satisfying the homogeneous equations are chosen, and their linear combination is used to satisfy the boundary conditions. To avoid the logarithmic singularity, the source nodes of FS are located outside of the solution domain S. This numerical algorithm is called the method of fundamental solutions (MFS). The MFS was first used by Kupradze in 1963 [21]. Since then, there have appeared numerous reports of the MFS for computation. The review of the MFS is given in Fairweather and Karageorghis [11], Golberg and Chen [12], and Chen et al. [3]. To celebrate the progress of the MFS, there held the first International Workshop on the Method of Fundamental Solutions (MFS 2007) (see [5]), Ayia Napa, Cyprus, June 11–13, 2007, and the second International Workshop held jointly with International Workshop on Trefftz method IV, National Sun Yat-sen University, Kaohsiung, Taiwan, March 15–18, 2011. Furthermore, the third and the fourth International Workshops on the MFS follow in HangZhou, China, October 11–13, 2015, and in Poznan, Poland, July 4–9, 2017, respectively. On the other hand, the Trefftz method (TM) [46] as boundary method has been fully developed in both theory and computation for several decades (see Li et al. [24,31,35]). In fact, the MFS is one of the TM using fundamental solutions (FS). In [24,35], only the particular solutions (PS) are used, and the method of

* 1

Corresponding authors. E-mail addresses: [email protected] (Z.-C. Li), [email protected] (H.-T. Huang), [email protected] (J.Y. Chiang). The work was supported in part by Ministry of Science and Technology, Taiwan under grant MOST 103-2632-M-214-001-MY3.

http://dx.doi.org/10.1016/j.apnum.2017.04.004 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

127

particular solutions (MPS) may be called. Both the MFS and the MPS are classified to TM, and they can be carried out by the collocation TM (CTM) in [35]. However, only a few papers of the MFS are involved in analysis. Even for the analysis, most of them are confined to the Dirichlet problems on the disk domains (see Kitagawa [20], Christiansen [6–8], Smyrlis and Karageorghis [44,45], and Li [22,23]). The analysis of the MFS in non-disk domains is reported in Bogomolny [2], Mathon and Johnston [40], Katsurada [16–18] and Li et al. [26–28]. So far, there seems to exist no analysis on Neumann problems by the MFS. Since the Neumann problems satisfying the consistency conditions have multiple solutions, and the stability analysis of the MFS is particularly interesting. Hence this paper is devoted to the Neumann problem in non-disk domains, to provide new error and stability analyses. The bounds of both errors and condition numbers are derived in detail. The optimal convergence rates in L 2 and H 1 norms in S are achieved, but the condition number grows exponentially as the number of fundamental functions increases. By the circulant matrix in Davis [9], the eigenvalues of the stiffness matrix from the MFS can be derived exactly for the Neumann boundary condition under the simple case that the source and the collocation nodes are uniformly located on concentric circles. When the stiffness matrix is singular or near singular, the numerical solution can be found by the truncated singular value decomposition (TSVD), see Hansen [14] and Li et al. [29]. In fact, the TSVD was used in Ramachandran [42] to remove the ill-conditioning of the MFS. In this paper, we adopt the TSVD to reduce the huge condition number √ Cond of the MFS down to O ( Cond), see (6.12). Next, we consider the Neumann problems for bounded simply-connected domains, and derive the bounds of the condition number Cond by following Li et al. [28]. Moreover, the error bounds are derived, by following Bogomolny [2] and Li [26], to achieve the optimal convergence rates. Numerical experiments are also carried out to verify the analysis made. The error and stability analyses of TM using other particular solutions, such as harmonic polynomials, have been explored systematically in [24,31,35]. Our recent efforts are devoted to develop the analysis of the MFS, to fill the gap between the advanced computation and the existing theory. A number of new and important results of accuracy and stability have been obtained. This paper reports one of them; other results can be found in [25–28,30,33,36,37,39,48]. This paper is organized as follows. In the next section, the algorithms of the MFS are described, the main results of the analysis in this paper are briefly addressed, and some guidance is provided for engineering applications. For the Neumann problems, the bounds of Cond are derived for disk domains in Section 3, and for the bounded simply-connected domains in Section 4. In Section 5, error analysis is made. In Section 6, numerical experiments are carried out, and in the last section, a few concluding remarks are made. 2. Method of fundamental solutions 2.1. Description of algorithms Consider Laplace’s equation with the Neumann boundary condition in 2D,

u =

∂ 2u ∂ 2u + 2 = 0 in S , ∂ x2 ∂y ∂u = g on , ∂ν

(2.1) (2.2)

where S is a bounded simply-connected domain with the boundary  = ∂ S, and ν is the exterior normal of  . For the existence of the solution of (2.1) and (2.2), we need the consistency condition  g = 0, so that the solution plus any constant is also a solution.  is smooth without corners of the boundary  (see Fig. 1). The assumption of g is described in Section 2.2 as

g ∈ H q (), q ≥

1 2

(2.3)

.

Denote in Fig. 1,

rmax = max r , S

rmin =

max

S in ( S in ⊆ S )

r,

(2.4)

where S in is the inscribed disk of S. Let the source (i.e., charge) nodes Q be located outside of S, then the fundamental solutions,

φ(r , θ) = ln | P Q |,

P ∈ S ∪ ∂ S,

(2.5)

are harmonic in S, where

P = {(x, y ) | x = r cos θ, y = r sin θ}.

(2.6)

A circle surrounding S is given by

 R = {(r , θ)| r = R ,

0 ≤ θ ≤ 2π },

R > rmax .

(2.7)

128

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

Fig. 1. The solution domain S with its circles.

Based on Bogomolny [2], the source nodes Q i may be simply located uniformly on  R :

Q i = {(x, y ) | x = R cos ih, y = R sin ih}, where R > rmax and h =

2π N

φi ( P ) = ln | P Q i |,

(2.8)

. We obtain the fundamental solutions (FS),

i = 1, 2, · · · , N ,

(2.9)

and represent the numerical solution of (2.1) and (2.2) as the linear combination

uN = c +

N 

c i φi ( P ),

(2.10)

i =1

where c i are the unknown coefficients to be sought, and c is an arbitrary constant. Since u N satisfies Laplace’s equation in S already, the coefficients c i can be sought by satisfying the boundary condition (2.2) only. We follow the Trefftz method (TM) in [24,35],2 to seek u N (i.e., c i ). Denote the energy by



I (u ) =

(

∂u − g )2 . ∂ν

(2.11)



Also denote V N as the set of (2.10). Then the numerical solution u N can be sought by

I (u N ) = min I ( v ).

(2.12)

v ∈V N

When the integral in (2.11) involves approximation, denote

 I (v ) = 



( 

∂v − g )2 , ∂ν

(2.13)



where  is the numerical approximation of  by some quadrature rules, such as the Gaussian rule or the trapezoidal rule. Hence, the numerical solution u˜ N ∈ V N is obtained by

 I (u˜ N ) = min  I ( v ).

(2.14)

v ∈V N

We may establish the collocation equations directly from (2.2), to yield N   ∂  ci φi ( P j ) = g ( P j ), ∂ν

P j ∈ .

i =1

2

In [24], for the Trefftz method, the boundary method or the boundary approximation method is called.

(2.15)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

On the other hand, let  be divided into small  j with the mesh spacings h j , i.e.,  =

βj

N

∂  ci φi ( P j ) = β j g ( P j ), ∂ν

129

M

j =1 

P j ∈ ,

i =1

j

. We obtain from (2.15)

(2.16)

√

where P j are the Gaussian nodes, the weights β j = O ( h), and h = max j h j . In computation, we may choose the number of collocation nodes to be equal to or larger than that of source nodes. The equations (2.16) are called the collocation Trefftz method (CTM) in [35]. The equations (2.16) are equivalent to (2.14), where the Gaussian rule or the trapezoidal  rule is chosen for  . This displays that the CTM is the TM involving integration approximation, and that the MFS is just the Trefftz methods in [35] by using FS. The overdetermined system is obtained from (2.16)

Fx = b,

(2.17)

m×n

where F ∈ R (m ≥ n), x ∈ R and b ∈ R . In this paper, we consider the rank deficient with rank(F) = t ≤ n. The singular value decomposition (SVD) of F is given by F = U V T , where matrices U ∈ R m×m and V ∈ R n×n are orthogonal, and matrix  ∈ R m×n is diagonal. The non-zero diagonal entries of  are given by t positive singular values: n

m

σ1 ≥ ... ≥ σt > 0, t ≤ n.

(2.18)

Denote U = (u1 , ..., um ) and V = (v1 , ..., vn ), where ui ∈ R

x=

t  uT b i

i =1

σi

m

and vi ∈ R are vectors. The solution of (2.17) is obtained by n

vi .

(2.19)

The modified condition number for (2.19) is defined in van Loan [47], by

σ1  Cond = .

(2.20)

σt

When t < n, the solution from (2.19) is called the TSVD solution. To distinguish the condition number with rank deficient from the traditional condition number

Cond =

σ1 , σn > 0, σn

(2.21)

we use the notation  Cond with a hat. Even for rank(F) = n, when σt +1 , ..., σn are infinitesimal, the TSVD solution (2.19) is solicited for better stability, see Hansen [14], Ramachandran [42] and Li et al. [29]. 2.2. Main results of the analysis and their applications Since the FS in (2.9) with Q i located outside of S are highly smooth, the functions as (2.10) are not well-suitable to approximate the singular solutions, where the derivatives are unbounded. In the error analysis in Section 5, we confine the smooth problems for the smooth boundary  without corners, and assume that the solution in (2.1) and (2.2) satisfy

u ∈ H p ( S )/ P 0 ( S ), p ≥ 2, p

(2.22) p

where H ( S )/ P 0 ( S ) is a quotient space, H ( S ) are the Sobolev spaces, and P 0 ( S ) is a constant space. The equation (2.22) implies

g ∈ H q (), q ≥

1 2

(2.23)

.

Condition (2.23) may be relaxed as g ∈ H q (), q > 0 for the solutions involving corner singularity by the MFS. Note that 1

for the interior Neumann problems, the solutions must satisfy u ∈ H 1 ( S )/ P 0 ( S ) and g ∈ H − 2 () from Babuska and Aziz [1]. Under (2.22), the optimal convergence in L 2 and H 1 norms in S can be achieved, and given in Theorem 5.1 as

u − u N k, S = O (

1 N p −k

) u p , S , k = 0, 1.

(2.24)

The above errors in the quotient space H k ( S )/ P 0 ( S ) mean that the differences of a constant solution are not included into (2.24). The stability is more severe for the MFS. Under the source nodes given in (2.8), the following bounds of Cond are given in Corollary 4.2 as

130

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145



Cond(F) =

 Cond(F) =

σmax (F) R ≤C σmin (F) rmin σmax (F) σmin-next (F)



≤C

N (2.25)

, R

rmin

N 2

(2.26)

,

where rmin is given in (2.4), and C is a bounded constant independent of N. From (2.24), the higher degree of smoothness the solution u has (i.e., the higher p is), the faster the convergence is reached. For the infinitely smooth solution with p = ∞, Eq. (2.24) indicates an exponential convergence. To deal with singularity problems, some renovations should be developed in the MFS, by following [24]. First, the local refinements of source nodes near the corners are adopted as in [27]. Next, a few leading singular particular solutions (SPS) in [34] may be added into the FS as in [27]. Numerical results for other problems in [28,32,39] show that this adding technique is very promising. Furthermore, we can employ the combined techniques of [24], where the SPS are chosen in the subdomains of S involving singularity, and the FS in the rest of S. The analysis and computation are explored in [25,36], to display a high efficiency to solve singularity problems with corners and the point at infinity in unbounded domains. Better choices of source f nodes Q i are important to the MFS. The choice (2.8) is simplest, and it is optimal for the disk domains. When S is arbitrary in shapes, other kinds of exterior contours of source nodes are worthy to be evaluated. First, the conformal transformation from a circle is suggested in Katsurada and Okamoto [19], to yield a non-circular contour of source nodes. Next, the greedy adaptive techniques of Schaback [43] may select good locations of source nodes. For the singularity problem of bihamonic equations, numerical experiments in [33] display an advantage in applications. For arbitrary solution domains, it is reasonable to locate the source nodes Q i on a larger fictitious contour to  with an equal distance. Several approaches are proposed in Chen et al. [4]. Furthermore, the selection of optimal location of Q i remains an outstanding open problem. Here, let us outline the key ideas of analysis, to help readers to easily follow the proof. Under (2.8), the eigenvalues of circulant matrices are known from the Davis theory [10] in Lemma 3.1, and the errors of the trapezoidal rule can be found from Lemmas 3.2 and 3.3. Hence, the eigenvalues of matrix A in (3.3) can be derived for disk domains in Section 3. Based on the matrix and numerical analysis, the eigenvalues of matrix F in (2.17) can also be derived for the simply-connected domains in Section 4. Note that the bounds in (2.25) and (2.26) are exponential via N. Since the MFS can be classified into the Trefftz method, we may derive the error bounds by following [24,35]. From (5.19) and (5.21), the errors of the MFS can be split into two parts: (a) the remainder R n of the harmonic polynomials P n , and (b) the differences between P n and a special solution u ∗N of the MFS. The bounds of (b) can also be derived by means of Lemmas 3.2 and 3.3, where Eq. (2.8) plays a key role, see [26]. The optimal convergence means that the error orders of the MFS solutions are the same as those of the best harmonic polynomials [35,38]. However, since the exponential Cond may damage the accuracy of the numerical solutions, to reduce the Cond is significant to the MFS in applications. 3. Stability analysis on disk domains Consider the Neumann problem (2.1) and (2.2) on disk domains with radius nodes, Eq. (2.15) may lead to N  i =1

ci

∂ φi ( P j ) = g ( P j ), ∂r

j = 1, 2, · · · , N ,

ρ . For the uniformly distributed collocation

(3.1)

where φi ( P j ) = ln | Q i P j |, Q i are given in (2.8), and the collocation nodes

P i = {ρ cos ih,

ρ sin ih}, i = 1, 2, · · · , N ,

(3.2)

with R > ρ . The equations (3.1) can also be denoted by

Ax = b,

(3.3)

where x = (c 1 , · · · , cn ) , b = ( g ( P 1 ), · · · , g ( P N )) , and A ∈ R T

T

N ×N

∂ ∂ 1 −−→ φi ( P j ) = cos(r , Q i P j ), ln | Q i P j | = ∂r ∂r |Q i P j|

. Below, we give the explicit entries of matrix A. We have

(3.4)

where

−−→ r = O P j = ρ cos jh  + ρ sin jh m , −−→ , Q i P j = (ρ cos jh − R cos ih) + (ρ sin jh − R sin ih)m

(3.5) (3.6)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

131

 are the unit vectors along x and y directions, respectively. Moreover, the equality exists: h = 2Nπ , and  and m

−−→

cos(r , Q i P j ) =

= =

1

ρ|Q i P j| 1

−−→ (r , Q i P j ) −−→ |r | | Q i P j |

{ρ cos jh(ρ cos jh − R cos ih) + ρ sin jh(ρ sin jh − R sin ih)}

ρ|Q i P j|



ρ 2 − ρ R cos( j − i )h .

(3.7)

Then the derivative (3.4) leads to

∂ ρ R cos( j − i )h − ρ 2  φi ( P j ) = − 2 ∂r ρ R + ρ 2 − 2R ρ cos( j − i )h =−

a cos( j − i )h − 1

ρ

(a2

+ 1 − 2a cos( j − i )h)

(3.8)

,

where a = ρR > 1. Hence, the matrix A is also the circulant matrix:

A = Circulant(a0 , a1 , · · · , a N −1 ),

(3.9)

where

a0 =

1

ρ (1 − a)

ak = −

(3.10)

,

a cos kh − 1

k = 1, 2, · · · , N − 1.

,

ρ (1 + a2 − 2a cos kh)

(3.11)

Since ai = a N −i , the matrix A is symmetric. We cite a lemma from Davis [9]. Lemma 3.1. Let matrix A = Circulant(a0 , a1 , · · · , a N −1 ) with the entries ai in (3.10) and (3.11). Then the eigenvalues of the symmetric matrix A are given by

λk =

N −1 

a j cos(kjh),

k = 0, 1, · · · , N − 1,

(3.12)

j =0

where h = 2Nπ . Also, the vector y0 = √1 (1, 1, · · · , 1) T is the eigenvector corresponding to the leading eigenvalue λ0 . N

To express the eigenvalues λk in (3.12) as O ( N μ ) or O (ak ), we will apply the following lemma from Davis and Rabinowitz [10] and Bogomolny [2]. Lemma 3.2. For the periodical function g = g (x) on [0, 2π ], denote the error by

2π N ( g ) = T N ( g ) −

g (x)dx,

(3.13)

0

where the trapezoidal rule, T N ( g ) = h

N −1 k=0

g (kh) with h = 2Nπ , is used. Then the equalities exist,

 N (sin mx) = 0, m = 1, 2, ...,  2π , if m = ν N , ν = 1, 2, ...,  N (cos mx) = 0, otherwise.

(3.14) (3.15)

We also cite a lemma from Li et al. [28]. Lemma 3.3. Let g (x) ∈ C ∞ [0, 2π ] satisfy the periodic boundary conditions,

g () (0) = g () (2π ),  = 0, 1, ...

(3.16)

132

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

Then the errors exist:

2π T N ( g (x) cos(kx)) −

g (x) cos(kx)dx

(3.17)

0

=π

∞ 

{αν N −k + αν N +k }, k = 0, 1, ..., N − 1,

ν =1

where

α =

2π

1

g (x) cos(x)dx,  = 0, 1, ...

π

(3.18)

0

Theorem 3.1. For the Neumann problem (2.1) and (2.2) on the disk S, the eigenvalues of matrix A from the MFS in (3.1) are given by

λ0 =

−N 1 × N, − N ρ (1 − a ) a

λk = −

N

{

1

2ρ ak

1

+

(3.19) 1

(

1 − a− N a N −k

+

1 a N +k

)},

k = 1, 2, · · · N − 1,

(3.20)

where a = ρR > 1. Proof. Form Lemma 3.1, the eigenvalues of the circulant matrix A are given by

λk =

N −1 

a j cos kjh,

k = 0, 1 , · · · , N − 1 ,

(3.21)

j =0

where a j are given in (3.10) and (3.11). First, we have from Lemma 3.3

λ0 =

N −1 

aj = −

j =0

=−

1

ρh

⎡ 2π  ⎣ 0

N −1  j =0

a cos jh − 1

(3.22)

ρ (1 + a2 − 2a cos jh) ⎤

a cos x − 1 1 + a2 − 2a cos x

dx⎦ + λ0 ,

where the explicit λ0 is given in (3.27) below. Since

a cos x − 1 1 + a2 − 2a cos x

1

a2 − 1

1

2

2

1 + a2 − 2a cos x

=− +

,

the equation (3.22) leads to

⎧ ⎫ 2π ⎨ ⎬ 2 1 a −1 1 dx + λ0 . λ0 = π− ρh ⎩ 2 1 + a2 − 2a cos x ⎭

(3.23)

0

There exists the integration formula [13, p. 366],

2π 0

⎧ ⎪ 2π ak ⎪ ⎨ , a 2 < 1, cos kx 1 − a2 dx = 2π 1 ⎪ 1 + a2 − 2a cos x ⎪ , a 2 > 1. ⎩ 2 (a − 1) ak

(3.24)

Hence we have from (3.23) and (3.24),

λ0 =

1

ρh



π−

a2 − 1 2

·

2π a2 − 1

 + λ0 = λ0 .

(3.25)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

133

Now, we show λ0 . From (3.24), the expansion exists,

1

=

1 + a2 − 2a cos x

¯ = where a2 > 1 and α

λ0 = −

2

ρ

aν N

ν =1

2

∞ 

{α¯  cos x + β¯ sin x},

(3.26)

=1

1

− 1 a

× 2π

( ≥ 1). Also from Lemma 3.3 and (3.23), we obtain ∞ 

α¯ ν N

(3.27)

ν =1

∞ N 1

=−

+

2 a2

1 a2 − 1

ρh

α¯ 0

=

−N 1 × N. − N ρ (1 − a ) a

Combining (3.25) and (3.27) gives the first desired result (3.19). Next, we have from (3.21) and Lemma 3.3

λk = −

N −1  j =0

=−

cos kjh(a cos jh − 1)

2π

1

(3.28)

ρ (1 + a2 − 2a cos jh) cos kx(a cos x − 1) 1 + a2 − 2a cos x

ρh 0

dx + λk , k = 1, 2, ..., N − 1,

where the explicit λk is given in (3.30) below. Since

1

cos kx cos x =

2

(cos(k + 1)x + cos(k − 1)x),

we have from (3.28) and (3.24),

⎧ ⎫ 2π 2π ⎨ ⎬ 1 a cos(k + 1)x + cos(k − 1)x cos kx dx + dx + λk − λk = 2 2 ρh ⎩ 2 1 + a − 2a cos x 1 + a − 2 cos x ⎭ 0 0    1 a 2π 1 1 2π 1 − 2 + λk = + k −1 + 2 ρh 2 a − 1 ak+1 a − 1 ak a =−

1

π

ρh

ak

+ λk = −

N 2ρ ak

+ λk .

(3.29)

From Lemma 3.3, we obtain from (3.26) and (3.29)

λk =

1

ρh +π

−

∞ aπ 

∞ 

2

(α¯ ν N −k + α¯ ν N +k )

ν =1

=

2π

1

ρ h a2 − 1 +

∞ 

(α¯ ν N −(k+1) + α¯ ν N +(k+1) + α¯ ν N −(k−1) + α¯ ν N +(k−1) )

ν =1



−

∞ a

2

1

1

1

1

( ν N −(k+1) + ν N +(k+1) + ν N −(k−1) + ν N +(k−1) ) a a a a

ν =1

1 1 ( ν N −k + ν N +k ) a a

ν =1

∞

=−

π  1 1 N 1 1 { + }=− + N +k }. { ρ h ν =1 aν N −k aν N +k 2ρ (1 − a− N ) a N −k a

Combining (3.29) and (3.30) gives desired result (3.20), and this completes the proof of Theorem 3.1.

(3.30)

2

In fact, we have used Mathematica for computing λk for three different cases: (1) Eigenvalues of matrix A in (3.9), (2) Eq. (3.12), and (3) Eqs. (3.19) and (3.20). The computed eigenvalues for three cases are perfectly matched with each other.

134

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

From (3.20) we have for k ≥ 1

|λk | =

N

{

1

2ρ ak

+

1

(

1

1 − a− N a N −k

+

1 a N +k

)} ≤ |λ1 | ≤ C N .

(3.31)

Let N be even without loss of generality, we then have N

|λk | ≥ |λ N | ≥ c 0 Na− 2 , k = 1, 2, ..., N ,

(3.32)

2

where C and c 0 (> 0) are two constants independent of N. In the following, C and c 0 are always the constants independent of N, but their values may be different in different contexts. Note that from (3.19) and (3.20) we have mink |λk | = |λ0 | > 0, to indicate the non-singularity of matrix A in (3.3) from the MFS. This result is interesting for Neumann problems. However, since

|λ0 | ≈

N 1

ρ aN

=

ρ ( )N ρ R

N

(3.33)

is much smaller than others in (3.32), |λ0 |  |λk | (k ≥ 1), for better stability we should employ the TSVD solution (2.19). By ignoring the components with respect to u0 (= y0 ), we have

x=

N −1 

uiT b

σi

i =1

vi ,

(3.34)

where t = N − 1 and n = N in (2.19). In this case the modified condition number,

 Cond =

max |λi | min-next|λi |

(3.35)

,

is given from (2.20), where min-next|λi | > 0 and min |λi | = 0. Then we have the following corollary from (3.31), (3.32) and (3.33). Corollary 3.1. Let the conditions in Theorem 3.1 hold. Then the bound exists,

Cond =

|λ1 | ≤C |λ0 |

 N R

ρ

(3.36)

.

Also by the TSVD in (3.34), the bound exists,

|λ1 |  Cond = ≤C |λ N | 2

 N R

2

(3.37)

.

ρ

4. Stability analysis for bounded simply-connected domains From Section 2, the MFS is classified to the Trefftz method by choosing the fundamental solutions as admissible functions. Hence, we may study the MFS as the TM [46], and the TM involving integration approximation (i.e. the CTM) in [35]. 4.1. Trefftz methods Consider the Neumann problem (2.1) and (2.2), where S is a bounded simply-connected domain. First, we study the stability of the TM (2.12) without integration approximation, and obtain the normal equation

Bx = b,

(4.1)

where x = (c 1 , · · · , c N ) , b ∈ R by T

1 2

N

is the known vector, and B ∈ R

N ×N

. The matrix B is symmetric and positive definite, given

x T Bx = I 0 ( v ),

where v ∈ V N and



I0(v ) =

( 

∂v 2 ) . ∂ν

(4.2)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

135

Hence, the maximal and the minimal eigenvalues of B can be obtained by, respectively

x T Bx

λmax (B) = max x=0

2I 0 ( v )

= max

xT x

x 2

x=0 v ∈V N

,

(4.3)

.

(4.4)

and

λmin (B) = min

x T Bx

x=0

xT x

= min

2I 0 ( v )

x=0 v ∈V N

x 2 ∂φ

When R > rmax , since the derivatives ∂ νi ( P ) are bounded on  , we have from the Schwarz inequality,



∂v ( )2 ds = ∂ν



⎧ ⎫ 2   N N N  ⎨  2 ⎬  ∂φ ∂φi  i  ci (s) ds ≤ c i2 (s) ds  ⎩ ⎭ ∂ν ∂ν i =1



N 

≤ CN

i =1

 i =1

c i2 ≤ C N x 2 .

(4.5)

i =1

From (4.2) and (4.3), we obtain



λmax (B) = max x=0 v ∈V N

2  ( ∂∂ νv )2 ≤ C N. x 2

(4.6)

The challenging work is to estimate the lower bound of λmin (B) in (4.4). Since the matrix B is nonsingular from Theorem 3.1, we have λmin (B) > 0. Denote the inscribed disk S 0 ⊆ S and its circular boundary 0 by

S 0 = {(r , θ)| r ≤ rmin , 0 ≤ θ ≤ 2π } ,

0 = {(r , θ)| r = rmin , 0 ≤ θ ≤ 2π } ,

(4.7)

where rmin is given in (2.4). We have the following lemma. Lemma 4.1. Let v ∈ H p ( S ) ( p > 32 ), then the bound exists,

I 0 ( v ) ≥ c 0 min v r 20,0 ,

(4.8)

v ∈V N

where v r = ∂∂vr , and c 0 (> 0) is a constant independent of N. Proof. For the harmonic functions v ∈ V N in S, we have from Oden and Reddy [41]

v 3 , S ≤ C v ν 0, ,

(4.9)

2

where v ν = ∂∂ νv . Then it follows from (4.2) that

I 0 ( v ) ≥ c 0 v 23 .

(4.10)

2 ,S

Moreover, since  v = 0 and S 0 ⊆ S, we have from [41]

v 3 , S ≥ v 3 , S 0 ≥ c 0 v r 0,0 . 2

(4.11)

2

Hence, we obtain from (4.10) and (4.11)

I 0 ( v ) ≥ c 0 v r 20,0 . This is the desired result (4.8), and completes the proof of Lemma 4.1.

(4.12)

2 



It is shown in [35] that even by the rough central rule, the integrals  v r2 and  v r2 are equivalent to each other, 0 0

v r 0,0  v r 0,0 , v ∈ V N ,

(4.13)

136

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145



where v 0,0 = { v 2 } 2 . The equation (4.13) means that there exist two positive constants C 1 and C 2 such that 1

0

 v r2 ≤

C1



 0

0

v r2 ≤ C 2

v r2 .

(4.14)

0

We have the following theorem. Theorem 4.1. Let (4.13) hold. Then the bounds exist,

rmin !2N

λmin (B) ≥ c 0 N

R

(4.15)

,

rmin ! N

λmin-next (B) ≥ c 0 N

R

(4.16)

.

Proof. We have

 "

 v r 20,0

v r2

=

=

0

=

#2 N ∂  c i φi ( P ) ∂r

N 

 ci c j

i , j =1

(4.17)

i =1

0

∂φi ( P ) ∂φ j ( P ) , P ∈ 0 . ∂r ∂r

0

Denote

P ∗ = (rmin cos h , rmin sin h) , where h =

2π N

 = 1, 2, · · · , N ,

. We have from (4.17) and (4.13) 2

v r 20,0 ≥ c 0 v r 0,0 N 

= c 0 hrmin

N  ∂φi ( P ∗ ) ∂φ j ( P ∗ )



ci c j

i , j =1

= c 0 hrmin

N 

(4.18)

=1

N 

ci c j

=1 i , j =1

∂r



∂r

∂φi ( P ∗ ) ∂φ j ( P ∗ ) = c 0 hrmin xT AT Ax, ∂r ∂r

where A is the symmetric and circulant matrix (3.9) with

x T A T Ax xT x

≥ min(λk2 (A)) ≥ λ20 (A) = c 0 N 2 k

rmin !2N R

ρ = rmin . From (3.19), the bound exists: (4.19)

.

Since v ∈ V N , v ∈ H p ( S )( p > 32 ), Lemma 4.1 holds. Hence, from Lemma 4.1, (4.4), (4.18) and (4.19), we obtain

λmin (B) ≥ c 0 ≥ c¯ 0 N

v r 20,0

rmin

x 2

≥ c 0 hrmin

!2N

x T A T Ax

x 2 (4.20)

,

R

where c¯ 0 (> 0) is also a constant independent of N. This is the desired result (4.15), and the proof for (4.16) is completed by noting

λmin-next (B) = min(λk2 (A)) ≥ λ2N (A) = c 0 N 2 k ≥1

rmin ! N

2

from (3.32). This completes the proof of Theorem 4.1.

R

,

2

From Theorem 4.1 and (4.6), we obtain the following corollary.

(4.21)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

137

Corollary 4.1. Let (4.13) hold. Then the bounds exist:

 Cond(B) ≤ C



 Cond(B) ≤ C

2N

R rmin R

(4.22)

, N

(4.23)

.

rmin

4.2. The collocation Trefftz methods In computation, we use (2.16) with M ≥ N, and obtain the over-determined system:

Fx = d,

(4.24)

where x = (c 1 , · · · , c N ) , d ∈ R , and F ∈ R T

1 2

M

M ×N

. We have

x T F T Fx = ˆI 0 ( v ),

where

ˆI 0 ( v ) =



( 

(4.25)

∂v 2 ) . ∂ν

(4.26)







Denote the discrete norms | v ν |0, = ( v 2ν ) 2 , where  is the numerical approximation of  by the Gaussian rule (or the trapezoidal rule). From [35], we assume that the equivalence exists: 1

v ν 0,  v ν 0, , v ∈ V N .

(4.27)

σi of F in (4.24), and

Denote the singular values

σmax = max σi , σmin-next > σmin = min σi > 0. i

i

We have the following theorem. Theorem 4.2. Let (4.13) and (4.27) hold. For the MFS in (2.16) with M ≥ N, the bounds exist:

√

σmax (F) = σ1 (F) ≤ C N , √

σmin (F) = σN (F) ≥ c0 N √

rmin ! N R

rmin ! 2

(4.28) (4.29)

,

N

σmin-next (F) ≥ c0 N

R

(4.30)

.

Proof. From (4.27), we have

ˆI 0 ( v )  I 0 ( v ), v ∈ V N ,

(4.31)

where ˆI 0 ( v ) and I 0 ( v ) are given in (4.26) and (4.2), respectively. From (4.31) and (4.6), we have

x T F T Fx

λmax (FT F) = max

xT x

x=0

≤ C max x=0 v ∈V N

I0(v )

x 2

= max x=0 v ∈V N

2 ˆI 0 ( v )

(4.32)

x 2

= C λmax (B) ≤ C¯ N ,

where C¯ is a constant independent of N. Next, from (4.31) and (4.15) we have

λmin (FT F) = min

x T F T Fx

x=0

≥ c 0 min

x=0 v ∈V N

I0(v )

x 2

xT x

= min

x=0 v ∈V N

2 ˆI 0 ( v )

(4.33)

x 2

= c 0 λmin (B) ≥ c¯0 N

rmin !2N R

,

138

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

where c¯ 0 (> 0) is also a constant independent of N. The desired results (4.28) and (4.29) follow from (4.32), (4.33) and the following equalities



1

σmax (F) = λmax (FT F)

2



1

σmin (F) = λmin (FT F)

,

2

.

Similarly, Eq. (4.30) follows from (4.16). This completes the proof of Theorem 4.2.

2

We have the following corollary from Theorem 4.2. Corollary 4.2. Let the conditions in Theorem 4.2 hold. Then the bounds exist:



Cond(F) =

 Cond(F) =

σmax (F) R ≤C σmin (F) rmin σmax (F) σmin-next (F)

N



≤C

(4.34)

, R

rmin

N 2

(4.35)

.

Corollary 4.2 indicates that the TSVD may greatly improve the stability of Neumann problems by the MFS. 5. Error estimates Since the MFS can be regarded as the TM choosing the FS, we may follow the analysis of TM in [24,35], and pay attention to the errors between the harmonic polynomials and the solutions of the MFS in [2,26,30]. To avoid multiple solutions, we may follow [15] by adding a constraint:

 

u  = d0 ,

(5.1)

D

where  the point is given by D ∈  and d0 is a constant (say, d0 = 1). Then the solution of (2.1) and (2.2) is unique due to  g = 0 and (5.1). Note that for v ∈ H 1 ( S ) satisfying (5.1), the norms v 1. S and | v |1, S are equivalent to each other: v 1, S  | v |1, S . In fact, Eq. (5.1) may not be used in our computation; see Section 6. Choose

uN =

N 

c i φi (r , θ), (r , θ) ∈ S ,

(5.2)

i =1

where

φi (r , θ) = ln

$

R 2 + r 2 − 2Rr cos(θ − ih), h =

2π N

.

(5.3)

Let V¯ N denote the set of the functions in (5.2) satisfying (5.1). Also denote



I (v ) =

( 

 ∂v ∂v − g )2 , ˆI ( v ) = ( − g )2 . ∂ν ∂ ν 

(5.4)

Then for the Neumann problem, the solutions of the TM and the CTM using FS are given by, respectively,

I (u N ) = min I ( v ),

(5.5)

ˆI (u˜ N ) = min  I ( v ).

(5.6)

v ∈ V¯ N

and v ∈ V¯ N

Below, we only derive the error bounds of u N by the TM (5.5), since the error bounds for (5.6) can be obtained similarly (see [35]). Denote the inter-product by,



[u , v ] B ∗ =

∂u ∂ v , ∂ν ∂ν

(5.7)



and the boundary norm by 1

v B ∗ = {[ v , v ] B ∗ } 2 =

∂v 0, . ∂ν

(5.8)

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

139

By following [24], the solution of u N also satisfies the minimum property:

u − u N B ∗ = min u − v B ∗ ,

(5.9)

v ∈ V¯ N

and the orthogonality property:

[u − u N , v ] B ∗ = 0, ∀ v ∈ V¯ N .

(5.10)

In Fig. 1, denote two disks S min and S max by

S min ⊂ S ⊂ S max ,

(5.11)

where

S min = {(r , θ)|r ≤ rmin , 0 ≤ θ ≤ 2π } , S max = {(r , θ)|r ≥ rmax , 0 ≤ θ ≤ 2π } ,

(5.12)

and rmax and rmin are given in (2.4). Denote the harmonic polynomial functions of order n by

P n = P n (r , θ) = where

α0 2

+

n 

r i (αi cos i θ + βi sin i θ),

(5.13)

i =1

αi and βi are coefficients, and denote the fundamental solutions by vN =

N 

c i ln | P Q i |,

(5.14)

i =1

where P ∈ ( S ∪ ∂ S ), and Q i are given in (2.8). We cite a lemma from Li et al. [30, p.140]. Lemma 5.1. Let (5.11) hold. When q ≥ 0 and

22q+1



R

−2N

rmax

≤ 1,

(5.15)

there exists a special FS, denoted by u ∗N ∈ V N , such that for N > 2n,



 2n− N  n ∂ R rmax ( P n (r , θ) − u ∗N ) q, ≤ C N q+1 P n (r , θ) 0, , ∂ν rmax rmin

(5.16)

where ν is the exterior normal to  , v q, is the Sobolev norm, and P n (r , θ) is given in (5.13). Lemma 5.2. Let  be a smooth boundary, condition (2.22) hold, and N satisfy



R rmax

2n− N 

rmax

n ≤

rmin

1 n

p − 12

.

(5.17)

Then for the Neumann problem, the solution u N by the MFS has the error bound,

u − u N B ∗ = O (

1 N

p − 32

) u p , S .

(5.18)

Proof. Let one true solution be

u = u (r , θ) = P n (r , θ) + R n (r , θ),

(5.19)

where P n (r , θ) is given in (5.13), and

R n = R n (r , θ) =

∞ 

r i (αi cos i θ + βi sin i θ),

(5.20)

i =n+1

and

αi and βi are the true coefficients. We have from (5.9) u − u N B ∗ ≤ u − u ∗N B ∗ ≤

∂ ∂ ( P n − u ∗N ) 0, + R n 0, . ∂ν ∂ν

(5.21)

140

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

When u ∈ H p ( S ) ( p ≥ 2)3 from (2.22), the bound exists [1],



∂ 1 R n 0, ≤ C R n 3 , S = O ( u p , S ). 2 p − 32 ∂ν n

(5.22)

When N is chosen to satisfy (5.17) and then N  n, we have from (5.16)



∂ 1 ( P N (r , θ) − u ∗N ) 0, = O ( P n (r , θ) 0, ). p − 32 ∂ν n

(5.23)

Since P n (r , θ) ≈ u, we have

P n (r , θ) 0, ≤ C u 0, ≤ C u 1 , S ≤ C u p , S .

(5.24)

2

Hence, Eq. (5.23) leads to



∂ 1 ( P n (r , θ) − u ∗N ) 0, = O ( ) u p , S . p − 32 ∂ν n

(5.25)

Combining (5.21), (5.22) and (5.25) and noting n  N yield the desired result (5.18), and this completes the proof of Lemma 5.2. 2 We need to estimate the bounds of negative norms. Define the negative norms,

 ∂w ·v ∂w −, = sup  ∂ ν . ∂ν v ∈ H  () v ,

(5.26)

Lemma 5.3. Let the conditions in Lemma 5.2 hold. Then the error bounds exist,

∂ 1 (u − u N ) − 1 , = C p −1 u p , S , 2 ∂ν N ∂ 1 (u − u N ) − 3 , = C p u p , S . 2 ∂ν N

3

Proof. Let t ∈ H 2 () and

(5.27) (5.28)

η ∈ H 3 ( S ) satisfy the following auxiliary problem,

η = 0 in S , ∂η = t on , ∂ν η = 1, D ∈ .

(5.29) (5.30) (5.31)

D

Then we have from (5.26) and (5.30)

∂(u − u N ) − 3 , ≤ 2 ∂ν

 

∂(u −u N ) ∂ν

·t

t 3 ,

 

=

∂(u −u N ) ∂ν

· ∂∂ ην

t 3 ,

2

.

(5.32)

2

Hence from (5.32), (5.7), (5.10) and Lemma 5.2



∂(u − u N ) ∂ η · = [u − u N , η ] B ∗ ∂ν ∂ν

(5.33)



= [u − u N , η − η N ] B ∗ ≤ u − u N B ∗ η − η N B ∗ ≤ C(

1 N

p − 32

u p , S ) × (

1 N

3 2

η 3, S ) =

C Np

u p , S × η 3, S .

However, from the auxiliary problem (5.29)–(5.31), we have

η 3, S ≤ C t 3 , .

(5.34)

2

The second desired result (5.28) follows from (5.32)–(5.34) and (5.26).

3

This assumption may be relaxed to u ∈ H p ( S ) ( p >

3 ) 2

for singularity problems with corners, and all error estimations of MFS in this paper hold.

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

141

Fig. 2. The Neumann problem of Laplace’s problem with the source circle  R .

Next, we prove the first result (5.27). From the Sobolev interpolation theory [1], the bound exists, 2

1

v − 1 , ≤ C { v 0, } 3 × { v − 3 , } 3 . 2

(5.35)

2

Hence, we have from (5.28), (5.35) and Lemma 5.2



2 1 ∂ ∂ ∂ (u − u N ) − 1 , ≤ C { (u − u N ) 0, } 3 × { (u − u N ) − 3 , } 3 2 2 ∂ν ∂ν ∂ν

≤ C{

1 N

p − 32

1

2

}3 × {

1

Np

} 3 u p , S = C

1 N p −1

u p , S .

This is the first desired result (5.27), and completes the proof of Lemma 5.3.

(5.36)

2

Theorem 5.1. Let the conditions in Lemma 5.2 hold. Then for the Neumann problem, the numerical solutions by the MFS (i.e., the TM using FS) have the optimal error bounds in both L 2 and H 1 norms in S,

u − u N k, S = O (

1 N p −k

) u p , S , k = 0, 1.

(5.37)

Proof. Since for the harmonic function (u − u N ), the bound exists:

u − u N k, S ≤ C

∂ (u − u N ) k− 3 , , 2 ∂ν

the desired results (5.37) follow from Lemma 5.3, and this completes the proof of Theorem 5.1.

(5.38)

2

For the quotient space H k ( S )/ P 0 ( S ), the differences of a constant solution are not included into (5.38). For the pure Neumann problem, the optimal rates in both L 2 and H 1 norms in S are derived in this section. Under (2.22), the polynomial convergence rates are obtained. When the Laplace’s solutions are infinitely smooth, the exponential convergence rates can also be achieved, as those by the MPS in [35]. Note that the proof here does not need the assumption of inverse inequalities in [24,35]. However, for the mixed boundary problems, the analytical approaches in [24,27,28] should be solicited, in terms of the assumption of inverse inequalities. This is a distinct feature in error analysis for Neumann problems. 6. Numerical experiments Consider the Neumann problem of Laplace’s equation (see Fig. 2):

u = 0, in S ,

∂u = g on ∂ S , ∂ν

(6.1)

142

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

Table 1 The errors and condition numbers by the MFS for R = 2 and M = 90. N

 0, S ,∞

 0, S

| |1, S

Cond

x

10 20 30 40

9.49(1) 9.32(−1) 2.70(−3) 7.03(−6)

6.31(1) 2.86(−1) 7.19(−4) 7.02(−7)

1.48(2) 1.99 6.41(−3) 8.97(−6)

3.88(2) 3.40(5) 1.19(8) 1.45(11)

3.13(3) 5.57(3) 3.46(3) 2.94(3)

Table 2 The errors and condition numbers by the MFS for R = 2 and M = 120. N

| 0, S ,∞

| 0, S

 1, S

Cond

x

10 20 30 40

9.49(1) 9.32(−1) 2.61(−3) 3.75(−6)

6.31(1) 2.86(−1) 6.68(−4) 3.68(−7)

1.48(2) 1.99 6.31(−3) 7.97(−6)

3.87(2) 3.40(5) 1.18(8) 8.84(10)

3.12(3) 5.57(3) 3.45(3) 2.94(3)

Table 3 The errors and condition numbers by the MFS for N = 20 and M = 90. R

√

3

2 3 4 6 8

 0, S ,∞

 0, S

| |1, S

Cond

x

1.85 9.32(−1) 3.23(−1) 2.43(−1) 2.07(−1) 1.98(−1)

4.37(−1) 2.86(−1) 1.25(−1) 9.47(−2) 7.93(−2) 7.50(−2)

3.58 1.99 8.56(−1) 7.09(−1) 6.52(−1) 6.41(−1)

2.18(4) 3.40(5) 7.65(8) 1.81(11) 4.02(14) 9.51(16)

2.16(3) 5.57(3) 1.16(5) 1.25(6) 4.06(7) 5.08(8)

 

where S = {(x, y ) − 1 < x < 1, 0 < y < 1}. Choose the smooth solution

u (x, y ) = cosh π (x + 1) cos(π y ).

(6.2)

Then the function g in (6.1) is given explicitly by

g = π sinh(2π ) cos(π y ), on A B ,

(6.3)

g = 0, on BC ∪ C D ∪ D A . The consistency condition is satisfied, since

 ∂S

1

 gd =

g = π sinh(2π )

cos π ydy = 0.

(6.4)

0

AB

Hence, the solutions of (6.1) exist, and the solution by adding an arbitrary constant is also a solution of (6.1). In order to provide the unique solution, we may add the constraint as in (5.1) (see [15]):

u (0, 0) = cosh π .

(6.5)

In Fig. 2, there are four corners A, B, C and D. Note that the smooth solution (6.2) is chosen, to reach the exponential convergence. Let M denote the number of uniform collocation nodes along ∂ S. We use the central rule, and always choose M > N so that Eq. (2.16) is an over-determined system. For R = 2, the errors in S and condition numbers are given in Tables 1 and 2 for M = 90 and M = 120, respectively. In the tables,  0, S , | |1, S and  0, S ,∞ (= max S |u |) are the Sobolev norms, where  = u − u N . From Tables 1 and 2 we can see the empirical rates:

 0, S ,∞ = O (0.58 N ),  0, S = O (0.54 N ), | |1, S = O (0.57 N ), N

Cond = O (2.0 ).

(6.6) (6.7)

From our computational experiments, when M is large enough, say M ≥ M 1 > N, different choices of M do not have an evident impact on errors and condition numbers. Table 3 lists the errors and the condition numbers for different R. From Table 3 (and Table 6), when R (> rmax ) increases, the condition number grows exponentially, to coincide with (4.34). Next, we may remove the condition (6.5). Once the coefficients c i are obtained from the MFS without condition (6.5), the general solutions are given in (2.10), where c is an arbitrary constant. In order to compare with the true solution (6.2),

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

143

Table 4 The errors and condition numbers by the MFS without (6.5) for R = 2 and M = 90. N

 0, S ,∞

 0, S

| |1, S

Cond

x

10 20 30 40

5.39(1) 7.62(−1) 2.26(−3) 6.86(−6)

2.46(1) 1.55(−1) 3.44(−4) 6.57(−7)

1.36(2) 1.97 6.38(−3) 8.97(−6)

3.88(2) 3.40(5) 1.19(8) 1.45(11)

3.13(3) 5.57(3) 3.46(3) 2.94(3)

Table 5 The errors and condition numbers by the MFS using TSVD for R = 2 and M = 90. N

 0, S ,∞

 0, S

| |1, S

ˆ Cond

x

10 20 30 40

5.39(1) 7.62(−1) 2.26(−3) 6.86(−6)

2.46(1) 1.55(−1) 3.44(−4) 6.57(−7)

1.36(2) 1.97 6.38(−3) 8.97(−6)

3.12(1) 2.62(3) 2.69(5) 3.00(7)

3.13(3) 5.57(3) 3.46(3) 2.94(3)

Table 6 The errors and condition numbers by the MFS using TSVD for N = 20 and M = 90. R

√

3

2 3 4 6 8

 0, S ,∞

 0, S

| |1, S

ˆ Cond

x

1.61 7.62(−1) 2.50(−1) 1.91(−1) 1.68(−1) 1.62(−1)

2.75(−1) 1.55(−1) 7.12(−2) 6.02(−2) 5.60(−2) 5.52(−2)

3.56 1.97 8.50(−1) 7.05(−1) 6.50(−1) 6.38(−1)

9.52(2) 2.62(3) 7.03(4) 8.75(5) 3.28(7) 4.35(8)

2.16(3) 5.57(3) 1.16(5) 1.25(6) 4.06(7) 5.08(8)

we obtain a constant c from (6.5) by

uN ( P ∗) = c +

N 

c i φi ( P ∗ ) = u ( P ∗ ) = cosh π , P ∗ = (0, 0).

(6.8)

c i φi ( P ∗ ), P ∗ = (0, 0).

(6.9)

i =1

Then we obtain

c = cosh π −

N  i =1

The errors can be easily evaluated and listed in Table 4. Compared Table 4 with Table 1, the errors are slightly smaller, but the condition numbers remain the same. Hence, we may not need the condition (6.5) in computation. To reduce the condition numbers, we use the TSVD in (3.34) to remove the component of u0 corresponding to λ0 . The results are listed in Tables 5 and 6. From Table 5, we can see that Eqs. (6.6) remain the same, and the modified condition number in (3.35) has the numerical rate,

 Cond = O (1.60 N ),

(6.10)

%

which is much smaller than that in (6.7). On the other hand, by noting the ratio

R rmin

=

%

2 0.5

= 2, the theoretical rate is

given from (4.35) by

   Cond

T heor y

= O (2 N ).

(6.11)

Both (6.7) and (6.10) are consistent with (6.11). Also from Tables 3 and 6, we see for different R,

√

Cond = O (( Cond)2 ),  Cond = O ( Cond),

(6.12)

to coincide with (4.34) and (4.35). The equations (6.12) indicate that by using the TSVD in (3.34), the stability can be improved significantly. From Tables 3 and 6, we can see that when R increases, the condition numbers Cond grow extremely quickly, but the errors decrease significantly. Such numerical results coincide with Corollary 4.2 and Lemma 5.1. For choosing R of the MFS, the balance between accuracy and stability must be maintained. A trial computation is necessary, as shown in Tables 1, 2, 3 and 6, so that good choices of R (i.e., R = 2) matching with M and N can be found, and used in Tables 4 and 5.

144

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

7. Concluding remarks To close this paper, let us address a few remarks. 1. For disk domains, the eigenvalues of the stiffness matrix from the MFS are expressed in Theorem 3.1 in terms of N ν and ak , which may easily lead to the sharp bounds of Cond in Corollary 3.1. The exponential growth rates of condition numbers are obtained. Note that the discrete matrix from the MFS is not singular for Neumann problems. For bounded simply-connected domains, the new and sharp bounds of Cond by the MFS are derived in Section 4, and the same exponential rates are obtained. 2. To reduce condition numbers significantly, as shown in (6.12), the TSVD in (3.34) is used by removing the component of y0 corresponding to λ0 . This is a distinct feature of Neumann problems by the MFS. In contrast, for Dirichlet problems and the mixed boundary problems, the improvement of stability by the TSVD is not so significant, see [29]. 3. Let us compare the results in this paper with other boundary conditions of Laplace’s equation in [26–28]. The optimal convergence in L 2 norms in S is derived in [26] for Dirichlet problems in both simply-connected domains and annular domains. In this paper, the error bounds in L 2 and H 1 norms in S are derived for Neumann problems without using the inverse inequalities in [24,27,28], to give the optimal polynomial convergence rates. For the Dirichlet problems in disk domains with radius ρ , there exist the sharp bounds in [28],

Cond ≤

N | ln R |

 N R

ρ

2

2

, ρ < R = 1,

(7.1)

 N

Cond ≤ C N

R

, ρ < R = 1.

ρ

(7.2)

For the mixed boundary value problems in the simply-connected domains, the following bounds are derived in [28],

 Cond(F) ≤ C N

2

R rmin

N 2

.

(7.3)

The bounds of (7.1) and (7.3) are slightly larger than those of (3.37) and (4.35) of  Cond(F), respectively. Note that the MFS for the Neumann problem in this paper is more interesting and intriguing due to its distinct features, and the new analysis and the numerical techniques may enhance the application of the MFS. 4. Recently, our efforts have been paid to establish the theoretical analysis of the MFS. This paper is devoted to the analysis for Neumann problems of Laplace’s equation. This paper and Li et al. [26–28] have formulated the basic approaches of accuracy and stability of the MFS. Based on this theoretical framework, the analysis of the MFS for other problems has been developed, and the results are reported for biharmonic equations in [33], for linear elastostatics in [30], for unbounded domain problems in [36], and for the Cauchy problems in [48]. Combined methods of the MFS with the MPS for singularity problems are reported in [25,32,36,39], and comparisons of the MFS with the MPS in [37,48]. For the MFS, the optimal errors can always be achieved, as those in the MPS in [35]. This explains the fact the MFS is very effective for wide engineering problems, to support numerical computation in many papers, such as [3,5,11,12,44,45]. However, since the condition numbers grow exponentially, to seek techniques for better stability is significant in real applications. Some progress has been made in this paper. Acknowledgements The authors are grateful to the anonymous reviewers for their valuable suggestions, and indebted to L.F. Lo for the computation in Section 6. References [1] I. Babuska, A.K. Aziz, Survey lectures on the mathematical foundations of the finite element method, in: A.K. Aziz (Ed.), The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, Academic Press, New York, 1972, pp. 3–359. [2] A. Bogomolny, Fundamental solutions method for elliptic boundary value problems, SIAM J. Numer. Anal. 22 (1985) 644–669. [3] C.S. Chen, Y.C. Hon, R.A. Schaback (Eds.), Scientific Computing with Radial Basis Functions, Department of Mathematics, University of Southern Mississippi, Hattiesburg, MS 39406, USA, 2005, preprint. [4] C.S. Chen, A. Karageorghis, Y. Li, On choosing the location of the sources in the MFS, Numer. Algorithms 72 (2016) 107–130. [5] C.S. Chen, A. Karageorghis, Y.S. Smyrlis (Eds.), The Method of Fundamental Solutions – A Meshless Method, The Proceedings of the First International Workshop on the Method of Fundamental Solutions (MFS2007) held in Ayia Napa, Cyprus, June 11–13, 2007, Dynamic Publishers, Inc., USA, 2008. [6] S. Christiansen, On Kupradze’s functional equations for plane harmonic problems, in: R.P. Gilbert, R. Weinacht (Eds.), Function Theoretic Method in Differential Equations, in: Research Notes in Mathematics, vol. 8, London, 1976, pp. 205–243. [7] S. Christiansen, Condition number of matrices derived from two classes of integral equations, Math. Methods Appl. Sci. 3 (1981) 364–392. [8] S. Christiansen, On the two method for elimination of non-unique solutions of an integral equation of an integral equation with logarithmic kernel, Appl. Anal. 13 (1982) 1–18. [9] P.J. Davis, Circulant Matrices, John Wiley & Sons, New York, 1979. [10] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second ed., Academic Press, 1984, p. 315.

Z.-C. Li et al. / Applied Numerical Mathematics 119 (2017) 126–145

145

[11] G. Fairweather, A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Adv. Comput. Math. 9 (1998) 69–95. [12] M.A. Golberg, C.S. Chen, The method of fundamental solutions for potential, Helmholtz and diffusion problems, in: M.A. Golberg (Ed.), Boundary Integral Methods: Numerical and Mathematical Aspects, WIT Press, 1998, pp. 103–176. [13] I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series and Products, Academic Press, New York, 1965. [14] P.C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, Philadelphia, 1997. [15] H.T. Huang, Z.C. Li, N. Yan, New error estimates of Adini’s elements for Poisson’s equation, Appl. Numer. Math. 50 (2004) 49–74. [16] M. Katsurada, A mathematical study of the charge simulation method II, J. Fac. Sci., Univ. Tokyo, Sect. 1A, Math. 36 (1989) 135–162. [17] M. Katsurada, Asymptotic error analysis of the charge simulation method in a Jordan region with an analytical boundary, J. Fac. Sci., Univ. Tokyo, Sect. 1A, Math. 37 (1990) 635–657. [18] M. Katsurada, Charge simulation method using exterior mapping functions, Jpn. J. Ind. Appl. Math. 11 (1994) 47–61. [19] M. Katsurada, H. Okamoto, The collocation points of the method of fundamental solutions for the potential problem, Comput. Math. Appl. 31 (1996) 123–137. [20] T. Kitagawa, Asymptotic stability of the fundamental solution methods, J. Comput. Appl. Math. 38 (1991) 263–269. [21] V.D. Kupradze, Potential methods in elasticity, in: J.N. Sneddon, R. Hill (Eds.), Progress in Solid Mechanics, vol. III, Amsterdam, 1963, pp. 1–259. [22] X. Li, Convergence of the method of fundamental solutions for solving the boundary value problem of modified Helmholtz equation, Appl. Comput. Math. 159 (2004) 113–125. [23] X. Li, On convergence of the method of fundamental solutions for solving the Dirichlet problem of Poisson’s equation, Adv. Comput. Math. 23 (2005) 227–265. [24] Z.C. Li, Combined Methods for Elliptic Equations with Singularities, Interface and Infinities, Academic Publishers, Dordrecht, Boston and London, 1998. [25] Z.C. Li, Combinations of method of fundamental solutions for Laplace’s equation with singularities, Eng. Anal. Bound. Elem. 32 (2008) 859–869. [26] Z.C. Li, The method of fundamental solutions for annular shaped domains, J. Comput. Appl. Math. 228 (2009) 355–372. [27] Z.C. Li, The method of fundamental solutions for Laplace’s equation with mixed boundary problems, in: C.S. Chen, A. Karageorghis, Y.S. Smyrlis (Eds.), The Method of Fundamental Solutions – A Meshless Method, Dynamic Publishers, Inc., 2008, pp. 29–49. [28] Z.C. Li, J. Huang, H.T. Huang, Stability analysis of method of fundamental solutions for mixed boundary value problems of Laplace’s equation, Computing 88 (2010) 1–29. [29] Z.C. Li, H.T. Huang, Y. Wei, Ill-conditioning of the truncation singular value decomposition and the Tikhonov regularization and their application to numerical partial differential equations, Numer. Linear Algebra Appl. 18 (2011) 205–221. [30] Z.C. Li, H.T. Huang, M.G. Lee, J.Y. Chiang, Error analysis of the method of fundamental solutions for linear elastostatics, J. Comput. Appl. Math. 251 (2013) 133–153. [31] Z.C. Li, H.T. Huang, Y. Wei, A.H.-D. Cheng, Effective Condition Number for Partial Differential Equations, Science Press, Beijing, China, 2013. [32] Z.C. Li, M.G. Lee, J.Y. Chiang, The collocation Trefftz methods for the Stokes equations with singularity, Numer. Methods Partial Differ. Equ. 29 (2013) 361–395. [33] Z.C. Li, M.G. Lee, J.Y. Chiang, Y.P. Liu, The Trefftz method using fundamental solutions for biharmonic equations, J. Comput. Appl. Math. 235 (2011) 4350–4367. [34] Z.C. Li, T.T. Lu, H.Y. Hu, A.H.-D. Cheng, Particular solutions of Laplace’s equations on polygons and new models involving mild singularities, Eng. Anal. Bound. Elem. 29 (2005) 59–75. [35] Z.C. Li, T.T. Lu, H.Y. Hu, A.H.-D. Cheng, Trefftz and Collocation Methods, MIT, Southampton, Boston, 2008. [36] Z.C. Li, T.T. Lu, G.Q. Xie, The method of fundamental solutions and its combinations for exterior problems, Eng. Anal. Bound. Elem. 36 (2012) 394–403. [37] Z.C. Li, L.J. Young, H.T. Huang, Y.P. Liu, A.H.-D. Cheng, Comparisons of fundamental solutions and particular solutions for Trefftz methods, Eng. Anal. Bound. Elem. 34 (2010) 248–258. [38] Z.C. Li, T.T. Lu, H.T. Huang, A.H.-D. Cheng, Error analysis of Trefftz method for Laplace’s equations and its applications, Comput. Model. Eng. Sci. 52 (2009) 39–81. [39] M.G. Lee, L.J. Young, Z.C. Li, P.C. Chu, Combined Trefftz method of particular and fundamental solutions for corner and crack singularity of linear elastostatics, Eng. Anal. Bound. Elem. 34 (2010) 632–654. [40] R. Mathon, R.L. Johnston, The approximate solution of elliptic boundary-value problems by fundamental solutions, SIAM J. Numer. Anal. 14 (1977) 638–650. [41] J.T. Oden, J.N. Reddy, An Introduction to the Mathematical Theory of Finite Elements, John Wiley & Sons, New York, 1976. [42] P.A. Ramachandran, Method of fundamental solutions: singular value decomposition analysis, Commun. Numer. Methods Eng. 18 (2002) 789–801. [43] R. Schaback, Adaptive numerical solution of MFS systems, in: C.S. Chen, A. Karageorghis, Y.S. Smyrlis (Eds.), The Method of Fundamental Solutions – A Meshless Method, Dynamic Publishers, Inc., USA, 2009, pp. 1–28. [44] Y.S. Smyrlis, A. Karageorghis, Some aspect of the method of fundamental solutions for certain harmonic problems, J. Sci. Comput. 16 (2001) 341–371. [45] Y.S. Smyrlis, A. Karageorghis, Numerical analysis of the MFS for certain harmonic problems, Modél. Math. Anal. Numér. 38 (2004) 495–517. [46] E. Trefftz, Konvergenz und Ritz’schen Verfahren, in: Proc. 2nd Int. Cong. Appl. Mech., Zürich, 1926, pp. 131–137. [47] C.F. van Loan, Generalizing the singular value decomposition, SIAM J. Numer. Anal. 13 (1976) 76–83. [48] L.P. Zhang, Z.C. Li, Y. Wei, J.Y. Chiang, Cauchy problems of Laplace’s equation by methods of fundamental solutions and particular solutions, Eng. Anal. Bound. Elem. 37 (2013) 765–780.