Analysis on the method of fundamental solutions for biharmonic equations

Applied Mathematics and Computation 339 (2018) 346–366

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Analysis on the method of fundamental solutions for biharmonic equations Fangfang Dou a, Zi-Cai Li b, C.S. Chen c, Zhaolu Tian d,∗ a

School of Mathematical Sciences, University of Electronic Science and Technology of China, China Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan 80424, China c Department of Mathematics, University of Southern Mississippi, Hattiesburg, MS 39406, USA d College of Data Science, Taiyuan University of Technology, Taiyuan, China b

a r t i c l e

i n f o

MSC: 65N10 65N30 Keywords: Error analysis Stability analysis Biharmonic equations Method of fundamental solutions Trefftz methods

a b s t r a c t In this paper, the error and stability analysis of the method of fundamental solution (MFS) is explored for biharmonic equations. The bounds of errors are derived for the fundamental solutions r2 ln r in bounded simply-connected domains, and the polynomial convergence rates are obtained for certain smooth solutions. The bounds of condition number are also derived to show the exponential growth rates for disk domains. Numerical experiments are carried out to support the above analysis, which is the first time to provide the rigorous analysis of the MFS using r2 ln r for biharmonic equations. © 2018 Elsevier Inc. All rights reserved.

1. Introduction of the MFS The method of fundamental solutions (MFS) was first used in Kupradze [14] in 1963. Since then, there have appeared numerous reports of the MFS for computation such as the reviews of the MFS in Fairweather and Karageorghis [9] and Golberg and Chen [11]. On the other hand, the Trefftz method (TM) [27] as boundary methods has been fully developed in theory and computation for several decades (see [24]), where only the particular solutions (PS) are discussed. In fact, the MFS is one of the Trefft Methods using fundamental solutions (FS). Both the MFS and the MPS can be classified into TM, and they are carried out by the collocation TM (CTM) in [24]. In contrast to the traditional methods [3,4,10,16,28], we focus on the error and stability analysis of the MFS for the biharmonic equation in this paper. In the stability analysis, the study of condition number is also crucial [17,29]. Bounds of condition number will be further investigated in this paper. Consider the biharmonic equations and clamped boundary conditions

2 u = 0, u = f, uν =

in S, on  ,

∂u = g, ∂ν

(1.1) (1.2)

on  ,

(1.3)

where  = ∂ 2 /∂ x2 + ∂ 2 /∂ y2 , S is the bounded simply connected domain, uν the outward normal derivative to  , √  its smooth boundary, and f, g the given functions with sufficient smoothness. Denote r = |P Q |, P = ρ eiθ , Q = Reiφ , and i = −1. ∗

Corresponding author. E-mail address: [email protected] (Z. Tian).

https://doi.org/10.1016/j.amc.2018.07.016 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

347

The general traditional fundamental solutions (FS) of biharmonic equations in 2D are well-known as

(ρ , θ ) = r2 ln r, r =



R2 + ρ 2 − 2Rρ cos(θ − φ ),

(1.4)

where Q is the source point outside the domain. Denote the FS at Q j = Reiφ j as

 j (ρ , θ ) = r2j ln r j , φ j (ρ , θ ) = ln r j ,  where r j = |P Q j | = R2 + ρ 2 − 2Rρ cos(θ − φ j ). Hence we may also choose the linear combination

vN =

N 

(c j  j (ρ , θ ) + d j φ j (ρ , θ )),

(1.5)

(1.6)

j=1

where {cj } and {dj } are the unknown coefficients to be determined by the boundary conditions (1.2) and (1.3). Denote VN the set of (1.6). We may solicit the Trefftz method [24] to seek the numerical solution uN by

I (uN ) = min I (v ),

(1.7)

v∈VN

where

I (v ) =

 

( v − f )2 + w2

 

( vν − g )2 ,

(1.8)

and w is the weight which may be chosen as w = 1/N (see [24]). The TM using the FS is exactly the very MFS (see [20,22]). The Almansi’s FS for biharmonic equations are obtained directly from the general solutions, u(ρ , θ ) = ρ 2 v + z, where v and z are harmonic functions, to give A (ρ , θ ) = ρ 2 ln r. We may choose the linear combination

vAN =

N 

(c j Aj (ρ , θ ) + d j φ j (ρ , θ )),

(1.9)

j=1

to replace (1.6), where Aj (ρ , θ ) = ρ 2 ln r j , and the coefficients cj and dj are also obtained from (1.7). Eq. (1.7) with (1.9) is called the MFS using Almansi’s FS: ρ 2 ln r. For biharmonic equation, the MFS using ρ 2 ln r was introduced in Karageorghis and Fairweather [9,12,13]. The analysis of the MFS for Laplace’s equations can be found in [5,18,20]. In Li et al. [22], a preliminary analysis of error and stability was given for the MFS using ρ 2 ln r. Note that the analysis for the MFS using r2 ln r is more challenging and important. This paper is organized as follows. In Sections 2 and 3, the preliminary lemmas and the main theorems are derived for error bounds, respectively. In Section 4, the bounds of condition numbers are derived for circular domains. In the last section, numerical experiments are reported to illustrate our results. 2. Preliminary lemmas Firstly let us cite lemmas in [1,20]. Lemma 2.1. For the periodical function g = g(x ) on [0, 2π ], denote the errors

N (g) = TN (g) −

 2π 0

g(x )dx,

where the trapezoidal rule, TN (g) = h

(2.1)

N−1 k=0

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

N (sin mx ) = 0, m = 1, 2, . . . ,



N (cos mx ) =

2π , 0,

if m = ν N, otherwise.

(2.2)

ν = 1, 2, . . . ,

(2.3)

Lemma 2.2. Let g(x) ∈ C∞ [0, 2π ] satisfy the periodical boundary conditions,

g( ) (0 ) = g( ) (2π ),  = 0, 1, . . . There exist the errors

TN (g(x ) cos(kx )) − where

α =

1

π

 2π 0

 2π 0

g(x ) cos(kx )dx = π

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

(2.4) ∞ 

ν =1

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

(2.5)

(2.6)

348

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

From the general biharmonic solutions, u(ρ , θ ) = ρ 2 v + z, where v and z are harmonic functions, the particular solutions of biharmonic equations can be represented as



u (ρ , θ ) = ρ

2

∞

a∗0  i ∗ + ρ (ai cos iθ + b∗i sin iθ ) 2



∞

+

i=1

a0  i + ρ (ai cos iθ + bi sin iθ ), 2

(2.7)

i=1

where a∗i , b∗i , ai and bi are the true coefficients. Denote two harmonic polynomials of degree n,

PnH (ρ , θ ) =

a∗0  i ∗ + ρ (ai cos iθ + b∗i sin iθ ), 2

(2.8)

a0  i + ρ (ai cos iθ + bi sin iθ ), 2

(2.9)

n

i=1

Pnh (ρ , θ ) =

n

i=1

and their residuals

RH n (ρ , θ ) =

∞ 

ρ i (a∗i cos iθ + b∗i sin iθ ),

(2.10)

ρ i (ai cos iθ + bi sin iθ ).

(2.11)

i=n+1

Rhn (ρ , θ ) =

∞  i=n+1

Then we have

u(ρ , θ ) = PnB + RBn ,

(2.12)

PnB = PnB (ρ , θ ) = ρ 2 PnH (ρ , θ ) + Pnh (ρ , θ ),

(2.13)

where h RBn = RBn (ρ , θ ) = ρ 2 RH n ( ρ , θ ) + Rn ( ρ , θ ).

Suppose that the weak solution satisfies u ∈ Hσ (S) with σ ≥ 5/2. There exist the bounds of the residuals [6],

RBn (ρ , θ )k,S ≤ C

1 uσ ,S , nσ −k

RHn (ρ , θ )k,S , Rhn (ρ , θ )k,S ≤ C

1 uσ ,S , nσ −k

(2.14)

where uσ , S are the Sobolev norms. Lemma 2.3. For ρ < R, there exists the expansion of (1.4),

(ρ , θ ) = R2 ln R + ρ 2 (1 + ln R ) + R2 − ρ2

∞  m=1

m

1 ρ m (m + 1 ) R

βm ρ m

∞  m=1

m

R

cos m(θ − φ )

(2.15)

cos m(θ − φ ),

where β1 = −(1 + 2 ln R ) and βm = 1/(m − 1 ) (m ≥ 2 ). Proof. From [5,18] there exists the expansion of the FS of Laplace’s equation,

ln r = ln R −



∞  1 ρ n cos n(θ − φ ), n R

ρ < R.

(2.16)

n=1

For biharmonic equations, the traditional FS, r2 ln r, is expressed as



