Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators S.Yu. Reutsky 1 Abstract: In ths paper a new meshless method for egenproblems wth Laplace and bharmonc operators n smply and multply connected domans s presented. The soluton of an egenvalue problem s reduced to a sequence of nhomogeneous problems wth the dfferental operator studed. These problems are solved usng the method of fundamental solutons. The method presented shows a hgh precson n smply and multply connected domans. The results of the numercal experments justfyng the method are presented. keyword: Method of fundamental solutons, Membranes and Plates, Free vbraton problem 1 Introducton The goal of ths paper s to present a new numercal technque for soluton of the followng egenvalue problems: 2 w+k 2 w = 0, x Ω R 2, (1) B 1 [w]=0, x Ω and 4 w k 4 w = 0, x Ω R 2, (2) w = 0, B 2 [w]=0, x Ω. Here Ω s a smply or multply connected doman of nterest wth boundary Ω. The boundary operator n (1) B 1 [...] wll be consdered of the two types: the Drchlet B 1 [w]=w and of the eumann type B 1 [w]= w/ n; for bharmonc operator n (2), B 2 [w] = w/ n or B 2 [w] = 2 w/ n 2. As a mechancal applcaton, ths corresponds to recoverng the free vbraton frequences of membranes and plates. Such problems often arse n engneerng applcatons. The usual approach for egenvalue problems wth a selfadjont operator s to use the Raylegh mnmal prncple. 1 Laboratory of Magnetohydrodynamcs, Tmurovtzev, 29 D, ap. 51, 61142, Kharkov, Ukrane In partcular, the statonary ponts of the functonal Z Z R(w)= w 2 dω/ w 2 dω. Ω concde wth egenfunctons of the Laplace operator. See [Courant (1943); Courant and Hlbert (1953); Morse and Feshbach (1953); Strang and Fx (1973)] for more detals and references. Then, usng an approxmaton for w wth fnte number of free parameters, one gets the same problem n a fnte-dmensonal subspace whch can be solved by a standard procedure of lnear algebra, e.g., see [Golub and Loan (1996); Strang (1976)]. However, a standard fnte dfferences method can produce good results when dealng wth a partcular type of shapes defned on rectangular grds, whle for other type of shapes the fnte element method or the boundary element method are more approprated. The method of fundamental solutons (MFS) [Farweather and Karageorghs (1998); Golberg and Chen (1998, 1997)] s the convenent tool n ths feld. The smlar technque s used n the boundary knot method (BKM) [Chen (2005); Chen and Tanaka (2002)]. Unlke the MFS, t employs nonsngular general solutons as the bass functons to avod the fcttous boundary outsde the physcal doman. In the framework of the boundary methods a general approach to solvng these problems s as follows. Frst, usng an ntegral representaton of w n the BEM, or an approxmaton over fundamental solutons n MFS, one gets a homogeneous lnear system A (k)q = 0 wth matrx elements dependng on the wave number k. The determnant of ths matrx must be zero to obtan the nontrval soluton: det[a (k)] = 0 (3) Ths equaton must be nvestgated analytcally or numercally to get the egenvalues. Ths technque s descrbed n [Karageorghs (2001); Chen, Ln, Kuo, and Chyuan (2001); Chen, Lu, and Hong (2003); Chen, Chen, Chen,
178 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 Lee, and Yeh (2004); Chen, Chen, and Lee (2005)] wth more detals. In the two latest papers there s a complete bblography on the subject consdered. In [Alves and Antunes (2005)] the MFS based procedure s appled to egenproblems n general smply connected shapes. The method presented n ths artcle uses the same MFS boundary technque. Ths s a mathematcal model of physcal measurements when the resonance frequences of a system are determned by the ampltude of response to some external exctaton. As a result, e.g., nstead of (1) we solve a sequence of nhomogeneous boundary value problems (BVP): 2 w+k 2 w = f (x), x Ω R 2, (4) B 1 [w]=0, x Ω, where f descrbes some source placed outsde the soluton doman. Let F (k) be some norm of the soluton w. As t wll be shown below, ths functon of k has sharp maxmums at the egenvalues and, under some condtons can be used for ther determnng. Certanly such behavour of F (k) near the egenvalues s a consequence of (3). Technques of numercal soluton of lnear BVPs lke (4) are well developed. It should be emphaszed that any Helmholtz (or bharmonc) equaton solver can be used n the framework of the method presented. However, the MFS technque seems to be a more sutable one for ths goal n the case of an arbtrary doman. The outlne of ths paper s as follows: for the sake of smplcty we begn by descrbng the 1D case n Secton 2. In Secton 3, we present the algorthm of MFS n applcaton to problem (1). Here we present numercal examples to llustrate the method presented for smply and multple connected domans. In Secton 4, the same technque s descrbed n applcaton to problem (2). Some generalzaton of the technque and the felds of ts development are dscussed n Secton 5. 2 One-dmensonal egenproblem To llustrate the method presented let us consder the wave equaton [Morse and Feshbach (1953)] 2 u t 2 = 2 u x 2 (5) wth the Drchlet condtons at the endponts of the nterval [0,1],.e., u(0,t) =u(1,t)=0. Consderng the free harmonc vbratons u(x,t)=e kt w(x), wegetthe followng 1D Sturm Louvlle problem on the nterval [0,1]: d 2 w dx 2 +k2 w = 0, w(0)=w(1)=0. (6) The well known soluton s: k n = nπ, w n = sn(nπx), n = 1,2,...,. Followng the boundary approach, let us consder the fundamental soluton Ψ(x,ξ,k)= 1 exp (k x ξ ), (7) 2k whch satsfes the homogeneous equaton everywhere except the sngular pont x = ξ. A general soluton of the homogeneous equaton n the nterval [0,1] can be wrtten n the form: w = q 1 Ψ(x,ξ 1,k)+q 2 Ψ(x,ξ 2,k). Here ξ 1,ξ 2 are two source ponts placed outsde the soluton doman [0,1], e.g., ξ 1 < 0, ξ 2 > 1; q 1,q 2 are free parameters. Usng the boundary condtons w(0)=w(1)= 0, one gets the lnear system: { q1 e A (k)q = ( kξ1) +q 2 e (kξ2) = 0 q 1 e ( k(1 ξ1)) +q 2 e (k(ξ2 1)) = 0 The wave numbers k n can be determned from the condton: det[a (k)] = 0. After smple transforms we get: exp (2k)=1, or k = nπ. Thus, MFS gets the exact soluton. ote that n multdmensonal cases such computatons are tme consumng and not so smple. As t s mentoned above, the method suggested s a mathematcal model of physcal measurements, when a mechancal or acoustc system s excted by an external source and resonance frequences can be determned usng an ncrease of ampltude of oscllatons near these frequences. So, nstead of (6) we solve the nhomogeneous problem: d 2 w dx 2 +k2 w = f (x), w(0)=w(1)=0. (8) The general soluton can be wrtten n the form: w = q 1 Ψ(x,ξ 1,k)+q 2 Ψ(x,ξ 2,k)+w p. (9) When the exctaton s performed by the pont source wth the same wave number k whch s placed at the pont
The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators 179 ξ 0 outsde the soluton doman, then f (x) =δ(x ξ 0 ) and the partcular soluton s: w p = Ψ(x,ξ 0,k)= 1 2k exp (k x ξ 0 ). (10) Usng agan the same homogeneous boundary condtons w(0)=w(1)=0, now we get an nhomogeneous lnear system for each k. Let us ntroduce the norm of the soluton as F (k)= 1 w(x n ) 2, F d (k)=f(k)/f(k 0 ), (11) where F d (k) s the dmensonless value, k 0 s a reference wave number and the ponts x n are randomly dstrbuted n [0,1]. In all the calculatons presented n ths secton we use = 5. Ths functon characterzes the value of the response of the system to the outer exctaton. ote that the rght hand sde f correspondng to (10) equals to zero dentcally nsde [0,1] and BVP (8) has a unque soluton w = 0forallk except k = k n - egenvalues when the soluton s not unque. In Fg. 1 the value F d as a functon of the wave number k s shown. The graph contans large sharp peaks at the postons of egenvalues. Generally speakng, ths resonance curve can be used to determne the egenvalues n the same way as det[a (k)] n the technque descrbed above. However, the graph F d (k) s a non smooth one, as t s shown n the lower part of Fg. 1 wth more detals. Ths can be explaned by the followng reasons. Problem (8), (9) wth w p gven n (10) has the exact soluton q 1 = 0, q 2 = e k(ξ 0 ξ 2 ) and so the total soluton w(x)=0, for x [0,1]. So, here we have F (k) whchsequaltozerowth machne precson accuracy when k s far from egenvalues; F (k) grows consderably n a neghbourhood of the egenvalues when the lnear system becomes almost degenerated. And a smoothng procedure s needed to get an approprate curve whch s convenent for applyng an optmzaton procedure. The followng two smoothng procedures are used n the paper. 2.1 smoothng by a dsspatve term The frst procedure conssts of ntroducng an addtonal dsspatve term n the governng equaton. And nstead of (8) we consder the problem: d 2 w dx 2 +( k 2 +εk ) w = f, w(0)=w(1)=0. (12) lnf d 15 10 5 F d 10 8 6 4 2 5 10 15 20 25 30 3 4 5 6 7 k Fgure 1 : Resonance curve n 1D egenproblem. lnf d 15 10 5 15 10 5 F d 5 10 15 20 25 30 k 3 4 5 6 7 k Fgure 2 : Resonance curve n 1D egenproblem. Smoothng by a frcton term. k
180 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 Here ε s a small parameter. Ths means that the ntal wave equaton (5) s changed by the equaton ttu 2 = 2 xxu ε t u whch descrbes vbraton of a homogeneous strng wth a frcton [Morse and Feshbach (1953)]. The fundamental soluton s: Ψ(x,ξ,k,ε)= 1 exp (χ x ξ ), 2χ χ = k 2 +εk. (13) ow the system w(0) =0, w(1)=0 wth w p gven n (10) has a unque non zero soluton for all real k. The resonance curve correspondng to ε = 10 6 s shown n Fg. 2. ow ths s a smooth curve wth separated maxmums at the postons of egenvalues. To fnd the egenvalues we use the followng algorthm through the paper. Let us look for the egenvalues on the nterval [a,b] Then: (A) step 0: Choose h > 0; f F (a) > F (a +h) goto step 5; step 1: x 1 = a; F1 = F (x 1 ); step 2: x 2 = x 1 +h; F2 = F (x 2 ); f x 2 > b stop; step 3: f F2 > F1 then [F1 = F2;x 1 = x 2 ]; step 4: goto step 2; fnd the maxmum pont x m of F (x) on [x 2 2h,x 2 ]; step 5: x 1 = a; F1 = F (x 1 ); step 6: x 2 = x 1 +h; F2 = F (x 2 ); f x 2 > b stop; step 7: f F2 < F1 then [F1 = F2;x 1 = x 2 ; goto step 6]; else goto step 2. ote that any unvarate optmzaton procedure can be used at step 4. In partcular, we appled Brent s method based on a combnaton of parabolc nterpolaton and bsecton of the functon near to the extremum(see [Press, Teukolsky, Vetterlng, and Flannery (2002)], Ch. 10 and [Brent (1973)], Ch. 5 ). The step s taken h = 0.01 through the paper. The data placed n Tab. 1 are obtaned by applyng ths technque wth ε = 0.1, 10 3, 10 6. Other parameters are: ξ 1 = 0.5, ξ 2 = 1.5, ξ 0 = 5. Here we place the relatve errors e r = k k (ex) /k (ex) (14) n the calculaton of the frst fve egenvalues. Table 1 : One dmensonal egenproblem. The relatve errors n calculatons of the egenvalues. Smoothng by the frcton term. k (ex) ε = 0.1 ε = 10 3 ε = 10 6 π 1.3 10 4 1.3 10 8 1.7 10 12 2π 3.2 10 5 3.1 10 9 1.6 10 12 3π 1.4 10 5 1.4 10 9 1.5 10 12 4π 7.9 10 6 7.9 10 10 9.7 10 13 5π 5.1 10 6 5.0 10 10 9.0 10 13 2.2 smoothng by a shft of wave numbers The second smoothng technque s as followng. Let us ntroduce the constant shft k between the exctng source and the studed mode,.e., nstead of (10), we take the partcular soluton n the form: w p = Ψ(x,ξ 0,k + k) 1 = 2(k + k) exp ((k + k) x ξ 0 ). (15) ow the lnear system w(0)=w(1)=0 has non zero solutons for all k except the egenvalues k n when the system becomes degenerate. However, due to teratve procedure of soluton and roundng errors we never solve the system wth the exact k n. And we observe degeneraton of the system as a consderable growth of the soluton n a neghbourhood of the egenvalues. The resonance curve correspondng to k = 1 s shown n Fg. 3. Some results of the calculatons we got usng the second smoothng technque are presented n Tab. 2. The values ξ 1, ξ 2, ξ 0 are the same as above. Below we wll name these procedures as ε procedure and k procedure. Table 2 : One dmensonal egenproblem. The relatve errors n calculaton of the egenvalues. Smoothng by shft of the wave numbers. k (ex) k = 0.1 k = 1 k = 10 π 1.4 10 11 9.1 10 12 7.8 10 12 2π 5.8 10 13 3.5 10 12 5.5 10 12 3π 6.4 10 12 1.3 10 12 3.5 10 12 4π 3.3 10 13 2.8 10 12 2.3 10 12 5π 5.3 10 12 3.5 10 12 5.9 10 13
The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators 181 lnf d 30 20 10 15 10 5 F d 5 10 15 20 25 30 k 3 4 5 6 7 k Fgure 3 : Resonance curve n 1D egenproblem. Smoothng by a shft between the wave numbers. 3 Helmholtz egenproblem Applyng the MFS to problem (4) we look for an approxmaton soluton n the form of a lnear combnaton: w(x q)=w p (x)+ q n Φ n (x), (16) where w p s the partcular soluton correspondng to f, and the tral functons Φ n (x)=h (1) 0 (k x ζ n ) (17) satsfy the homogeneous PDE. Ths s the so-called Kupradze bass [Kupradze (1967)]. The sngular ponts ζ n are located outsde the soluton doman. The free parameters q n should be chosen to satsfy the boundary condton B 1 [w(x q)] = 0, x Ω. In partcular the unknowns q n are taken as a soluton of the mnmzaton problem: mn q c =1 { B 1 [w p (x )] + q n B 1 [Φ n (x )]} 2 (18) Here the ponts x, = 1,..., c are dstrbuted unformly on the boundary. We take c approxmately twce as large as the number of free parameters. The problem s solved by the standard least squares procedure. ote that we get (18) as a result of dscretzaton of the ntegral condton: Z mn {B 1 [w(x q)]} 2 ds w Ω More detals of ths technque can be found, e.g., n [Farweather and Karageorghs (1998); Golberg and Chen (1998)]. As a partcular soluton correspondng to the exctng source we take the same fundamental soluton w p (x)=φ ex (x,ζ ex,k) H (1) 0 (k x ζ ex ) (19) wth ζ ex placed outsde the soluton doman. When dealng wth problems n multply connected domans, the same tral functons can be used. And the source pontsshould be placed also nsde each hole. As an alternatve approach one can use the specal tral functons assocated wth each hole: Φ s,1 (x)=h (1) 0 (kr s), Φ s,2n+1 (x)=h n (1) (kr s )cosnθ s, Φ s,2n (x)=h n (1) (kr s )snnθ s. (20) Here r s = x x s,θ s s the local polar coordnate system wth the orgn at x s. Ths s so-called Vekua bass [Vekua (1957); Hafner (1990)], or multpole expanson. It s proven that every regular soluton of the 2D Helmholtz equaton n a doman wth holes can be approxmated wth any desred accuracy by lnear combnatons of such functons f the orgn x s of a multpole s nsde every hole. In ths case nstead of (16) we use: w(x q,p s )=w p (x)+ + S M s=1 m=1 q n Φ n (x) p s,m Ψ s,m (x), (21) where S s the number of holes and M s the number of terms n each multpole expanson. When the ε smoothng procedure s appled, then nstead of (4) we consder the problem: 2 w h + ( k 2 +εk ) w h = 0, x Ω, B 1 [w h (x)] = B 1 [w p (x)], x Ω. (22)
182 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 wth some small ε > 0. ote that ths problem has a unque nonzero soluton for all real k. Then the tral functons (17) should be also modfed: Φ n (x)=h (1) 0 (χ x ζ n ), χ(k,ε)= k 2 +εk. (23) Applyng the k procedure we modfy the partcular soluton whch should be taken n the form: w p (x)=φ ex (x, k) H (1) 0 ( k x ζ ex ), k = k + k. (24) 3.1 numercal examples Here the results of the numercal experments are gven to llustrate the method presented. In all the cases consdered below the resonance curve F (k) s computed usng t testng ponts x t,l Ω: F (k)= 1/ t t l=1 w(x t,l) 2. In all the calculatons we use 15 testng ponts dstrbuted nsde Ω wth the help of RUF generator of pseudorandom numbers from the Mcrosoft IMSL Lbrary. To get the egenvalues we look for the maxma of F (k) usng the Brent s procedure mentoned. Example 1) A crcular doman wth the radus r = 1 subjected to Drchlet or eumann boundary condton s consdered. The exctng source s placed at the poston ζ ex =(5,5); the sngular ponts ζ n of the fundamental solutons (17) are located on the crcle wth the radus R = 2. The results shown n Tab. 3 correspond to ε = 10 6. Here we place the relatve errors (14) n the calculaton of the frst 5 egenvalues. The lne n a cell ndcates that the soluton process faled wth these parameters. The exact egenvalues k (ex) are the roots of the equaton J n (k)=0 (Drchlet) or J n (k)=0 (eumann). Example 2) The role of the parameter ε s shown n Tab. 4. We solve the same problem as above wth Drchlet condton. Here we fx the number of free parameters = 25 and vary the parameter ε. The parameter ε coarsens the system. For a large ε we can calculate all the egenvalues k, = 1,...,10 but the precson s not very hgh. When ε decreases, the precson n determnng of k ncreases but t fals for large. The fgures Fg. 4, Fg. 5, Fg. 6 correspond to the data placed n Tab. 4. For ε = 10 2 the resonance peaks are spread because the frcton. When ε decreases the peaks become more sharp and narrow. Besdes for ε = 10 8 the peaks correspondng to k, > 2 are placed on the rsng Table 3 : Crcular doman wth the radus r = 1. The relatve errors n calculatons of the egenvalues. ε procedure; ε = 10 6. Drchlet condton = 15 = 20 = 25 1 8 10 11 8 10 12 7 10 12 2 2 10 3 5 10 11 2 10 11 3 3 10 9 1 10 9 1 10 9 4 2 10 3 4 10 11 1 10 11 5 6 10 7 2 10 3 1 10 9 eumann condton 1 2 10 9 2 10 9 2 10 9 2 4 10 9 2 10 9 2 10 9 3 9 10 12 1 10 11 6 10 12 4 7 10 8 9 10 10 8 10 10 5 2 10 6 6 10 10 3 10 10 Table 4 : Crcular doman wth the radus r = 1. Drchlet condton. The relatve errors n calculatons of the egenvalues. ε procedure wth varyng ε, = 25. ε = 10 2 ε = 10 4 ε = 10 6 ε = 10 8 1 6.4 10 6 6.0 10 10 7.3 10 12 4.9 10 11 2 2.4 10 6 1.9 10 10 2.0 10 11 4.3 10 11 3 3.2 10 6 1.4 10 9 1.0 10 9 4 9.0 10 7 1.6 10 10 1.3 10 11 5 1.1 10 6 1.6 10 9 1.4 10 9 6 6.5 10 7 1.5 10 10 7 4.9 10 7 4.8 10 10 8 2.7 10 6 1.1 10 9 9 4.9 10 7 5.9 10 9 10 5.2 10 6 sharply part of the resonance curve. As a result the algorthm (A) jumps over the egenvalues and one should decrease the step parameter h to capture the maxma. As t s shown n Tab. 4, for ε = 10 8 the algorthm fnds k 1 and k 2 wth h = 0.01. When h s reduced to 0.001 then the algorthm also gves the egenvalues k 3 and k 4.Toget k, = 1,...,10 one should take h = 0.0001. However, the algorthm becomes hghly expansve n the CPU tme. Example 3) ext, we consder the case when Ω s the unt square wth the same Drchlet or eumann boundary condton.ths problem has an analytcal soluton:
The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators 183 500 F d 400 300 200 100 3 4 5 6 Fgure 4 : Crcular membrane wth the radus 1. Drchlet condtons. ε - procedure wth ε = 10 2. 500 F d 400 300 200 100 3 4 5 6 Fgure 5 : Crcular membrane wth the radus 1. Drchlet condtons. ε - procedure wth ε = 10 4. 500 F d 400 300 200 100 3 4 5 6 Fgure 6 : Crcular membrane wth the radus 1. Drchlet condtons. ε - procedure wth ε = 10 8. k (ex) = π 2 + j 2,, j = 1,2,..(Drchlet condton) and, j = 0,1,2,..(eumann condton). In Tab. 5, we show k k k the results of calculaton of the frst 5 egenvalues wth ε = 10 6. The placement of the sngular ponts ζ n and the exctng source are the same as above. Example 4) For the next example, we consder an annular case of the doubly connected doman between the two crcles: Ω = { (x 1,x 2 ) r1 2 x2 1 +x2 2 2} r2 The nner and outer rad of an annular doman are r 1 = 0.5 and r 2 = 2, respectvely. We take Drchlet condton on the outer boundary and eumann on the nner one. The sngular ponts are dstrbuted at the crcles wth the rad a = 5(outsde the doman) and b = 0.3 (nsde the hole). The number of the sngular ponts on each auxlary contour s equal to. The exctng source s placed at ζ ex =(10,10). In Tab. 6 we present the relatve errors (14) n calculaton of the frst 5 egenvalues of the problem descrbed wth ε = 10 5. The values k (ex) are obtaned numercally as the roots of the equaton: J n (r 1 k)y n (r 2 k) J n (r 2 k)y n (r 1 k)=0. Table 5 : Square wth the sde a = 1. The relatve errors n calculatons of the egenvalues. ε procedure; ε = 10 6. Drchlet condton = 15 = 20 = 25 1 1 10 6 3 10 8 1 10 9 2 1 10 5 9 10 8 1 10 8 3 8 10 5 3 10 8 8 10 9 4 3 10 4 1 10 6 3 10 9 5 3 10 3 4 10 5 6 10 7 eumann condton 1 4 10 7 5 10 8 8 10 12 2 1 10 6 3 10 8 3 10 9 3 4 10 5 1 10 7 3 10 10 4 1 10 4 6 10 6 5 10 9 5 5 10 4 2 10 5 6 10 7 Table 6 : Annular doman. The relatve errors n calculatons of the egenvalues. ε procedure; ε = 10 5. k (ex) = 15 = 20 = 25 1 1.3339427880 5 10 11 2 10 11 2 10 11 2 1.7388632616 6 10 8 7 10 12 5 10 12 3 2.4753931967 7 10 11 8 10 12 4 3.1645013237 7 10 8 5 10 11 5 3.2899912986 7 10 11
184 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 Example 5) In ths example, doubly connected regon wth the nner regon of vanshng maxmal dmenson s consdered. The geometry of the problem s the same as n Example 3. However, here we consder the case of very small nner holes. In partcular, we take r 1 = 10 1,10 2,10 3 wth the same fxed r 2 = 2. ow, the Kupradze type bass functons (17) are unft to approxmate the soluton n a neghbourhood of the hole. Here we use a combned bass whch ncludes the tral functons (17) wth the sngular ponts placed on an auxlary crcular contour outsde the soluton doman and a multpole expanson wth the orgn at the center of the hole. Thus, we look for an approxmate soluton n the form: w(x q,p)=w p (x)+ q n Φ n (x)+ M m=1 p m Ψ m (x). The data presented n Tab. 7, Tab. 8, Tab. 9 correspond to the number of sources on the outer auxlary crcular contour = 50. The number of terms n multpole expanson M vares from M = 11(r 1 = 10 1 ) to M = 5(r 1 = 10 3 ). The exctng source s placed at the poston ζ ex =(10,10). We use the k procedure wth the shft k = 1. We would lke to draw the readers attenton to the fact that the method presented can separate very close egenvalues: k (ex) 4 = 3.1900833197 and k (ex) 5 = 3.2126996563(see data correspondng to r 1 = 10 1 ). Here the step n the algorthm (A) s taken h = 0.001. The detaled dscusson of Vecua bass for Helmholtz equaton can be found n [Hafner (1990)]. 4 Egenproblems wth bharmonc operator Accordng to the technque proposed, nstead of (2) let s consder BVP 4 w k 4 w = f, x Ω R 2, (25) w = 0, B 2 [w]=0, x Ω. In applcaton to ths problem, the MFS technque s smlar to the one consdered n the prevous secton. The tral functons now are of the two types: the fundamental solutons of the Helmholtz operator 2 +k 2 : Φ (1) n (x)=h (1) 0 (k x ζ n ) (26) consdered above and the fundamental solutons of the modfed Helmholtz operator 2 k 2 : Φ (2) n (x)=h (1) 0 (k x ζ n )= 2 π K 0 (k x ζ n ), (27) Table 7 : Crcle wth a small hole. Drchlet boundary condton. The outer radus: r 2 = 2. The relatve errors n calculaton of the frst ten egenvalues. k procedure wth k = 1. r 1 = 0.1, = 50, M = 11 k (ex) e r 1 1.5322036536 1.9 10 8 2 1.9301625755 5.8 10 9 3 2.5680354360 1.6 10 9 4 3.1900833197 1.3 10 11 5 3.2126996563 7.4 10 9 6 3.5522743165 3.7 10 10 7 3.7941712382 1.2 10 11 8 4.2101115868 9.0 10 12 9 4.3857419081 4.4 10 12 10 4.8805392651 1.0 10 11 Table 8 : Crcle wth a small hole. Drchlet boundary condton. The outer radus: r 2 = 2. The relatve errors n calculaton of the frst ten egenvalues. k procedure wth k = 1. r 1 = 0.01, = 50, M = 7 k (ex) e r 1 1.3709447159 2.5 10 8 2 1.9160005377 5.4 10 9 3 2.5678112121 1.6 10 9 4 2.9632630840 5.3 10 9 5 3.1900809955 2.9 10 12 6 3.5082790790 2.3 10 12 7 3.7941712738 1.0 10 9 8 4.2086222910 7.6 10 12 9 4.3857419733 1.1 10 11 10 4.5543927267 1.3 10 9 where H (1) 0 s the Hankel functon and K 0 s the modfed Bessel functon of the second knd and of order zero. So, an approxmate soluton s sought n the form of the lnear combnaton: w(x q 1,q 2 )=w p (x)+ q 1,n Φ (1) n (x)+ q 2,n Φ (2) n (x). (28) where w p (x) s a partcular soluton correspondng to the external source. We take t n the same form as n the
The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators 185 Table 9 : Crcle wth a small ole. Drchlet boundary condton. The outer radus: r 2 = 2. The relatve errors n calculaton of the frst ten egenvalues. k procedure wth k = 1. r 1 = 0.001, = 50, M = 5 k (ex) e r 1 1.3148533741 2.0 10 8 2 1.9158544900 5.4 10 9 3 2.5678111892 1.5 10 9 4 2.8883437835 2.8 10 9 5 3.1900809955 1.1 10 10 6 3.5077982552 3.0 10 11 7 3.7941712738 1.2 10 11 8 4.2086221329 5.9 10 12 9 4.3857419733 1.2 10 12 10 4.4650868082 3.6 10 10 Table 10 : A crcular plate wth the radus: r = 1;the relatve errors n calculaton of the frst eght egenvalues. k procedure, k = 0.1. w = w/ n = 0 = 20 = 25 = 30 1 3 10 9 3 10 9 3 10 9 2 7 10 7 3 10 9 2 10 9 3 1 10 5 8 10 8 6 10 10 4 6 10 10 6 10 10 6 10 10 5 7 10 5 9 10 7 1 10 8 6 1.2 10 5 1 10 7 2 10 9 7 9.4 10 4 8 10 6 8 10 8 8 1.3 10 4 2 10 6 2 10 8 w = 2 w/ n 2 = 0 1 2 10 9 3 10 9 3 10 9 2 5 10 5 3 10 6 1 10 7 3 3 10 4 2 10 5 8 10 7 4 4 10 9 5 10 9 4 10 9 5 2 10 3 9 10 5 4 10 6 6 3 10 4 1 10 5 6 10 7 7 9 10 3 3 10 4 2 10 5 8 9 10 2 6 10 5 3 10 6 prevous secton,.e. (19). The free parameters are determned from the boundary condtons. We apply the same ε and k smoothng procedure. When the ε procedure s appled the governng equaton should be replaced by the followng one: 4 w(x) ( k 4 +εk 2) w(x)= f (29) and so the arguments of the tral functons Φ (1) n (x), (x) should be modfed. Applyng the k procedure Φ (2) n we modfy the external source and take t n the form (24). 4.1 numercal examples Example 6) A crcular plate wth the radus r = 1 subjected to the boundary condtons: a) w = w/ n = 0 (clamped boundary) and b) w = 2 w/ n 2 = 0 s consdered. The exctng source s placed at the poston ζ ex =(5,5); the sngular ponts ζ n of the fundamental solutons (26), (27) are located on the crcle wth the radus R = 2. Remark that now the number of free parameters s 2. The data presented n Tab. 10 are obtaned usng k procedure wth k = 0.1. Here we place the relatve errors (14). The exact egenvalues k (ex) are the roots of the equaton J n (k)i n (k) J n (k)i n (k)=0(condtons a)) or J n (k)i n (k) J n (k)i n (k)=0 (condtons b)). Example 7) ext, we consder a square plate wth the sde a = 1 subjected to the boundary condtons w = 2 w/ n 2 = 0. Ths problem has an analytcal soluton: k (ex) = π 2 + j 2,, j = 1,2,.. The results placed n Tab. 11 are obtaned usng k procedure wth k = 0.1. Example 8) A rectangular 1.2 0.9 plate subjected to the boundary condtons w = w/ n = 0 (clamped boundary) s consdered. The results placed n Tab. 12 are obtaned usng k procedure wth k = 0.1. In ths case, the analytc soluton s not avalable. The results obtaned n [Chen, Chen, Chen, Lee, and Yeh (2004)] and [Kang and Lee (2001)] are used for comparson. These data are placed n the last two columns of the table. ote that usng ε procedure wth ε = 0.01 and = 56, we get the followng egenvalues: k 1 = 5.95263, k 2 = 7.70983, k 3 = 9.12854, k 4 = 10.27133, k 5 = 11.96763, k 6 = 12.49617 5 Concludng remarks In ths paper, a new meshfree method for egenproblems wth Laplace and bharmonc operators s proposed. Ths s a mathematcal model of physcal measurements, when a mechancal or acoustc system s excted by an external source and resonance frequences can be determned usng the growth of ampltude of oscllatons near these
186 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 Table 11 : A square plate. The relatve errors n calculaton of the frst sx egenvalues. k procedure, k = 0.1. = 20 = 25 = 30 = 40 1 3.3 10 6 2.0 10 7 1.8 10 8 1.7 10 8 2 6.9 10 4 1.5 10 5 3.2 10 8 1.7 10 8 3 7.9 10 5 3.7 10 6 1.7 10 8 4 9.2 10 5 5.5 10 7 1.7 10 8 5 4.5 10 2 1.5 10 5 2.2 10 8 6 1.4 10 3 1.2 10 7 Table 12 : A rectangular plate 1.2 0.9 wth clamped boundary. = 35 = 42 = 49 I II 1 5.9515 5.9529 5.9527 5.952 5.952 2 7.7125 7.7116 7.7104 7.703 7.703 3 9.1333 9.1319 9.1316 9.129 9.131 4 9.9466 9.9510 9.9493 9.947 9.955 5 10.2692 10.2717 10.2742 10.266 10.27 6 11.9501 11.9552 11.9565 11.95 11.95 7 12.3849 12.3719 12.3710 Table 13 : The BKM soluton. Crcular doman wth Drchlet condtons. The relatve errors n calculatons of the egenvalues. k procedure; k = 0.1. = 10 = 14 = 20 = 30 1 2 10 4 2 10 6 4 10 9 7 10 9 2 3 10 4 4 10 7 1 10 10 1 10 8 3 9 10 5 2 10 8 1 10 8 4 4 10 7 4 10 9 5 1 10 6 8 10 9 frequences. The method shows a hgh precson n smply and multply connected domans. The dea can be extended qute smply to the 3D case. The method presented s based on the MFS soluton of the problem. However, t can be combned wth other boundary technques. The BKM mentoned n Secton 1 seems to be perspectve n ths connecton. For example, f the BKM s appled to Helmholtz equaton, the approxmaton soluton s looked for n the form: w(x q)=w p (x)+ q n J 0 (k x ζ n ) cf. (16). Here the source ponts ζ n can be placed nsde the soluton doman. To test BKM n the framework of the method presented we solve the same problem as the one descrbed n Example 1 wth Drchlet condton. The half of the source ponts ζ n, n = 1,...,1/2 are placed unformly on the boundary Ω. The rest source ponts ζ n, n = 1/2 + 1,..., are dstrbuted nsde Ω wth the help of the generator of pseudorandom numbers. The data presented n Tab. 13 are obtaned usng k procedure wth k = 0.1. The parameters of the exctng source are the same as above n Example 1. It should be noted that the BKM and the MFS, as well as the all methods of the Trefftz type n general, have a narrow feld of applcaton. It s restrcted by the cases when there exsts a representatve set of known exact solutons of PDEs under consderaton,.e. by the problems posed by lnear PDEs wth constant coeffcents. See, however, [Reutsky (2002)], where a Trefftz type technque s developed for PDEs wth varyng coeffcents. Besdes the Trefftz type technques produce the systems of equatons wth unsymmetrc fully populated matrces. As a result, the MFS s hghly ll condtoned. In some cases one can overcomes ths drawback by the use of matrces of the specal block crculant structure and an effcent matrx decomposton technque [Tsangars, Smyrls, and Karageorghs (2004)]. However, takng n mnd further applcatons of the method presented n the paper to egenproblems wth PDEs of general type n rregular domans, one should combne t wth meshless methods based on the local approxmaton of the soluton lke the Meshless Local Petrov-Galerkn Method [Atlur (2004), Han and Atlur (2003), Han and Atlur (2004)]. The comparson between global and local approxmaton, e.g. BEM and FEM, and they combnaton see n [Grannell and Atlur (1978)]. Comparng the method wth the technque based on computatons of the determnant of the system, the followng crcumstances should be taken nto account. Snce the MFS s hghly ll condtoned, the determnant s very small. Indeed, let us consder agan the same egenvalue problem whch s descrbed n Example 1,.e. Helmholtz equaton n the crcle wth the radus 1 and Drchlet boundary condton. We take the number of the sources equal to the number of the collocaton ponts on the boundary. Thus, we get a square matrx of the problem
The method of fundamental solutons for egenproblems wth Laplace and bharmonc operators 187 Table 14 : Crcular doman wth Drchlet condtons. The number of the source ponts = 30; ε procedure. ε = 10 1 ε = 10 4 ε = 10 6 e r F(k ) e r F(k ) e r F(k ) 1 4 10 4 0.701 4 10 10 0.701 5 10 12 0.701 2 2 10 4 0.652 1 10 10 0.654 6 10 11 0.654 3 9 10 5 0.509 9 10 10 0.516 1 10 9 0.516 A(k,) and can calculate the determnant deta(k,). Placng the sources on the crcle wth radus 2 and takng k = 1weget: deta(1,20) = 3 10 47, deta(1,30) = 4 10 117, deta(1,40) = 3 10 217. The wave number k = 1 s not the egenvalue of the problem. Ths s the background value between extremums and one looks for the mnma of deta(k,) on such background. So, usng ths technque one operates wth values of the order 10 50 10 500, see [Alves and Antunes (2005); Chen, Chen, and Lee (2005)] for more detaled nformaton. At the same tme let us calculate the norm functon F(k,) whch s used to obtan the egenvalues n the method presented. We get for ε = 0.0001: F(1,20)= 2.13 10 5, F(1,30)=2.13 10 5, F(1,40)=2.13 10 5. We present the values of the norm functon F (k) when k s close to egenvalue n Tab. 14. Here the number of the sources s fxed = 30 and the smoothng parameter ε s vared. e r s the relatve error n determnng of the approxmated egenvalue k and F (k ) denotes the value of the norm functon at ths approxmated egenvalue. So, n the framework of the method presented we always deal wth the values whch can be handled on PC wth a sngle precson. The method s easy to program and not expensve n the CPU tme. The all calculatons presented n the paper were performed usng 366 MHz PC. References Alves, C. J. S.; Antunes, P. R. S. (2005): The Method of Fundamental Solutons appled to the calculaton of egenfrequences and egenmodes of 2D smply connected shapes. CMC: Computers, Materals & Contnua, vol. 3, pp. 00 00. Atlur, S. (2004): The Meshless Method (MLPG) for Doman & Boundary Dscretzatons. Tech Scence Press, Forsyth. Brent, R. P. (1973): Algorthms for mnmzaton wthout dervatves. Prentce-Hall, Englewood Clffs, J. Chen, J. T.; Chen, I. L.; Chen, K. H.; Lee, Y. T.; Yeh, Y. (2004): A meshless method for free vbraton analyss of crcular and rectangular clamped plates usng radal bass functon. 535 545. Eng. Anal. Bound. Elem., vol. 28, pp. Chen,J.T.;Chen,I.L.;Lee,Y.T.(2005): Egensolutons of multply connected membranes usng the method of fundamental solutons. Eng. Anal. Bound. Elem.,vol. 29, pp. 166 174. Chen, J. T.; Ln, J. H.; Kuo, S. R.; Chyuan, S. W. (2001): Boundary element analyss for the Helmholtz egenvalues problems wth a multply connected doman. Proc. R. Soc. Lond. A, vol. 457, pp. 2521 2546. Chen, J. T.; Lu, L. W.; Hong, H. K. (2003): Spurous and true egensolutons of Helmholtz BIEs and BEMs for a multply connected problem. Proc. R. Soc. Lond. A, vol. 459, pp. 1897 1924. Chen, W. (2005): Symmetrc boundary knot method. Eng. Anal. Bound. Elem., vol. 26, pp. 489 494. Chen, W.; Tanaka, M. (2002): A meshless, exponental convergence, ntegraton-free, and boundary-only RBF technque. Computers and Mathematcs wth Applcatons, vol. 43, pp. 379 391. Courant, R. (1943): Varatonal methods for the soluton of problems of equlbrum and vbratons. Bull. Amer. Math. Soc., vol. 43, pp. 1 23. Courant, R.; Hlbert, D. (1953): Methods of Mathematcal Physcs. Wley Interscence, ew York. Farweather, G.; Karageorghs, A. (1998): The method of fundamental solutons for ellptc boundary
188 Copyrght c 2005 Tech Scence Press CMC, vol.2, no.3, pp.177-188, 2005 value problems. 69 95. Advances n Comp. Math., vol. 9, pp. Golberg, M. A.; Chen, C. S. (1997): Dscrete Projecton Methods for Integral Equatons. Computatonal Mechancs Publcatons, Southampton. Golberg, M. A.; Chen, C. S. (1998): The Method of Fundamental Solutons for Potental, Helmholtz and Dffuson Problems, pp. 103 176. WIT Press, 1998. Golub, G. H.; Loan, C. F. V. (1996): Matrx Computatons. Johns Hopkns Unversty Press, Baltmore. Grannell, J. J.; Atlur, S.. (1978): Boundary Methods (BEM) and Combnaton of BEM-FEM. GIT-ESM- SA-78-16. Reutsky,S.Y.(2002): A boundary method of Trefftz type wth approxmate tral functons. Engneerng Analyss wth Boundary Elements, vol. 26, pp. 341 353. Strang, G. (1976): Lnear Algebra and Its Applcatons. Academc Press. Strang, G.; Fx, G. (1973): An Analyss of the Fnte Element Method. Prentce-Hall. Tsangars, T.; Smyrls, Y. S.; Karageorghs, A. (2004): A Matrx Decomposton MFS Algorthm for Bharmonc Problems n Annular Domans. CMC: Computers, Materals, & Contnua, vol. 1, pp. 245 258. Vekua, I. (1957): ew Methods for Solvng Ellptc Equatons. orth-holland, Amsterdam. Hafner, C. (1990): The Generalzed Multpole Technque for Computatonal Electromagnetcs. Artech House Books, Boston. Han, Z. D.; Atlur, S.. (2003): Truly Meshless Local Petrov-Galerkn (MLPG) Solutons of Tracton & Dsplacement BIEs. CMES: Computer Modelng n Engneerng & Scences, vol. 4, pp. 665 678. Han, Z. D.; Atlur, S.. (2004): A Meshless Local Petrov-Galerkn (MLPG) Approach for 3-Dmensonal Elasto-dynamcs. CMC: Computers, Materals, and Contnua, vol. 1, pp. 129 140. Kang,S.W.;Lee,J.M.(2001): Free vbraton analyss of arbtrary shaped plates wth clamped edges usng wave-type functons. J. Sound. Vbr., vol. 242, pp. 9 26. Karageorghs, A. (2001): The method of fundamental solutons for the calculaton of the egenvalues of the Helmholtz equaton. Appled Math. Letters., vol. 14, pp. 837 842. Kupradze, V. D. (1967): On approxmate soluton of problems n mathematcal physcs. Russan Math. Surveys, vol. 22, pp. 58 108. Methods of The- Morse, P. M.; Feshbach, H. (1953): oretcal Physcs. McGraw Hll. Press, W. H.; Teukolsky, S. A.; Vetterlng, W. T.; Flannery,B.P.(2002): umercal Recpescpes n C++. Cambrdge Unversty Press.