(ρ , θ ) = (R + ρ 2

2

n=1

From the trigonometric identity

2 cos θ cos nθ = cos(n + 1 )θ + cos(n − 1 )θ , we have

2 cos(θ − φ )



∞  1 ρ n cos n(θ − φ ) n R n=1





∞  1 ρ n − 2Rρ cos(θ − φ )) ln R − cos n(θ − φ ) . n R

(2.17)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

= =

349



∞  1 ρ n (cos(n + 1 )(θ − φ ) + cos(n − 1 )(θ − φ ) ) n R n=1 ∞ 

m=2

1 m−1

ρ m−1 R

∞ 

cos m(θ − φ ) +

m=0

m+1

ρ 1 m+1 R

cos m(θ − φ ),

where we have used transformation n ± 1 = m. Hence we obtain from (2.17)





∞  1 ρ n (ρ , θ ) = (R2 + ρ 2 − 2Rρ cos(θ − φ )) ln R − cos n(θ − φ ) n R n=1

= (R2 + ρ 2 − 2Rρ cos(θ − φ )) ln R − (R2 + ρ 2 )



∞  1 ρ n + 2Rρ cos(θ − φ ) cos n(θ − φ ) n R

(2.18)





∞  1 ρ n cos n(θ − φ ) n R n=1

n=1

= [R2 + ρ 2 − 2Rρ cos(θ − φ )] ln R

ρ m ρ

∞  1 1 + R2 − cos m(θ − φ ) − R2 cos(θ − φ ) m−1 m R R + ρ2

m=2 ∞ 

1 m+1

m=1

1 m

−

ρ m R

= R2 ln R + ρ 2 (1 + ln R ) + R2 − ρ2

m

∞ 

ρ 1 m (m + 1 ) R

m=1

∞  m=1

cos m(θ − φ ) + Rρ

βm ρ m m

R

ρ

R

cos m(θ − φ )

cos m(θ − φ ).

This is the desired result (2.15), and completes the proof of Lemma 2.3.



The expansion (2.15) is also given by Chen et al. [7]. Lemma 2.4. For R > ρ , and R = e−1 , there exist the integration formulas:

1 = 2π (1 + ln R )

ρ

2

ρ

m+2

cos mθ sin mθ



2π 0

 =−

r ln rdφ − R 2

Rm (m(m + 1 ))

+

π

βm R

2

 2π 0

2

ln rdφ ,

0





 2π

2π

0



cos mφ r ln r dφ sin mφ 2





cos mφ ln r dφ , sin mφ

where r = |P Q |, P = ρ eiθ , Q = Reiφ and i =

(2.19)

m ≥ 1,

(2.20)

√ −1.

Proof. By using the orthogonality of trigonometric functions, we have from (2.15)

 2π 0

r 2 ln rdφ =

 2π 0

(ρ , θ )dφ

= 2π (R2 ln R + ρ 2 (1 + ln R ))  2π = R2 ln rdφ + 2π ρ 2 (1 + ln R ),

(2.21)

0

 2π where we have used 2π ln R = 0 ln rdφ from (2.16). This gives the first result (2.19). Next, from (2.16) there exists equality for m ≥ 1

 2π 0



π cos mφ ln r dφ = − sin mφ m

ρ m R

From (2.15) and (2.22) for m ≥ 1, we have



cos mθ . sin mθ

(2.22)

350

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

 2π 0

  ρ m βm π R2 ρ m cos mθ 1 cos mθ 2 − πρ sin mθ sin mθ m R m (m + 1 ) R

   2π

1 ρ m cos mθ cos mφ = −βm R2 ln r dφ − π ρ 2 , sin mφ sin mθ m (m + 1 ) R 0



cos mφ r ln r dφ = sin mφ

then give

ρ

2

m+2

cos mθ sin mθ

 =−

Rm m ( m + 1 )



π

2π 0



cos mφ r ln r dφ + βm R2 sin mφ 2

This is the second result (2.20), and completes the proof of Lemma 2.4.

 2π 0



(2.23)



cos mφ ln r dφ . sin mφ

(2.24)



Based on Lemma 2.4, we seek (1.6) with special coefficients c¯i and d¯i , to better approximate the major term ρ 2 PnH of u in (2.12). Choose the source points Qj to be located uniformly on the circle {(ρ , θ )|ρ = R, 0 ≤ θ ≤ π }, denoted as Q j = {R cos jh, R sin jh} with h = 2π /N. By the trapezoidal rule we have the approximation from (2.19):



ρ

2

N N   h ≈ r 2j ln r j − R2 ln r j 2π (1 + ln R )



h = 2π (1 + lnR )

j=1

N 



j=1

 j ( ρ , θ ) − R2

j=1

N 



φ j (ρ , θ ) ,

(2.25)

j=1

where  j (ρ , θ ) = r 2j ln r j , φ j (ρ , θ ) = ln r j , r j = |P Q j | and P = {ρ cos θ , ρ sin θ }. Similarly, from (2.20), we have

ρ

m+2

cos mθ sin mθ



≈−

Rm m ( m + 1 )

π



h

N  j=1

cos m jh  j (ρ , θ ) sin m jh



+ βm R

2

N  j=1

cos m jh φ j (ρ , θ ) sin m jh



, m ≥ 1.

(2.26)

From (2.13), (2.25), and (2.26), we have

ρ 2 PnH (ρ , θ ) =

a∗0 2  m+2 ∗ ρ + ρ (am cos mθ + b∗m sin mθ ) 2 n

m=1

≈

N 

(ck∗ k (ρ , θ ) + dk∗ φk (ρ , θ )),

(2.27)

k=1

where the coefficients are given by

ck∗ =

n h  m h a∗0 − R m(m + 1 )(a∗m cos mkh + b∗m sin mkh ), 4π (1 + ln R ) π

(2.28)

m=1

dk∗ = −

n R2 h R2 h  a∗0 + βm Rm m(m + 1 )(a∗m cos mhk + b∗m sin mkh ), R = e−1 . 4π (1 + ln R ) π

(2.29)

m=1

Similarly from (2.9), we find a better approximation (see [18])

Pnh (ρ , θ ) =

 a0  m + ρ (am cos mθ + bm sin mθ ) ≈ dkh φk (r, θ ), 2 n

N

m=1

k=1

(2.30)

where

dkh =

n h  h a0 − mRm (am cos mkh + bm sin mkh ), R = 1. 4π ln R π

(2.31)

m=1

Hence from (2.27)–(2.31), when R = e−1 and R = 1, we have found a better approximation for PNB in (2.13)

PnB ≈ N (PnB ; ρ , θ ) :=

N 

ci∗ i (ρ , θ ) +

i=1

N 

d¯i φi (ρ , θ ),

(2.32)

i=1

where d¯i = dkh + dk∗ , and dkh and dk∗ are the coefficients in (2.31) and (2.29), respectively. Based on Lemmas 2.1 and 2.2, the errors of the trapezoidal rule exist only for the terms of ρ m+2 cos mθ and m ρ cos mθ (m ≥ 1). Since the errors of harmonic polynomials ρ m cos mθ have been derived in [18,22,23] already, we will focus on estimation on the bounds of the following errors

ρ 2 PnH (ρ , θ ) − N (ρ 2 PnH ; ρ , θ ) =

n  i=1

a∗i (ρ i+2 cos iθ − N (ρ i+2 cos iθ ; ρ , θ )).

(2.33)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

From (2.26), we have

N ( ρ

m+2 m



cos mθ ; ρ , θ ) = −

N Rm m ( m + 1 )h 

π

cos m jh  j (ρ , θ ) − R

j=1

2

N 

351

βm cos m jh φ j (ρ , θ ) , m ≥ 1.

(2.34)

j=1

Denote the errors

E (ρ m+2 cos mθ ) = ρ m+2 cos mθ − N (ρ m+2 cos mθ ; ρ , θ ),

(2.35)

where  N is also the discrete approximation from the trapezoidal rule. Below we estimate the bounds of (2.35) by the Fouries series. For the circle ρ = {(r, θ )| r = ρ , 0 ≤ θ ≤ 2π } we have

E (ρ m+2 cos mθ ) =

∞

a0,m  + (ak,m cos kθ + bk,m sin kθ ), 2

(2.36)

k=1

where

ak,m = bk,m =

1

π 1

 2π 0

 2π

π

0

E (ρ m+2 cos mθ ) cos kθ dθ ,

k = 0, 1, . . . ,

(2.37)

E (ρ m+2 cos mθ ) sin kθ dθ ,

k = 1, 2, . . . .

(2.38)

From Lemma 2.1,

bk,m = 0 for all k, m.

(2.39)

To estimate the bounds of ak, m , we have the following lemma. Lemma 2.5. For the Fourier expansion (2.36), the coefficients are given by



ak,m (m ≥ 1 ) =

 k Rm m(m + 1 ) ρR



R2 k

(βm + βk ) −

if k = ±m + ν N,

0,

where ν = 1, 2, . . . , β1 = −(1 + 2 ln R ), βm =

1 m−1



, if k = ±m + ν N,

k(k+1 )

(2.40)

(m ≥ 2 ), and

⎧

 ⎨ (ρ /R )k R2 ρ2 (1 + βk ) − , ak,0 = (1 + ln R ) k k (k + 1 ) ⎩ 0,

ρ2

if k = ν N,

ν = 1, 2, . . . ,

(2.41)

otherwise,

Proof. First, for m ≥ 1 we have from (2.26),

1

 2π

(ρ m+2 cos mθ ) cos kθ dθ  N  1 Rm m ( m + 1 ) h  2π − − cos jmh j (ρ , θ ) cos kθ dθ π π 0 j=1 

ak,m =

π



0

−R2

N  2π  j=1

0

βm cos jmhφ j (ρ , θ ) cos kθ dθ

(2.42)

:= T1 − T2 . In (2.42) the second term T2 is just the integration approximation of the first term

T1 =

 2π

1

π

0

(ρ m+2 cos mθ ) cos kθ dθ

by the trapezoidal rule. We may express T1 (m ≥ 1) as

T1 = =

1

π 1

π

Rm m ( m + 1 ) 2 Rm m ( m + 1 )

 2π 0



0

2π

r 2 ln r cos kθ dθ + βm R2

ρ k  R2 βk − R

k

ρ2

k (k + 1 )

 + βm

 2π

R2 k

0

 ×

 ln r cos kθ dθ  2π 0

where we have used (2.15) and (2.16). Hence, from Lemma 2.2 we obtain

ak,m = T1 − T2

cos mφ dφ

cos mφ cos kφ dφ ,

352

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

1

=

R m (m + 1 ) m

π

ρ k  R2 R

k

(βk + βm ) −



ρ2

k (k + 1 )

N (cos mφ cos kφ ),

(2.44)

where

N (g) = TN (g) −

 2π 0

g(x )dx

and TN is the trapezoidal rule. Then we have for k = ±m + ν N,

ak,m = R m(m + 1 ) m

ρ k R 2 R

k

(βm + βk ) −

k (k + 1 )

otherwise ak,m = 0. This is the result (2.40). Next, for m = 0 there exists the equality from (2.19)

1

ak,0 =

 2π

π

0

ρ

2

h cos kθ dθ − 2π (1 + ln R )



ρ2



 2π 0

N 

,

(2.45)

(ρ , θ ) − R

2

i=1

N 

φi (ρ , θ ) cos kθ dθ

i=1

:= T1 − T2 .

(2.46)

For T1 we have form (2.16) and (2.15), similarly

T1 =

1

π

 2π

ρ 2 cos kθ dθ

  2 π  2 π  2π 1 2 2 r ln r cos k θ d θ − R ln r cos θ d θ dφ 2π 2 (1 + ln R ) 0 0 0    2π

2 (ρ /R )k R ρ2 21 = β − +R cos kφ dφ . 2π (1 + ln R ) 0 k k k (k + 1 ) k 0

=

(2.47)

Hence we obtain

ak,0 =

2  (ρ /R )k R ρ2 (1 + βk ) − N cos kφ , 2π (1 + ln R ) k k (1 + k )

and from Lemma 2.1 for k ≥ 1,



 k ak,0 = (1+1ln R ) ρR ak,0 = 0,



R2 k

(1 + βk ) −

ρ2



k(k+1 )



a0,0

 2π 0

ρ

2

ν = 1, 2, . . . ,

(2.49)

otherwise.

For m = 0 and k = 0 we have

1 = 2π

if k = ν N,

,

(2.48)

h dθ − 2π (1 + lnR )

 2π 0



N 

i (ρ , θ ) − R

2

i=1

N 





φi ( ρ , θ ) d θ ,

i=1

to give

a0,0 =

ρ 2  ( 1 ) = 0.

(2.50)

Combining (2.49) and (2.50) gives the third desired results (2.41). This completes the proof of lemma 2.5.



Lemma 2.6. Let N > 2n and n ≥ m. The Fourier coefficients in (2.40) and (2.41) have the bounds,

|ak,m | ≤ CRm−k ρ k ,

(2.51)

where k = ±m + ν N, ν = 1, 2, . . . , and C is a constant independent of N. Proof. From Lemma 2.5, we have

|ak,m | ≤ CR

m−k

ρ

k



m m (m + 1 ) (m + 1 )(|βm | + |βk | ) + , m ≥ 1. k k (k + 1 )

(2.52)

Firstly let k = ν N − m for ν ≥ 1. Since k = ν N − m ≥ N − m ≥ m, it is easy to see that the constants in the brackets of (2.52) are bounded by noting β1 = −(1 + 2 ln R ) and βm = 1/(m − 1 ) (m ≥ 2 ). Detailed proof is given in [19]. Secondly, let k = ν N + m. Since ν N + m > ν N − m, Eq. (2.51) also holds. This completes the proof of Lemma 2.6. 

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

353

3. Error bounds We are ready to give a theorem. Theorem 3.1. Let 22q+1 (R/ρ )−2N ≤ 1 hold, there exists the bound for q ≥ 0,

ρ 2 PnH (ρ , θ ) − N (ρ 2 PnH (ρ , θ ); ρ , θ )2q,ρ ≤ CNq

R 2n−N r n max rmax

rmin

ρ 2 PnH , (r, θ )0, ,

(3.1)

where rmax = minS r, rmin is the maximal radial of Sρ ⊆S, Sρ is a disk, and C is a constant independent of N. Proof. From (2.36), (2.39), and Lemmas 2.5, 2.6, under 22q+1 (R/ρ )−2N ≤ 1 we have by the arguments in [5,18]

E (ρ m+2 cos mθ )2q,ρ ≤ C

∞  k=1

R 2(m−N )

≤ CN 2q

∞ 

a2k,m k2q ≤ C

ρ

R2(m−k ) ρ 2k k2q

k=±m+ν N

R2m .

(3.2)

For 0 < ρ0 ≤ ρ˜ ≤ ρ , from (2.33) and the Schwarz inequality, we have

ρ 2 PnH (ρ , θ ) − N (ρ 2 PnH (ρ , θ ); ρ , θ )2q,ρ =   ≤C

n 

(

) ρ˜

a∗m 2

2m+2

n 

1

ρ˜ 2m+2 m=1

m=1

n  m=1

a∗m E (ρ m+2 cos mθ )2q,ρ

E (ρ m+2 cos mθ )2q,ρ .

(3.3)

From (2.27), there exists the bound



ρ 2 PnH (ρ , θ )20,ρ˜ = ρ˜ 4

a ∗ 2 0

2

+

∞ 

ρ˜ 2m ((a∗m )2 + (b∗m )2 ) ≥ ρ˜ 4

m=1

n 

ρ˜ 2m (a∗m )2 ,

(3.4)

m=1

and from (3.2) n  m=1

1

ρ

˜ 2m+2

≤ CN 2q

E (ρ m+2 cos mθ )2q,ρ ≤ C

R 4n−2N ρ 2n ρ

ρ˜

≤ CN

n 

N 2q

R 4(m−N ) ρ 2m

ρ ρ˜ m=1 R 4n−2N r 2n max 2q rmax

.

rmin

(3.5)

The desired result (3.1) follows from (3.3)–(3.5), and this completes the proof of Theorem 3.1. Theorem 3.2. Let

PnB

=

ρ 2 PnH

+ Pnh



in (2.12). There exists the bound

PnB − N (PnB ; ρ , θ )q,ρ ≤ CNq Proof. From [5,18], we have



R 2n−N r n max rmax

n 

Pnh − N (Pnh ; ρ , θ )q,ρ ≤ C

rmin

12 (am )2 ρ˜ 2m

m=1

PnB 0, .

(3.6)

  n 1   m E ( ρ cos m θ )   , ρ˜ 2m q,ρ m=1

(3.7)

and from (3.3) and (3.7)

  PnB − N (PnB ; ρ , θ )q,ρ ≤ ρ 2 PnH − N (ρ 2 PnH ; ρ , θ )

q,ρ

≤ CN

q

≤ CN q This is the desired result (3.6).

+ Pnh − N (Pnh ; ρ , θ )q,ρ

 R 2n−N r 2n  n  max rmax

rmin

R 2n−N r 2n max rmax

rmin

(

) ρ˜

a∗m 2

2m+2

+

a2m

ρ˜

2m

12 

m=1

PnB 0, .

(3.8)



We need the norm of normal derivatives. Define the norm

 ∞  ∂ 2 2q+2 2m−2 2 2  f (θ )q,ρ = π ρ m ρ ( am + bm ) , q ≥ 0, ∂ν m=1

(3.9)

354

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

where

f (θ ) =

∞

a0  i + ρ (ai cos iθ + bi sin iθ ). 2 i=1

According to the arguments in Theorems 3.1 and 3.2, we have the following theorem. Theorem 3.3. There exists the bound for q ≥ 0,

   ∂ B   ∂ν Pn − N (PnB ; ρ , θ ) 

Denote



[u, v]H =

q,ρ

≤ CN q+1

R 2n−N r n max rmax

rmin

PnB 0, .

(3.10)

 

uv + w2

with the norm vH =



uν vν ,



(3.11)

[u, v]H . The MFS (1.7) for (1.1)–(1.3) reads

u − uN H = min u − vH .

(3.12)

v∈VN

Since there exists the orthogonality [u − uN , v]H = 0, Eq. (3.12) can be rewritten as



[uN , v]H =





f u + w2



g · vν ,

∀v ∈ VN .

(3.13)

We have a main theorem. Theorem 3.4. Suppose that u ∈ Hp (S) (p ≥ 5/2), and N is chosen such that

R 2n−N r n max rmax

rmin

1 . n p+1/2

≤

(3.14)

When w = 1/N there exists the error bound

u − uN H = O



1



N p−1/2

u p,S .

(3.15)

Proof. For the residuals of the solution in (2.12), there exist the bounds

RBn q,S ≤ C

1 u p,S , n p−q



∂ B 1 R  ≤ C p−q−1 u p,S . ∂ν n q,S n

(3.16)

Choose u¯ (∈ VN ) = N (PnB ; ρ , θ ). From (3.12), we have

u − uN H ≤ u − u¯ H ≤ PnB − N (PnB ; ρ , θ )H + RBn H .

(3.17)

From (3.16) and w = 1/N,

RBn H ≤ RBn 0, + w

∂ B 1 R  ≤ C p−q u p,S . ∂ν n 0, n

(3.18)

Denote z = PnB − N (PnB ; ρ , θ ). From the embedding theorem [2,15,25]:



∂q z ≤ C zk+q+ 1 ,S , k, q ≥ 0, 2 ∂ν q k,

we have

    ∂   −  ( ρ , θ ) H = z H ≤ z 0, + w  z  ∂ν 0, ≤ z 1 ,S + wz 3 ,S ≤ z 1 ,Smax + wz 3 ,Smax , 2 2 2 2 PnB

B N Pn ;

(3.19)

where Smax = {(r, θ )| r = rmax , 0 ≤ θ ≤ 2π }. For z = 0 with the Dirichlet condition u = f on  (or the Neumann condition uν = g on  ), there exists the bound (Oden and Reddy [25], and Babuska and Aziz [2])

z 12 +q,Smax ≤ C  f q,∂ Smax , (or z 12 +q,Smax ≤ C gq−1,∂ Smax ), q ≥ 0.

(3.20)

Also for (1.1)–(1.3), there exists the bound (see Grisvard [15])

  ∂  z 12 +q,Smax ≤ C zq,∂ Smax +  z ∂ν q−1,



∂ Smax

, q ≥ 0.

(3.21)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

355

Hence we obtain from (3.11) and (3.21)

     ∂  ∂  zH ≤ C z0,∂ Smax +  z + w z1,∂ Smax +  z ∂ν −1,∂ Smax ∂ν 0,∂ Smax

     ∂  ∂  ≤ C z0,∂ Smax + wz1,∂ Smax +  z + w z . ∂ν 0,∂ Smax ∂ν 0,∂ Smax

(3.22)

From Theorems 3.2 and 3.3, when w = 1/N we have for q = 0, 1,

zH ≤ CN

R 2n−N r n max rmax

rmin

PnB 0, .

(3.23)

Since

PnB 0, ≤ C u0, ≤ C u p,S ,

(3.24)

From (3.23) and (3.24), we have

PnB − N (PnB ; ρ , θ )H ≤ CN

R 2n−N r

max rmax

rmin

Combining (3.17), (3.18) and (3.25) gives





R 2n−N rmax n u − uN H ≤ C N + rmax

rmin

Hence from (3.14) we have

u − uN H

1 + Nn =O u p,S n p−1/2

 =O



u p,S .

(3.25)

 u p,S . p−1/2 1

n

(3.26)

 u  , p,S p−1/2 1

N

where we have used N ≤ Cn. This completes the proof of Theorem 3.4.

(3.27) 

4. Stability analysis for circular domains 4.1. Approaches for seeking eigenvalues Consider only the circular domains in this paper. For the non-circular domains, the stability analysis may follow [20,26]. Firstly let us cite the lemmas in [8]. Lemma 4.1. Let matrix A = Circulant(a0 , a1 , . . . , aN−1 ) with real entries ai . Then the eigenvalues of A are real, given by

λk =

N−1 

a j cos(k jh ),

k = 0, 1, . . . , N − 1,

j=0

where h = 2π /N. Also, the vector y0 =

√1 (1, 1, . . . , 1 )T N

is the eigenvector corresponding to the leading eigenvalue λ0 .

Denote the circular domain S◦ = {(, θ )|(0 ≤  ≤ ρ , 0 ≤ θ ≤ 2π )}. From (1.6), we have

uN = uN (  , θ ) =

N 

c j  j (, θ ) + d j φ j (, θ ),

j=1

where cj and dj are the unknown coefficients. Then we have 2N collocation equations on ∂ S:

u N ( ρ , θk ) =

N  j=1

w

c j  j ( ρ , θk ) +

N 

d j φ j (ρ , θk ) = f (ρ , θk ), k = 1, 2, . . . , N,

(4.1)

j=1

N N   ∂ ∂ ∂ uN ( ρ , θ ) = w cj  j ( ρ , θk ) + w dj φ j (ρ , θk ) = wg(ρ , θk ), k = 1, 2, . . . , N, ∂ρ ∂ρ ∂ρ j=1 j=1

(4.2)

where θk = kh, h = 2π /N, and w is a weight constant with w = 1/N. Hence the 2N coefficients cj and dj can be obtained by (4.1) and (4.2) if the system of linear algebraic equations is nonsingular. Denote (4.1) and (4.2) as the form of matrix and vector:

Ax = b,

(4.3)

356

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

where the vectors x = {c1 , . . . , cN , d1 , . . . , dN }T and b = { f1 , . . . , fN , wg1 , . . . , wgN }T , and the matrix A ∈ R2N × 2N is decomposed as



A=



A12 (φ ) . wA22 (Dφ )

A11 () wA21 (D)

(4.4)

In (4.4), the matrices A11 (), A12 (φ ), A21 (D), A22 (Dφ ) ∈ RN × N are given by

⎡

A11 () = ⎣

⎡ A12 (φ ) = ⎣

1 (θ1 ) .. .

... .. . ...

.. .

..

N (θ1 ) .. .

1 (θN ) φ1 ( θ 1 ) . . .

⎤

⎦,

N (θN ) ⎤ φN ( θ 1 ) .. ⎦,

(4.5)

.,

.

φ1 ( θ N ) . . . ⎡∂ ∂ρ 1 (θ1 ) ⎢ .. A21 (D) = ⎣ . ∂ ∂ρ 1 (θN ) ⎡∂ ∂ρ φ1 (θ1 ) ⎢ .. A22 (Dφ ) = ⎣ . ∂ ∂ρ φ1 (θN )

φN ( θ N ) ∂ ∂ρ N (θ1 )

... .. . ...

⎤

⎥ ⎦, ∂  (θ ) ∂ρ N N ⎤ ∂ φ (θ ) ∂ρ N 1 ⎥ .. ⎦. . ∂ φ (θ ) ∂ρ N N .. .

... .. . ...

The matrices A12 (φ ) and A22 (Dφ ) result from the Dirichlet and Neumann problems of Laplace’s equations on circular domains, given in [20,23], respectively. All four sub-matrices A11 (), A12 (φ ), A21 (D) and A22 (Dφ ) are circulant. Denote the eigen-matrix (see Davis [8], p.32)

⎡

1 1 ⎢ 1 F∗ (∈ C N×N ) = √ ⎢ . ⎣ N .. 1

... ... .. . ...

1

ω .. .

ωn−1

ω

⎤

1 n−1

.. .

ω

where ω = e2π i/N = cos 2π /N + i sin 2π /N with i =

⎥ ⎥, ⎦

(4.6)

√ −1. Based on Davis [8], we have

A11 () = F∗ 11 ()F, A12 (φ ) = F∗ 12 (φ )F,

(4.7)

A21 (D) = F∗ 21 (D)F, A22 (Dφ ) = F∗ 22 (Dφ )F, where F is unitary with FF∗ = F∗ F = I, and I is the identity matrix. In (4.7), the matrices 11 (), 12 (φ ), 21 (D) and 22 (Dφ ) are diagonal matrices, and denoted by

⎡

12 (φ ) = ⎣

λ0 (φ )

O ..

⎦,

.

λN−1 (φ )

O

⎡ λ0 (Dφ ) ⎣ 22 (Dφ ) = ⎡ 11 () = ⎣

..

(4.8)

⎤

O

⎦,

.

(4.9)

λN−1 (Dφ ) ⎤

O

λ0 ()

O

..

⎦,

.

λN−1 ()

O

⎡ λ0 (D) ⎣ 21 (D) =

O ..

⎤ ⎦.

.

λN−1 (D)

O We have from Li et al. [18,20,23]

λ0 (φ ) = N ln R + ε ,

⎤

ρ N ε≈− , R

(4.10)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

λk (φ ) ≈ −



N 1 ρ k 2

λ0 (Dφ ) = − λk (Dφ ) ≈ −

k

N

R

ρ N

+

N−k

,

,

R

357

k = 1, 2, . . . , N − 1, (4.11)

ρ R

k N ρ

2ρ

ρ N−k  1

+

R

ρ N−k 

,

R

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

The new eigenvalues λk () and λk (D) will be sought in the next subsection. By using matrix F∗ in (4.6), we obtain from (4.7)



A=

 =



A12 (φ ) wA22 (Dφ )

A11 () wA21 (D) F∗ O

O F∗



11 () w21 (D)

 12 (φ ) F w22 (Dφ ) O



O . F

(4.12)

Since the eigenvalues of similar matrices are the same, the eigenvalues of A are just those of matrix , denoted by



=

  12 (φ ) diag(λi ()) = w22 (Dφ ) w diag(λi (D))

11 () w21 (D)



diag(λi (φ )) . w diag(λi (Dφ ))

(4.13)

By a permutation transformation P, we have

⎡ λ () 0 ⎢wλ0 (D(φ )) ⎢ ⎢ ⎢ ⎢ ⎢ T P P = ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤

λ0 (φ )

wλ0 (Dφ ) ..

.

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥. ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

O ..

.

λi () wλi (D)

O

λi (φ ) wλi (Dφ ) ..

. ..

(4.14)

.

Hence the eigenvalues λ± of  can be found from the zero determinants: i

 λi (φ )  = 0, wλi (Dφ ) − λ

 λi () − λ wλi (D)

det 

i = 0, 1, 2, . . . , N − 1.

(4.15)

4.2. Eigenvalues λk () and λk (D) Firstly, consider the matrix in (4.5)

⎡

A11 () = ⎣

1 (θ1 ) .. .

1 (θN )

... .. . ...

N (θ1 ) .. .

⎤ ⎦.

(4.16)

N (θN )

Since i (θ j ) = i (ρ , θ j ), the matrix A11 () is circulant, denoted by

A11 () = circulant(b0 , b1 , . . . , bN−1 ),

(4.17)

where the entries are

b j = r 2j ln r j = (R2 + ρ 2 − 2Rρ cos jh ) ln



R2 + ρ 2 − 2Rρ cos jh.

(4.18)

Based on Lemma 4.1, the eigenvalues of A11 () are given by

λk =

N−1 

b j cos(k jh ),

k = 0, 1, . . . , N − 1.

(4.19)

j=0

Lemma 4.2. When ρ < R, there exist the approximations of λk () for large N:

λ0 () ≈ N (R2 ln R + ρ 2 (1 + ln R )) +

ρ N R 2 R

N−1

−

ρ2 N+1



,

(4.20)

358

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

N 2

λ1 () ≈ −

N λk () ≈ 2k

 ρ2 ρ

(1 + 2 ln R )R2 +

ρ k R 2 k−1

R

ρ N−k

R

N 2 (N − 1 )

ρ N−1 R2 N−2

R

−

ρ2 N

 ,

(4.21)



ρ2

−

N + 2 (N − k ) R

2

+

k+1



R2 ρ2 − , k = 2, 3, . . . , N − 1. N−k−1 N−k+1

(4.22)

Proof. Firstly for λ0 (), we have from (4.19)

λ0 () = where r j =



N−1 

bj =

j=0

N 1 2 hr j ln r j , h

(4.23)

j=1

R2 + ρ 2 − 2Rρ cos jh. Also from (2.15) there exists the expansion

(R2 + ρ 2 − 2Rρ cos x ) ln



R2 + ρ 2 − 2Rρ cos x

= R ln R + ρ (1 + ln R ) − (1 + 2 ln R )Rρ cos x ρ m ρ m ∞ ∞   1 1 + R2 cos mx − ρ 2 cos mx. m (m − 1 ) R m (m + 1 ) R 2

2

m=2

(4.24)

m=1

Then, from Lemma 2.1, we have

  1 2π 2 (R + ρ 2 − 2Rρ cos x ) ln R2 + ρ 2 − 2Rρ cos x dx + λ0 () h 0 1 = 2π (R2 ln R + ρ 2 (1 + ln R )) + λ0 (), h

λ0 () =

where the increments are obtained directly from (4.24) and (2.3)

2π λ0 () = h

≈N =



R2

 m=ν N 2

R N (N − 1 )

m

1 ρ m (m − 1 ) R

ρ N R

ρ N R 2

N−1

R

−

− ρ2



ρ2 N+1

m



− ρ2

m=ν N

1 ρ m (m + 1 ) R



ρ N 

1 N (N + 1 ) R

.

Combining (4.25) and (4.24) gives

λ0 () ≈ N{R2 ln R + ρ 2 (1 + ln R )} +

(4.25)

(4.26)

ρ N R 2 N−1

R

−



ρ2 N+1

.

(4.27)

This is the first result (4.20). Next for λ1 (), we have from Lemma 2.1, (4.19) and (4.24),

λ1 () =

N−1 

b j cos jh =

j=1

N N 1 1 2 hb j cos jh = hr j ln r j cos jh h h j=1

j=1

  1 2π 2 = (R + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x × cos xdx + λ1 () h 0 =

π

− (1 + 2 ln R )Rρ −

h

=−

"

ρ3 ! 2R

+ λ1 ()

#

N ρ2 ρ ((1 + 2 ln R )R2 + )( ) + λ1 (). 2 2 R

(4.28)

Also the increments are given directly from Lemma 2.2 and (4.24),

λ1 () =

π h

R2

 m=±1+ν N

m

ρ 1 m (m − 1 ) R

−

ρ2

 m=±1+ν N

m

ρ 1 m (m + 1 ) R

N 2 1 ρ 1 ρ ≈ R ( )N−1 − ρ 2 ( )N−1 2 (N − 1 )(N − 2 ) R N (N − 1 ) R

!

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

=

"

359

#

N ρ R2 ρ2 ( )N−1 − . 2 (N − 1 ) R (N − 2 ) N

(4.29)

Combining (4.26) and (4.29) gives the second result (4.21). Finally, for λk () (k = 2, . . . , N − 1 ), we have similarly

λk () =

N 1 hb j cos(k jh ) h j=1

  1 2π 2 = (R + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x cos kxdx + λk () h 0 = =

π

k

ρ 1 R k (k − 1 ) R 2

h

R2 N ρ k

2k

k−1

R

−

−ρ

2



ρ2 k+1

k 

ρ 1 k (k + 1 ) R

+ λk ()

+ λk (),

(4.30)

where λk () is obtained from Lemma 2.2 and (4.24)



π

λk () =

R2

h

N ≈ 2

m

 m=±k+ν N

ρ 1 m (m − 1 ) R

N−k

R2 ρ (N − k )(N − k − 1 ) R

N−k

N ρ = 2 (N − k ) R



− ρ2

m=±k+ν N

m

ρ 1 m (m + 1 ) R

ρ N−k ρ2 − (N − k )(N − k + 1 ) R  ρ2

(4.31)



R2 − . N−k−1 N−k+1

Combining (4.30) and (4.31) gives the desired result (4.22), and completes the proof of Lemma 4.2.



∂ (r 2 ln r ). The following lemma follows from (4.24). We need the expansions of the normal derivatives ∂ρ

Lemma 4.3. For ρ < R, there exists the expansion,

 ∂ 2 (R + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x ∂ρ = 2ρ (1 + ln R ) − (1 + 2 ln R )R cos x + R

∞  m=2

The matrix in (4.5) is expressed by

⎡ ⎢

A21 (D) = ⎣

∂ ∂ρ 1 (θ1 )

.. .

∂ ∂ρ (θN )

... .. . ...

∂ ∂ρ N (θ1 )

1 m−1

ρ m−1 R

cos mx − ρ

∞  m=1

m

ρ m+2 m (m + 1 ) R

cos mx.

(4.32)

⎤

⎥ ⎦, ∂  (θ ) ∂ρ N N .. .

(4.33)

which is also circulant:

A21 (D) = circulant (d0 , d1 , . . . , dN−1 ), where

dj =

 ∂ 2 (R + ρ 2 − 2Rρ cos jh ) ln R2 + ρ 2 − 2Rρ cos jh . ∂ρ

(4.34)

Based on Lemma 4.1, the eigenvalues of A21 (D) are given by

λk =

N−1 

d j cos k jh, k = 0, 1, . . . , N − 1.

(4.35)

j=0

Lemma 4.4. When ρ < R, there exist the approximations of λk (D) for large N

N λ0 (D) ≈ 2Nρ (1 + ln R ) + N−1

λ1 (D) ≈ −



ρ N−1 1

 (N + 2 )(N − 1 ) 2 R − ρ , R ( N + 1 )N

R



2

N 3 N (1 + 2 ln R )R2 + ρ 2 + 2R 2 2 (N − 2 )

ρ N−2 1 R

R

R2 −

 (N + 1 )(N − 2 ) 2 ρ , N (N − 1 )

(4.36)

(4.37)

360

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

λk (D) ≈

 (k + 2 )(k − 1 ) 2 ρ R R k (k + 1 )



N−k−1 ρ N 1 (N − k + 2 )(N − k + 1 ) 2 + R2 − ρ , 2 (N − k − 1 ) R R (N − k )(N − R + 1 ) N 2 (k − 1 )

ρ k−1 1

R2 −

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

(4.38)

Proof. Firstly for λ0 (D), we have from (4.32) similarly

λ0 (D) =

N−1 

dj =

j=0

N 1 hd j h j=1

 !  1 2π ∂ = (R2 + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x dx + λ0 (D) h 0 ∂ρ 2π = 2ρ (1 + ln R ) + λ0 (D) = 2Nρ (1 + ln R ) + λ0 (D), h

(4.39)

where

  2π 1 ρ m+2 ρ R ( )m−1 − ρ ( )( )m h m−1 R m (m + 1 ) R

λ0 (D) =

m=ν N

m=ν N

!

1 ρ N+2 ρ ≈N R ( )N−1 − ρ ( )N N−1 R N (N + 1 ) R " N ρ 1 (N + 2 )(N − 1 ) 2 # = ( )N−1 R2 − ρ . N−1 R R N (N + 1 )

(4.40)

Combining (4.39) and (4.40) gives the first result (4.36). Next for λ1 (D), we have from Lemma 2.2, (4.34) and (4.35),

λ1 (D) =

N−1 

N 1 hd j cos jh h

d j cos jh =

j=0

(4.41)

j=1

!  ∂ (R2 + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x cos xdx + λ1 (D) ∂ρ 0 π 3 ρ = {−(1 + 2 ln R )R − ρ ( )} + λ1 (D), =

1 h

 2π

h

2 R

where λ1 (D) is obtained from Lemma 2.2 and (4.32),

(m + 2 ) ρ m ! ( ) h m (m + 1 ) R m=±1+ν N m=±1+ν N ! N R ρ N+1 ρ ≈ ( )N−2 − ρ ( )N−1 2 N−2 R N (N − 1 ) R N ρ N−2 1 2 (N + 1 )(N − 2 ) 2 ! = ( ) R − ρ . 2 (N − 2 ) R R N (N − 1 )

λ1 (D) =

π

R



1 ρ ( )m−1 − ρ m−1 R



(4.42)

Combining (4.41) and (4.42) gives the second result (4.37). Finally, for λk (D) (k = 2, 3, . . . , N − 1 ), we have similarly

λk (D) =

N 1 hd j cos k jh h j=1

!  ∂ (R2 + ρ 2 − 2Rρ cos x ) R2 + ρ 2 − 2Rρ cos x cos kxdx + λk (D) ∂ρ 0 ! π R ρ k+2 ρ k = ( )k−1 − ρ ( ) + λk (D) h k−1 R k (k + 1 ) R N ρ k−1 " 2 (k + 2 )(k − 1 ) 2 # = ( ) R − ρ + λk (D), k = 2, 3, . . . , N − 1, 2 ( k − 1 )R R k (k + 1 ) 1 = h

 2π

where

λk (D) =

π h

 m=±k+ν N

R ρ ( )m−1 − ρ m−1 R

 m=±k+ν N

m+2 ρ ( )m m (m + 1 ) R

(4.43)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

361

!

N R ρ ρ (N − k + 2 ) ρ ( )N−k−1 − ( )N−k 2 N−k−1 R (N − k )(N − k + 1 ) R " N ρ 1 (N − k + 2 )(N − k − 1 ) 2 # = ( )N−k−1 × R2 − ρ . 2 (N − k − 1 ) R R (N − k )(N − k + 1 ) ≈

Combining (4.43) and (4.44) gives the last result (4.38), and this completes the proof of Lemma 4.4.

(4.44) 

4.3. Bounds of condition number Since the eigenvalues λj (φ ), λj (), λj (Dφ ) and λj (D) are provided in (4.11) and Lemmas 4.2 and 4.4, we can derive the approximate eigenvalues of the MFS. Lemma 4.5. Let ρ < R and w = 1/N. When R = 1 and R = 1/e, there exist the asymptotes of the leading eigenvalues1

λ+0  O(N ), λ−0  O(1 ). Proof. From (4.14), the



A0F

λ0 (φ ) = wλ0 (Dφ )

λ±0

(4.45)

result from the eigenvalues of the 2 × 2 matrix

  N ln R λ0 () ≈ − ρ1 ( ρR )N wλ0 (D)



N[R2 ln R + ρ 2 (1 + ln R )] . 2ρ (1 + ln R )

(4.46)

The term − ρ1 ( ρR )N ≈ 0 for ρ < R and large N. When R = 1 and R = 1/e, the eigenvalues of A0F are given by

λ+0 ≈ N ln R, λ−0 ≈ 2ρ (1 + ln R ).

(4.47) 

Lemma 4.6. Let ρ < R and w = 1/N. For even N there exist asymptotes:

λ+1  O(N ), λ−1  O(1 ), λ+k  O

k  N ρ k

,

R

λ±N−k = λ±k ,

(4.48)

λ−k  O

k = 2, 3, . . . ,

1 k2

ρ k 

, k = 2, . . . ,

R

N , 2

(4.49)

N − 1. 2

(4.50)

Proof. From (4.11) and Lemmas 4.2 and 4.4, we have the matrix for λ± 1 (1 )

AF



λ1 (φ ) = wλ1 (Dφ )  N ρ =−

2

 λ1 () wλ1 (D)

1

1 Nρ

R

 2 (1 + 2 ln R )R2 + ρ2   1 (1 + 2 ln R )R2 + 32 ρ 2 Nρ

= cNB(1) , where

 B (1 ) =

(4.51)



1ρ a2 ,c = − a4 2R

a1 a3

(4.52)

ρ2

a2 = (1 + 2 ln R )R2 +

a1 = 1,

2

,



1 1 3 a3 = , a4 = (1 + 2 ln R )R2 + ρ 2 . Nρ Nρ 2

(4.53)

Firstly, the determinant of B(1) is given by

det(B(1) ) = a1 a4 − a2 a3 = =

1

1 Nρ

ρ

N



 3 ρ2 (1 + 2 ln R )R2 + ρ 2 − (1 + 2 ln R )R2 + 2

2

> 0.

The notation a≈O(b) denotes that there exist two constants C1 and C2 such that C1 b ≤ |a| ≤ C2 b, b > 0.

(4.54)

362

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

Denote λ ± of eigenvalues of B(1) in (4.52), we have





 1 ( a1 + a4 ) ± ( a1 − a4 )2 + 4 a2 a3 . 2

λ± =

(4.55)

When N is large

1 N

1 N

|a4 | ≤ C ,

a1 = 1 ,

|a2 a3 | ≤ C ,

(4.56)

we have

a1  |a4 |,

4 |a2 a3 |  ( a1 − a4 )2 .

λ+1

We conclude that

λ−1

>

> 0, to obtain from (4.55)

det(B )

λ+1  O(1 ), λ−1 =

(4.57)

O

λ+

1

N

.

(4.58)

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

(4.59)

This is the first result (4.48). Next, we examine the matrix (k )

AF



 λk () , wλk (D)

λk (φ ) = wλk (Dφ )

The entries of AF(k ) in (4.59) are given from (4.11) and Lemmas 4.2 and 4.4,

λk (φ ) ≈ −

N 2

λk (Dφ ) ≈ − N λk () ≈ 2k ≈

αk

k 1 ρ k

N 2ρ

R

1 N−k

+

k ρ

+

R

ρ k R 2 k−1

R

R −

k−1

k

k−1 1

N ρ 2k − 1 R

R

ρ N−k 

R2 N ρ k 2k

ρ N−k 



ρ2 k+1 −

ρ2

N1 2k

≈ −αk

≈ −αk

N 2ρ

ρ k

ρ k R

R

,

,

N−k

N ρ + 2 (N − k ) R



(4.61) R2 ρ2 − N−k−1 N−k+1

,

k+1

where αk = 1 for k =

R

(4.62)

AF(N−k ) = AF(k ) ,

R2 −

and αk = 2 for k =

N 2

k = 1, 2, . . . ,

and λ± = λ± , k = 1, 2, . . . , N−k k Denote the matrix

AF(k ) = αk

N 2

ρ k R

N 2

N 2.

1 =− , k

N , 2

(4.64)

− 1.

× B∗ , B∗ =



a∗1 a∗3

a∗2 , a∗4

a∗3

1 =− , Nρ

a∗2

(4.65)



k−1 2 1 = R2 − ρ , k (k − 1 ) k+1

a∗4

1 1 = N (k − 1 ) ρ

 (k + 2 )(k − 1 ) 2 R − ρ . k (k + 1 )

2

)(k−1 ) Since (k+2 < 1, the bracket of a∗4 in (4.66) is positive for ρ < R, to give k(k+1 )

R2 −

(4.63)

For even N there exist equalities

where the entries are given from (4.59)–(4.63),

a∗1





(k + 2 )(k + 1 ) 2 ρ k (k + 1 )

 ρ N−k−1 1 N (N − k + 2 )(N − k + 1 ) 2 2 + R − ρ 2 (N − k − 1 ) R k (N − k )(k − N + 1 )  ρ k−1 1 N (k + 2 )(k − 1 ) 2 ≈ αk R2 − ρ , 2 (k − 1 ) R R k (k + 1 )

λk (D) ≈

(4.60)

(k + 2 )(k − 1 ) 2 ρ  O ( 1 ). k (k + 1 )

(4.66)

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

363

Hence we conclude a∗1 < 0, a∗2 > 0, a∗3 < 0 and a∗4 > 0. The determinant of AF is given by

det(B∗ ) = a∗1 a∗4 − a∗2 a∗3 =− =

1 1 (k + 2 )(k − 1 ) 2 R2 − ρ N ρ k (k − 1 ) k (k + 1 )





+



k−1 2 1 1 R2 − ρ N ρ k (k − 1 ) k+1



(4.67)

2ρ 1 O . Nk2 (k + 1 ) Nk3

Next, we obtain the eigenvalues of B∗ from (4.66)

 λ±∗ = (a∗1 + a∗4 ) ± (a∗1 − a∗4 )2 + 4a∗2 a∗3 $ 1

1

1

1

2 =O

+O

k

±

Nk

O

k

Similarly, we have

λ+∗  O

1

k

λ−∗ =

,

det(B∗ )

λ

+ ∗

+O

O

1

Nk2

Nk

+O

1

Nk2

.

(4.68)

.

(4.69)

Then the eigenvalues of AF(k ) have the asymptotes from (4.65) and (4.69)

λ+k  O

k  N ρ k

,

R

λ−k  O

1 k2

ρ k  R

.

(4.70)

The desired results (4.49) follow. This completes the proof of Lemma 4.6.



Theorem 4.1. Suppose that ρ < R, w = 1/N and even N(  1). When R = 1 and R = 1/e, there exists the asymptote

N  R 2 Cond = O N 3 . ρ

(4.71)

Proof. From Lemmas 4.5 and 4.6, we obtain

λmax |λ±i | = |λ+0 |  O(N ), λmin |λ±i | = |λ−N |  O 2

1 N2

ρ N/2  R

,

(4.72)

to give

N/2  λmax |λ±i | O(N ) R Cond = = = O N3 . ρ λmin |λ±i | |λ−N |

(4.73)

2

This is the desired result (4.71), and completes the proof of Theorem 4.1.



When N is odd, we can prove Cond = O(N 3 (R/ρ )N/2 ), where x is the maximal integer ≤ x. Note that Cond in (4.71) is similar to Cond = O(N (R/ρ )N/2 ) for Laplace’s equation in [20]. By the above arguments and [20,26], the bounds of Cond can be derived for the simply-connected domains, and extended for other boundary conditions, such as the simply support condition, and the mixed type of different boundary condition. Details are omitted. 5. Numerical experiments Consider the biharmonic Eq. (1.1) in the rectangular domain S = {(x, y )| − 1 < x < 1, −1 < y < 1}, and we choose the following exact solution

u(x, y ) = exp(x ) cos y + (x2 + y2 ) exp(y ) cos x.

(5.1)

Consider the mixed type of the clamped and simply support boundary conditions on  , then Eq. (1.3) is replaced by

uν = g on 1 , uνν = g∗ on 2 ,

(5.2)

where 1 ∪ 2 =  , and 1 ∩ 2 = ∅. The admissible functions (1.6) are expressed as

uN = uN ( ρ , θ ) =

N  j=1

=

N  j=1

c j (r 2j ln r j ) +

N 

d j ln r j

j=1

(c j  j (ρ , θ ) + d j φ j (ρ , θ )),

(5.3)

364

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366 Table 1 Errors and condition numbers for the biharmonic equation with the clamped boundary condition by the MFS, where  = u − uN .

R=2

R=5

N

M

B

0, S

1, S

Cond

Cond_eff

x

11 21 31 41

5 10 15 20

2.25(−1) 2.16(−3) 2.89(−5) 6,03(−7)

8.77(−2) 4.46(−4) 4.33(−6) 3.72(−8)

6.08(−1) 6.89(−3) 9.07(−5) 1.94(−6)

1.54(4) 1.28(6) 4.10(7) 1.32(9)

86.0 5.95(3) 1.89(5) 6.24(6)

27.3 23.6 19.6 16.6

11 21 31

5 10 15

2.86(−2) 1.46(−6) 3.10(−11)

1.11(−2) 5.35(−7) 1.71(−11)

8.73(−2) 5.89(−6) 8.48(−11)

2.22(7) 1.64(11) 4.41(14)

1.69(2) 6.68(5) 1.79(9)

2.33(3) 3.14(3) 2.59(3)

where cj and dj are the unknown coefficients, and r j = |P Q j |. Choose the uniform source nodes Q j = (R, ϕ j ), where ϕ j = 2π j/N. The energy I (v ) in (1.7) is replaced by

I ∗ (v ) =





( v − f )2 + w2



1

(vν − g)2 + (w∗1 )2



2

(vνν − g∗ )2 ,

(5.4)

where the weights w = 1/N and w∗1 = 1/N 2 as in [24]. Next, we formulate the collocation equations directly from (1.2) and (5.2), and obtain the collocation equations at nodes (ρ k , θ k ) ∈  for the mixed type of boundary conditions:

u N ( ρk , θ k ) =

N  



c j  j ( ρk , θ k ) + d j φ j ( ρk , θ k ) = f ( ρk , θ k ) , ( ρ k , θ k ) ∈  ,

(5.5)

j=1

 N  ∂ ∂ ∂ u N ( ρk , θ k ) = w cj  j (ρk , θk ) + d j φ j (ρk , θk ) = wg(ρk , θk ), (ρk , θk ) ∈ 1 , ∂ν ∂ν ∂ν j=1

 N  ∂2 ∂2 ∂2 w 2 2 u N ( ρk , θ k ) = w 2 cj  ( ρ , θ ) + d φ ( ρ , θ ) = w2 g∗ (ρk , θk ), (ρk , θk ) ∈ 2 , j j j k k k k 2 2 ∂ν ∂ν ∂ν j=1 w

(5.6)

(5.7)

where ν is the normal of  1 and  2 , and w = 1/N. For the collocation equations (5.5)– (5.7), the number of collocation nodes is often chosen to be larger than the number of source nodes. Then we obtain the over-determined system of linear algebraic equations

Fx = b, where matrix by

Cond =

(5.8)

F ∈ Rm × n (m ≥ n),

x ∈ Rn

and

b ∈ Rm .

The traditional and the effective condition numbers in [21] are defined

σmax b , Cond_eff = , σmin σmin x

(5.9)

where σ max and σ min are the maximal and the minimal singular values of the matrix F, respectively, and x is the Euclidean norm. The boundary errors are given by

B = u − uN B =





I (uN )or I∗ (uN ),

where I (v ) and I∗ (v ) are defined in (1.8) and (5.4), respectively. When 2 = ∅, the mixed type is reduced to the purely clamped boundary conditions. For the clamped boundary conditions, the errors and condition numbers are listed in Table 1, where M denotes the number of uniform collocation nodes along each edge of ∂ S. Then n = 2N and m = 8M in (5.8). From Table 1, we can see the asymptotes for R = 2

B = O(0.65N ), 0,S = O(0.66N ), 1,S = O(0.61N ),

(5.10)

Cond = O(1.4N ), Cond_eff = O(1.3N ),

(5.11)

and for R = 5

B = O(0.37N ), 0,S = O(0.37N ), 1,S = O(0.35N ),

(5.12)

Cond = O(2.3N ), Cond_eff = O(2.2N ).

(5.13)

The errors decrease exponentially via the number of FS used; the exponential convergence ratesresult from highly √ the √ smooth solution (5.1). On the the hand, the condition numbers grow also exponentially. By noting R/rmin = R = 2 (or

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

365

Table 2 Errors and condition numbers for the biharmonic equation with the mixed type of the clamped and simply supported boundary conditions by the MFS, where  = u − uN .

R=2

R=5

N

M

B

0, S

1, S

Cond

Cond_eff

x

11 21 31 41

5 10 15 20

1.87(−1) 2.01(−3) 2.62(−5) 5.77(−7)

2.27(−1) 1.05(−3) 6.33(−6) 5.71(−8)

1.07 8.23(−3) 1.08(−4) 2.32(−6)

1.83(4) 1.46(6) 4.87(7) 1.52(9)

1.37(2) 6.48(3) 2.24(5) 7.20(6)

20.3 24.1 19.6 16.6

11 21 31

5 10 15

2.48(−2) 1.28(−6) 4.87(−11)

1.19(−2) 7.39(−7) 3.59(−11)

8.91(−2) 6.93(−6) 1.28(−10)

3.35(7) 2.00(11) 4.72(14)

2.59(2) 8.14(5) 1.92(9)

2.28(3) 3.14(3) 2.59(3)

√

5) the Cond in (5.11) and (5.13) coincides with (4.71) very well. Let us scrutinize (4.71) more carefully. Choose R = 5, and find the numerical growth of Cond between N1 = 31 and N1 = 11 from Table 1

 

T1 =

Cond

4.41(14 ) N1 =31 = = 1.94(7 ).  2.27(7 ) Cond

(5.14)

N1 =11

On the other hand, the theoretical growth of Cond is evaluated from (4.71) by

T2 =

N 3 1

N2

×R

N1 −N2 2

=

31 3 11

×5

31−11 2

= (2.82 )3 × 510 = 2.19(7 ).

(5.15)

The relative error of T1 is given by

|T1 − T2 | T2

=

|1.94(7 ) − 2.19(7 )| = 11.5%, 2.19(7 )

(5.16)

to show the remarkable consistency between theory and computation of stability. The effective condition number is a better criteria for stability in numerical PDE, see [21]. We obtain the ratio from Table 1 at R = 5 and N = 31

Cond 4.41(14 ) = = 2.46(5 ), Cond_eff 1.79(9 )

(5.17)

to show that the effective condition number is much smaller than the traditional condition number. The bounds of the effective condition may be derived by following [21]; details are omitted. Next, we still choose the true solution (5.1), but use the mixed type of the clamped and simply supported boundary conditions:

u = f, uν = g, on x = ±1, u = f, uνν = g∗ , on y = ±1,

(5.18)

where ν is the exterior normal of ∂ S. The errors and condition numbers of the MFS are listed in Tables 2. The asymptotes of errors and conditions are similar to (5.10)–(5.13). 6. Conclusions To close this paper, let us summarize the new results in this paper. 1. The error and stability analysis are explored for biharmonic equations by the MFS using r2 ln r2 . The error bounds are derived as

u − uN H = O



1 N p−1/2



u p,S ,

to show the polynomial convergence rates. This paper is the advanced study of [22], where a preliminary analysis was made for the MFS using simple FS: ρ 2 ln r2 . 2. For the circular domains, the condition number grows exponentially, as Cond O(N 3 (R/ρ )N/2 ). For the bounded simplyconnected domains and under other boundary conditions, similar bounds of condition number can also be derived by following Ref. [20,26]. Details are omitted. 3. Numerical experimenters are carried out for the simply support conditions and the mixed type of boundary conditions, to coincide the analysis made in this paper.

366

F. Dou et al. / Applied Mathematics and Computation 339 (2018) 346–366

Acknowledgements The authors would like to thank Professor J.T. Chen and reviewers for valuable comments and suggestions, and also indebted to Mr. T. C. Wu for assisting the computation in Section 5. The third author acknowledges HPC at The University of Southern Mississippi supported by the National Science Foundation under the Major Research Instrumentation (MRI) program via Grant # ACI 1626217. The fourth author acknowledges the support of the China Scholarship Council under Grant # 201706935029. The first author thanks the support of the National Natural Science Foundation of China (No. 11501086) and the Fundamental Research Funds for the Central Universities (No. ZYGX2016J137). References [1] E.A. Atkinson, An Introduction to Numerical Analysis (Sec. Ed.), John Wiley & Sons, New York, 1989, p. 178. [2] 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, p. 32. 3–359. [3] Y. Bai, Y. Wu, X. Xie, Uniform convergence analysis of a higher order hybrid stress quadrilateral finite element method for linear elasticity, Adv. Appl. Math. Mech. 8 (3) (2016) 399–425. [4] Y.H. Bai, Y.K. Wu, X.P. Xie, Superconvergence and recovery type a posteriori error estimation for hybrid stress finite element method, Sci. China Math. 59 (2016) 1835–1850. [5] A. Bogomolny, Fundamental solutions method for elliptic boundary value problems, SIAM J. Numer. Anal. 22 (1985) 644–669. [6] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zhang, Spectral methods, in: Fundamentals in Single Domains, Springer, New York, 2006. [7] J.T. Chen, C.S. Wu, Y.T. Lee, K.H. Chen, On equivalence of Trefftz method of fundamental solutions for Laplace and biharmonic equations, Comp. Math. Appl. 53 (2007) 851–879. [8] P. Davis, Circulant Matrices, John Wiley & Sons, New York, 1979. [9] G. Fairweather, A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Adv. Comp. Math. 9 (1998) 69–95. [10] A. Ge, M. Feng, Y. He, Stabilized multiscale finite element method for the stationary Navier−Stokes equations, J. Math. Anal. Appl. 354 (2) (2009) 708–717. [11] 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. [12] A. Karageorghis, G. Fairweather, The Almansi method of fundamental solutions for solving biharmonic problems, Int. J. Numer. Methods Eng. 26 (1988) 1665–1682. [13] A. Karageorghis, G. Fairweather, The method of fundamental solutions for the numerical solution of biharmonic equation, J. Comp. Phys. 69 (1998) 434–459. [14] V.D. Kupradze, Potential methods in elasticity, in: J.N. Sneddon, R. Hill (Eds.), Progress in Solid Mechanics, III, Amsterdam, 1963, pp. 1–259. [15] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Advanced Publishing Program, Boston and London, 1985. 328–338 [16] Y. Jiang, B. Wang, Y. Xu, A fast Fourier−Galerkin method solving a boundary integral equation for the biharmonic equation, SIAM J. Numer. Anal. 52 (5) (2014) 2530–2554. [17] J. Le, H. Jin, X.G. Lv, Q.S. Cheng, A preconditioned method for the solution of the Robbins problem for the Helmholtz equation, Anziam J. 52 (2010) 87–100. [18] Z.C. Li, The Method of fundamental solutions for annular shaped domains, J. Comp. Appl. Math. 228 (2009) 355–372. [19] Z.C. Li, Analysis of Method of Fundamental Solution for Biharmonic Equations, Technical Report, Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, 2012. [20] Z.C. Li, J. Huang, H.T. Huang, Stability analysis of method of fundamental solutions for mixed boundary value problems of Laplace’s equations, Computing 88 (2010) 1–29. [21] Z.C. Li, H.T. Huang, Y. Wei, A.H.D. Cheng, in: Effective Condition Number for Partial Differential Equations, second ed., Science Press, Beijing, China, 2015. [22] Z.C. Li, M.G. Lee, J.Y. Chiang, Y.P. Liu, The Trefftz method using fundamental solutions for biharmonic equations, J. Comp. Appl. Math. 235 (2011) 4350–4367. [23] Z.C. Li, M.G. Lee, H.T. Huang, J.Y. Chiang, Neumann problems of 2D Laplace’s equation by method of fundamental solutions, Appl. Numer. Math. 119 (2017) 126–145. [24] Z.C. Li, T.T. Lu, H.Y. Hu, A.H.D. Cheng, Trefftz and Collocation Methods, WIT, Southsampton, Boston, 2008. [25] J.T. Oden, J.N. Reddy, An Introduction to the Mathematical Theory of Finite Elements, John Wiley & Sons, New York, 1976, pp. 182–196. [26] Z. Tian, Z.C. Li, H.T. Huang, C.S. Chen, The method of fundamental solutions for the modified Helmholtz equation, Appl. Math. Comput. 305 (2017) 262–281. [27] E. Trefftz, Konvergenz und ritz’schen verfahren, in: Proceedings of the Second International Congress of Applied Mechanics, Zurich, 1926, pp. 131–137. [28] H.Y. Wu, Y. Duan, Multi-quadric quasi-interpolation method coupled with FDM for the Degasperis−Procesi equation, Appl. Math. Comput. 274 (2016) 83–92. [29] C. Wang, T.Z. Huang, C. Wen, A new preconditioner for indefinite and asymmetric matrices, Appl. Math. Comput. 219 (2013) 11036–11043.