Application of the multiquadric method for numerical solution of elliptic partial differential equations

ii Application of the Multiquadric Method for Numerical Solution of Elliptic Partial Differential Equations Maithili S h a r a n

Centre for Atmospheric Sciences Indian Institute of Technology Hauz Khas, New Delhi-110016, India and E. J. K a n s a

Environmental Science Program Directorate Mail Stop L-200 Lawrence Livermore National Laboratory L ivermore, California 94551-9900 and Suman Gupta

Government Girls Senior Secondary School I Madangir, New Delhi, India Transmitted by Michael Golberg

ABSTRACT

We have used the multiquadric (MQ) approximation scheme for the solution of elliptic partial differential equations with Dirichlet a n d / o r Neumann boundary conditions. The scheme has the advantage of using the data points in arbitrary locations with an arbitrary ordering. Two-dimensional Laplace, Poisson, and biharmonic equations describing the various physical processes have been taken as the test examples. The agreement is found to be very good between the computed and exact solutions. The method also provides an excellent approximation with a curved boundary. © Elsevier Science Inc., 1997

APPLIED MA THEMA TICS AND COMPUTATION 84:275-302 (1997) © Elsevier Science Inc., 1997 655 Avenue of the Americas, New York, NY 10010

0096-3003/97/$17.00 PII S0096-3003(96)00109-9

276

M. SHARAN ET AL.

1. INTRODUCTION In general, the mathematical description of physical processes leads to partial differential equations (PDEs). In some cases, an exact solution of PDEs can be obtained using analytical tools such as reflection methods and superposition, separation of variables and integral transforms, etc. The approximate analytical solution of PDEs can be found by using perturbation methods [1-3], the method of successive approximation [4], and orthogonal functions [5]. The solution of equations that more closely model experimental and practical situations is possible by the methods of numerical analysis. The equations can be solved by using low order,finite difference, finite element and boundary element methods. When the boundary conditions and domains are irregular, the finite difference scheme becomes complex; it can be overcome to a greater extent by using the method of finite elements. The length scale, A, varies as h = If/V f[. In very steep gradient regions, low-order schemes require appropriately fine discretization that in higher dimensions can be costly to process. Also, such schemes require nearly uniform mesh spacings to minimize the introduction of spatial truncation errors to avoid wave distortions, etc. It is difficult to use finite difference methods on arbitrarily scattered data points. In 1968, Hardy [6] introduced the multiquadric (MQ) method for the construction of approximate two-dimensional surfaces of field data. The importance and usefulness of the MQ method is recognized by its successful applications in geology, geography, geophysics, hydrology, mapping, surveying, photogrammetry, remote sensing and signal processing, and digital terrain models. Hardy [6] has presented an excellent review summary on the various applications of the MQ method since 1968. Based on a comparison of different scattered data schemes, Franke [7] has concluded that MQ performs the best in accuracy and ease of implementation even against various finite element schemes. Kansa [8] has explored MQ scattered data schemes for interpolation of data points to approximate surfaces and estimating partial derivatives. Kansa [9] has pointed out that the interpolation scheme is intimately connected with the solution of PDEs, and he successfully applied this novel technique for the solution of PDEs in computational fluid dynamics. Several theoretical papers have been written showing the theoretical justification of the excellent performance of MQ in interpolation and approximation. First, Micchelli [1] showed that when polynomial terms are appended to the MQ expansion, the set of linear equations that arise are conditionally positive definite and always have a solution for distinct knots. Baxter [11] showed that as the MQ shape parameter goes to infinity, the

Solution o/ Elliptic PDEs

277

univariate MQ interpolant tends toward a sinc function, thus showing the relation between MQ and a class of well-known spectral methods. Madych and Nelson [12-14] showed that MQ interpolation converges exponentially as the density of the knots increases and is based upon a variational principle whose solution has minimal semi-norm error. Dubal et al. [15] and Dubal [16] have used the MQ scheme for the solution of a three-dimensional, highly nonlinear elliptic partial differential equation describing the initial state of the collision of two black holes. The papers of Dubal first pointed out the numerical observation that MQ, in the solution of PDEs, does indeed converge exponentially as the density of the knots increases. Although MQ has been proven to be exponentially convergent for interpolation and approximation, no theoretical studies have yet been developed for PDEs. The MQ method replaces the classical finite difference and finite element spatial discretization schemes by a exponentially convergent, grid-free scattered data approximation scheme. Finite difference and finite element methods use low-order polynomial basis functions. MQ in contrast is very high order [17]. They have discussed the advantages/disadvantages associated with finite difference, finite element, and MQ methods. In the present paper, we have applied the MQ scheme to solve the two-dimensional elliptic Laplace; Poisson, and biharmonic equations have been taken as the test examples. T h e method is also applied successfully to a problem with a curved boundary. 2.

APPLICATION OF T H E MQ METHOD Consider an elliptic PDE over the interior domain, ( t 2 / F ) , of the form. V2f = s(x)

(1)

where V 2 is the Laplacian operator, x is the spatial position s is a known function of x, and f is the unknown function of x. Equation (1) is subject to Dirichlet a n d / o r Neumann boundary conditions over the boundary, F f = fi(x) on r l n • Vf = f2(x) on F 2

(2.a)

(2.b)

where n is the outward unit normal, F 1 and F 2 are the boundaries of the domain, fl and f2 Are known functions of x.

278

M. SHARAN ET AL.

It is assumed that in an MQ approximation, any function can be expanded as a finite series of MQ basis functions [6, 9]. The expansion used in the Hardy expansion with an appended constant, i.e.,

N

/(~)

= ~ +

E

%~(x - xp

(3)

j=2

where N is the total number of data points under consideration, a's are the coefficients and the functions g are defined by, 3 , ( x - x j) = gj - gl,

(4)

where gj : g ( x -

Xj) : [ ( X - Xj)2 "~- R2] 1/2,

(5)

and (x - Xj) 2 is the square of the Euclidean distance in ~ ~ ( d > 0) and R~ > 0 is the shape parameter of the flh basis function. Carlson and Foley [18] have used a constant value of ( R j) 2 in (5) and obtained a heuristic recipe for its value for different functions. They observed that the root mean squared errors of many of the bivariate functions [7] were reduced several orders of magnitude when the optimal shape parameter was used. Kansa [8, 9] has found numerically that the following power law function increases accuracy up to five orders of magnitude for many monotonic functions. ( R j ) 2 : (Un)2[(Rm)2//(Rn)2]

[
(0

where R n and R m are input parameters. Later, Hagen and Kansa [19] studied the role of the parameter ( R j) 2 in the multiquadric function using a simple two-parameter optimization procedure developed by Marquadt [20]. Very recently, Golberg et al. [21] used the method of cross-validation to estimate the optimal shape parameter of elliptic PDE problems in two and three dimensions and observed exponential convergence. Spatial partial derivatives of the function in (3) are obtained by differentiating the spatial basis functions. Substitution of (3) and its partial derivatives in (1), together with the boundary conditions (2a and b), yields a system of linear algebraic equations which can be solved for the coefficients,

{aJ.

Solution of Elliptic PDEs

279

We illustrate the method for a two-dimensional problem. The derivatives at (x~, Yi) are given by

af/~((

N ~, y,) = E at[ ( x

--

xj)/gij

--

(x

--

xl)/gil],

j=2 N

j=2

at[( y - y~)/g,j- ( y -

~1)/g,1].

N 2

a2f/~z2( x~, y,) = E

2

j=2

--1//gi1{1 -- (X-- Xl)2//g21}], N

32f/ay2(x,,

Yi) = ~ aj[1/g,j(1-

( y-

yj) 2 /gi,}2

j=2

- 1 / g i l { 1 - ( Y-

2 Yl) 2 /gil}],

(7)

where

(8) Substituting the expansion of f, (3) and its derivatives (7), in the twodimensional version of (1), we find for each of the interior nodes, i, N

j=l

Pi3ai = bi

P/l = O~

-[(xi-

xl) 2 + ( y , +

y l ) 2 ] / ( g i l ) 3 for 2 _< j _< N,

(9)

and the source terms are represented as

b, = s( x,, y~).

(10)

280

M. SHARAN E T AL.

If a point i lies on the boundary, F, it m a y satisfy either a Diriehlet or a N e u m a n n b o u n d a r y condition. Case I: W i t h the Dirichlet condition f = fl, the m a t r i x elements for the i-th equation (9) will be P~I ~ 1, P i j = ]/ij = g i j -

for 2 < j < N,

gxj,

bl = f2( x~, yj).

(11)

Case II: N e u m a n n condition n • V f = f2. (a) If the point / lies on F where the n o r m is in the x direction, the m a t r i x elements for (9) are given by Pil = 0,

Po = ( x i -

x~)/g,j-

( x,-

xl)/g,1 ,

2
b, = f2( x~, y,).

(12)

(b) If the point i lies on F where the n o r m is in the y direction, the m a t r i x elements for (9) are given by Pil = O, P~j = ( y, -

y~)/g,j

- ( ~ - yx)/g,1

b, = f2( xi, Yi)-

2 _< j _< N

(13)

For all the knots (9) forms a system of linear algebraic equations of the form PA = B

(14)

where P = {po} is a square m a t r i x of order N, and A = (al, a2, . , . , a~) T and B = (bl, b 2 , . . . , bN) T are column vectors of length N, where A is the vector of unknown expansion coefficient, and B is the column vector of known source terms. T h e elements of P and B are defined by the relations (10-13).

Solution of Elliptic PDEs

281

Nonlinear PDEs are solved by replacing (3) everywhere with a nonlinear expansion of f in the governing set of PDEs, Dubal [15] solved for the expansion coefficients of the MQ expansion, (3), by Picard iteration, and Makroglou [22] solved for the expansion coefficients by a combination of Picard and Newton-Raphson schemes. Here we have solved the linear PDEs, (14), successfully using the Gauss elimination method with total pivoting. The mean square error ( E ) is computed from N

E= (l/N)

~'~ (fi - j~i)2

(15)

i=1

where f is the computed value of f using the MQ method and f is the known value of f at the point i. It may be known either analytically or computationally by means of some other methods such as finite differences, finite elements etc. We have used the double precision FORTRAN on ICL 3980 computing machine at IIT Delhi. 3.

EXAMPLES

In the following, we present the various examples for determining the solution of elliptic PDEs using the MQ method. 3.1. Laplace Equation in a Square Consider the Laplace equation, f ~ + fyy = 0

(16)

in a unit square domain 1~ (0 < x < 1, 0 < y < 1) with the conditions (i) f = 0 along the edges: x = 0; y = 1 and x = 1, (ii) f = T O along y = 0.

(17)

In (16), subscripts " x x " and " y y " denote the second partial derivatives with respect to x and y, respectively. TO is a constant (16), with the boundary conditions in (17), describes the steady-state temperature distribution in a square plate, one side of which is maintained at a constant temperature To, while the other three sides are maintained at 0. An analytical solution of (16) with the boundary conditions (17) is given by

282

M. SHARAN ET AL.

Carslaw and Jeager [23] oo

f ( x , y,) = (4To/qr) ~

[ 1 / ( 2 n + 1)]sin{(2n + 1)~rx}

Xsinh{(1 - y ) ( 2 n + 1)~-}cosech{(2n + 1)p}.

(18)

The problem described in (16), with the boundary conditions, (17), becomes a particular case of (1-2) with s( x, y) = 0 and f = To, with only Dirichlet conditions. We use the MQ method discussed in a previous section. The computation has been performed with TO = 100 for different number of knots (N). For 81 knots at a uniform spacing h x = A y = 0.125, we found the mean square error = 0.682 with a maximum error 13.4411. In the next case, we have removed the knots lying on the lines y = 0.375 and y -- 0.625, resulting in N -- 63. Further, the knots are reduced to 49 in the domain by deleting the points on the vertical lines x = 0.375 and y = 0.625. By removing the additional knots from the line y = 0.75, we have 42 points in the domain. The deletion of these points does not seriously degrade the accuracy of the MQ solution. Also, we have obtained the solution for scattered data points 30, 33, and 35, respectively, and for 25 (i.e., 5 x 5) uniform knots. Table 1 provides the information about the number of boundary points, number of interior points, the value of the input parameters R m and R~, mean square error and the maximum absolute error for each of the test run. Table 1 shows that mean square errors are not significantly different in all the cases. As an illustration, the computed results for N = 42 are given in Table 2. Referring to Fig. 1, there is a reasonably good agreement between the computed and

TABLE 1 VALUES OF INPUT PARAMETERS, MEAN SQUARE ERROR (E), AND MAXIMUM ABSOLUTE ERROR (MAE) SN

NT

Ns

N1

Rm

R,~

E

MAE

1 2 3 4 5 6 7 8

81 63 49 42 25 30 33 35

32 28 24 22 16 18 19 15

49 35 25 20 9 12 14 20

0.22 0.12 0.14 0.14 0.40 1.30 0.068 0.700

11 12 12 10 50 10 50 10

0.692 1.570 0.832 1.050 1.880 2.680 3.53 3.35

3.44 4.02 3.65 3.68 3.88 4.71 5.23 6.0

Solution of Elliptic PDEs

283

TABLE 2 COMPUTEDAND EXACTSOLUTIONOF THE LAPLACEEQUATION IN A SQUARE: N = 42, Rm = 0.14, R~ = 10 Index

X

Y

Computed

Exact

Error

*1 *2 *3 *4 *5 *6 *7 *8 9 10 11 12 13 "14 "15 16 17 18 19 20 "21 *22 23 24 25 26 27 *28 *29 30 31 32 33 34 *35 *36 *37 *38 *39 *40 "41 *42

0.000 0.125 0.250 0.500 0.750 0.875 1.000 0.000 0.125 0.250 0.500 0.750 0.875 1.000 0.000 0.125 0.250 0.500 0.750 0.875 1.000 0.000 0.125 0.250 0.500 0.750 0.875 1.000 0.000 0.125 0.250 0.500 0.750 0.875 1.000 0.000 0.125 0.250 0.500 0.750 0.875 1.000

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.250 0.250 0.250 0.250 0.250 0.250 0.250 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.875 0.875 0.875 0.875 0.875 0.875 0.875 1.000 1.000 1.000 1.000 1.000 1.000 1.000

0.0000 100.0000 100.0000 100.0000 100.0000 100.0000 0.00000 0.0000 51.9729 68.8883 74.9900 68.8720 51.9493 0.0000 0.0000 26.9965 44.2271 54.1377 44.2031 26.9641 0.0000 0.0000 11.0281 19.4364 26.1726 19.4055 10.9899 0.0000 0.0000 1.8004 3.1869 4.4374 3.1680 1.7750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 100.1664 99.9100 99.9363 99.9100 100.1664 0.0000 0.0000 48.2909 66.8953 75.4269 66.8953 48.2909 0.0000 0.0000 26.2587 43.2028 54.0529 43.2028 26.2587 0.0000 0.0000 10.0708 18.2028 25.0000 18.2028 10.0708 0.0000 0.0000 1.7091 3.1478 4.4316 3.1478 1.7091 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 -0.1664 0.0900 0.0637 0.0900 -0.1664 0.0000 0.0000 3.6820 1.9930 - 0.4369 1.9767 3.6583 0.0000 0.0000 0.7379 1.0243 0.0848 1.0002 0.7054 0.0000 0.0000 0.9573 1.2336 1.1726 1.2026 0.9191 0.0000 0.0000 0.0914 0.0391 0.0058 0.0202 0.0659 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

*Denotes boundary point.

284

M. S H A R A N E T AL. 100

811

60

: \ • \

50 40

'l *""1

'~ ' ,

30

XX

~,

,

"""'",.~. :

10

".

0 0.0

0.2

0.4

0.6

0.8

1.0

Y FIC. 1. Variation of temperature along the y axis for different values of x. The continuous lines from the bottom denote the MQ solution of x = 0.5, 0.25, and 0.125, respectively. The symbols denote the corresponding exact values. In this example, N = 42, R m = 0.14, R n = 10, and E = 1.06.

e x a c t values, c o n s i d e r i n g t h a t o v e r t h e u n i t s q u a r e t h e v a l u e s v a r y i n g 0 t o 100, e x c e p t at t w o p o i n t s w h e r e t h e m a x i m u m e r r o r is a b o u t I3.688[. T h e m a x i m u m e r r o r is r e l a t i v e l y h i g h e r a n d is i n t r o d u c e d a t t h e grid p o i n t n e a r t h e p o i n t of i n t e r s e c t i o n of t h e sides t e m p e r a t u r e s T o = 100 a n d T -- 0. T h u s , t h e large e r r o r is c a u s e d b y t h e d i s c o n t i n u i t y in t e m p e r a t u r e s T O -- 100 a n d T = 0. It r e m a i n s of t h e s a m e o r d e r e v e n w i t h t h e i n c r e a s e in t h e n u m b e r of k n o t s , see T a b l e 1.

3.2. Poisson Equation in a Elliptic Domain Consider the Poisson equation,

Lx + fyy =

-2

(19)

Solution of Elliptic PDEs

285

in an ellipse with major and minor semiaxes a and b. Equation (19) is subject to the condition that f vanishes on the surface of the ellipse.

(20)

x 2 / a 2 + y2/f12 = 1.

Since the solution of (19) with (20) is expected to be symmetric with respect to x and y, it is sufficient to solve (19) in the first quadrant of the ellipse and accordingly, we consider (19) with the following Neumann and Dirichlet boundary conditions

(i) L=oat x = 0 , (ii) f y = 0 a t

y=0,

(21)

(iii) f = 0 on the curved surface, F.

Equation (19) provides a function f from which the angle of twist of a cylindrical shaft of an elliptical cross-section under torsion can be calculated [24]. An analytical solution of (19) with (20) is

/ ( x , y) = - [

+

-

+

(22)

The solution of (19) with (21) can be found as a particular case of (1-2) with s(x, y ) = - 2 , f l = 0 a n d f2 = 0 - F ° r a = 10, b = 8 , and N = 28 (Figure 2), the results are presented in Table 3. Figure 3 is a representation of the knots used in the curved boundary problem. The function varies from 39 to 0 in domain with E = 2,92 × 10 -s. The maximum error is found to be 3.574 x 10 -4 which is 0.0025% of the exact functional value. The computed values are close to the exact values of the function. The MQ method provides accuracy up to five significant decimal places. Computations have also been done with different sets of data points and the results are found to be very good. For N = 87, the maximum error is found to be [6.68 × 10-51 which is 0.00027~. For a = 2, b = 1, the function varies from 0.8 to 0. Again, we obtain an accurate solution with N = 28, not only at the interior points, but also on the boundary points with E = 2.96 × 10 -1° and maximum error -- [4.8 X 10-~[. The MQ method provides an excellent approximate solution. The curved boundary is accounted for in this method in a natural way through the location of the boundary point, whereas, in the finite difference method, an interpolation technique is required for expressing boundary points in terms of knots.

286

M. SHARAN ET AL. 10 9 °•i

8

7 f 6

$ \

4

"'A

3 2 1

0 0.0

0.2

0.4

0.6

0.8

1.0

X

FIG. 2. Computed (continuous line) and exact (symbol) solution for the biharmonic equation versus x for y = 0 , y = 0 . 5 , and A y = 0 . 5 . Here, q = 2 0 , N = 2 5 , R ~ = 2 . 2 , R , = 0.15, and E = 3.89 • 10 -9.

3.3. Poisson Equation in a Rectangle with a Constant Right-hand Side Consider the Poisson equation

•f=x + fuy = Q

(23)

where Q is a constant. Equation (23) describes the steady state transport of oxygen in a slab of tissue with a constant oxygen consumption rate [25]. This also governs the steady state temperature distribution in a rectangular region with heat loss (production) at a constant rate Q ( - Q) [23]. If Q = q / D where D is the flexural rigidity of a thin plate subject to a normal load q per unit area, (34) arises by splitting a biharmonic equation (see example 3.5)) into two successive Poisson equations.

Solution of Elliptic PDEs

287 TABLE 3

COMPUTED AND EXACT SOLUTION OF THE POISSON EQUATION

IN AN ELLIPSE: a

=

10, b = 8, N = 28, R m = 55, R n = 4

Index

X

Y

Computed

*1 *2 *3 *4 *5 *6 *7 *8 *9 *10 11 12 13 14 15 16 17 18 19 20 21 *22 *23 *24 *25 *26 *27 *28

0.000 0.000 0.000 0.000 0.000 2.000 4.000 6.000 8.000 10.000 2.000 4.000 6.000 8.000 2.000 4.000 6.000 8.000 2.000 4.000 6.000 9.798 8.660 7.141 2.000 4.000 6.000 8.000

0.000 1.600 4.000 6.400 8.000 0.000 0.000 0.000 0.000 0.000 1.600 1.600 1.600 1.600 4.000 4.000 4.000 4.000 5.600 5.600 5.600 1.600 4.000 5.600 7.838 7.332 6.400 4.800

39.0247 37A635 29.2681 14.0484 0.0000 37.4637 32.7808 24.9759 14.0490 0.0000 35.9026 31.2197 23.4148 12.4879 27.7072 23.0244 15.2195 4.2927 18.3412 13.6584 5.8536 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

Exact 39.0244 37.4634 29.2683 14.0488 0.0000 37.4634 32.7805 24.9756 14.0488 0.0000 35.9024 31.2195 23.4146 12.4878 27.7073 23.0244 15.2195 4.2927 18.3415 13.6585 5.8537 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

Error 0.0003 0.0001 -0.0001 -0.0004 0.0000 0.0003 0.0003 0.0003 0.0002 0.0000 0.0001 0.0002 0.0002 0.0001 -0.0001 0.0000 0.0000 0.0000 - 0.0002 -0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

*Denotes boundary point.

T h e b o u n d a r y conditions for (23) are

f = 0 along th e four edges ( x = + a, y = + b)

(24)

of t h e rect an g l e w i t h sides a and b. T h e solution (23) w i t h b o u n d a r y conditions (24) will be s y m m e t r i c w i t h respect to x an d y, a n d thus, we solve (23) in t h e region (0 < x < a,

288

M. SHARAN ET AL.

x

x

x

x

x

x

x

x

Y

X

FIG. 3.

Location of knots for tile Poisson equation on an ellipse.

0 < y < b) with the b o u n d a r y conditions (i)

f~ = 0 a t x = 0 ,

(ii)

f=0at

(iii)

fv=Oat

(iv)

f=0at

x = a, y=0, y=

(25)

a.

An analytical solution of (23-24) is given by Carslaw and Jeager [23].

f( x, y) = {16Qa2/qra}~{( ×cosh{(2n + -

-1)"cos}(2n

+

1)~ry/2a}/(2n +

Q(a 2 - x2)/2

1)rrx/2a} 1)3cosh{(2n +

1)~'b/2a} (26)

The solution of (23) with (25) can be found similar to (12) with s( x, y) -- Q, fl = 0, and f2 = 0. T h e domain has been normalized to a unit square with a and b. For the computation, we have taken q = 20, D = 2.1978 × 104, a = 20, and b = 20. W e have performed a series of calculations with 25 (i.e., 5 × 5) knots with a uniform spacing A x = A y = 0.25 and with 35 scattered d a t a points. Values of the input parameters, mean square error, and m a x i m u m absolute error are summarized in Table 4. Table 5 shows t h a t the c o m p u t e d solution compares well with the exact solution for 35 scattered d a t a points with E = 2.03 × 10 -~. T h e function varies from - 0 . 1 1 to 0.

289

Solution of Elliptic PDEs TABLE 4 VALUES OF THE INPUT PARAMETERS, MEAN SQUARE ERROR (E), AND MAXIMUM ABSOLUTE ERROR (MAE) FOR DIFFERENT NFOR ~72f = Q

S.N

q

N

R.~

R~

E

MAE

1 2 3 4

20 20 1 1

25 35 25 25

0.5 0.8 0.5 0.7

0.2 0.1 0.2 0.1

1.58 10 -7 223 10 -7 3.98 10 -l° 1.28 10 -9

8.2 1 0 - 4 1.35 1 0 . 4 4.1 10 -5 1.46 10 -4

The m a x i m u m error is found to be 1.35 × 10 -3, which is less t h a n 2.3%. We have also obtained a good solution with N = 35 for q = 1.

3.4.

Poisson Equation in a Rectangle with Right Hand as a Function of Space Coordinates Consider a Poisson equation in a rectangular plate,

(-a

< x < a, - b < y < b):,

+

= g( x, y)

(27)

where g is a known function of x and y. I m p o r t a n t applications occur (a) in fluid dynamics with f = stream function and g = vorticity, (b) in electrostatics, with f = electric potential and g represents the ratio of charge density to dielectric constant, (c) in heat transfer, with f = temperature and g is the l o s s / g e n e r a t i o n of heat, (d) in mass transfer with f = concentration and g is the s o u r c e / s i n k term. Here we considere g of the form of (26). Equation (27) arises as the second Poisson equation in the splitting of the biharmonic equation (see, example (E)) describing the transverse deflection of a thin plate. Again, with the b o u n d a r y conditions, (24), the solution of (27) is expected to be symmetric in x and y and thus, we consider the region (0 < x < a, 0 < y < b ) with the b o u n d a r y conditions (25).

290

M. S H A R A N E T AL.

TABLE 5 COMPUTED AND EXACT SOLUTION OF THE POISSON EQUATION fxx + fyy = Q / D IN WHIC~ q IS A CONSTANT. q = 2 0 , N = 35, R m = 0.9, R n = 0.1 Index

X

Y

*1 *2 *3 *4 *5 *6 *7 *8 *9 *10 *11 "12 "13 "14 "15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

0.000 0.000 0.000 0.000 0.000 0.250 0.500 0.810 0.180 0.500 1.000 1.000 1.000 1.000 0.800 0.130 0.310

0.000 0.220 0.500 0.700 1.000 0.000 0.000 0.000 1.000 1.000 0.000 1.000 0.500 0.700 1.000 0.260 0.420 0.580 0.730 0.180 0.780 0.910 0.5.70 0.820 0.620 0.330 0.370 0.680 0.930 0.110 0.210 0.950 0.440 0.570 0.940

0.770 0.220 0.930 0.860 0.420 0.510 0.680 0.720 0.640 0.840 0.970 0.270 0.530 0.150 0.690 0.480 0.210 0.590

*Denotes b o u n d a r y point.

Computed -0.10644 -0.10222 -0.08387 -0.06006 0.00000 -0.10045 -0.08334 -0.04167 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.09909 -0.08330 -0.03505 -0.05261 - 0.01642 -0.01545 -0.01801 -0.05946 -0.02430 -0.03827 -0.06183 -0.03185 - 0.00481 -0.01555 -0.07978 -0.10062 -0.00766 -0.07151 -0.07358 -0.01061

Exact

Error

-0.10727 -0.10282 - 0.08348 -0.05871 0.00000 -0.10151 -0.08348 -0.04048 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.09961 -0.08349 - 0.03496 -0.05216 -0.01599 -0.01538 - 0.01816 -0.05961 -0.02435 - 0.03824 -0.06152 -0.03147 -0.00477 - 0.01554 -0.07963 -0.10126 - 0.00763 -0.07163 -0.07332 -0.01058

0.00082 0.00060 -0.00039 -0.00135 0.00000 0.00107 0.00014 -0.00119 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00052 0.00018 -0.00009 -0.00046 -0.00044 - 0.00007 0.00015 0.00015 0.00005 -0.00003 -0.00030 -0.00039 -0.00004 -0.00001 -0.00015 0.00064 -0.00003 0.00012 -0.00026 -0.00003

Solution of Elliptic PDEs

291 TABLE 6

VALUES OF THE INPUT PARAMETERS, MEAN SQUARE ERROR (E), AND MAXIMUM ABSOLUTE ERROR ( M A E ) FOR DIFFERENT N FOR V 2 f = Q

S.N

q

N

Rm

R~

1 2 3 4 5 6

20 20 20 1 1 1

25 30 35 25 30 35

1.52 2.8 2.8 1.5 2.8 2.8

0.15 0.1 0.1 0.3 0.1 0.1

E 4.86 5.06 2.11 1.30 1.27 5.28

MAE

10 - 4

10 -~ 10 -5 10 -6 10 -7 10 - s

5.65 1.25 1.48 2.94 6.60 7.38

10 -2 10 -2 10 -2 10 -3 10 -4 10 -4

A n a n a l y t i c a l s o l u t i o n of (27) w i t h b o u n d a r y conditions (24) is o b t a i n e d in t h e form

f( x, y) = ( Q l a ) E [ ( - 1 ) n l ( a n ) ~ ] { 2

-

+ a n y{sinh a J c o s h

[2 + a n b[ {sinh{ a n b } / c o s h { a n b} ] cosh{

hanb}

a n y/} cosh{ a n b}.

(28) W e n o r m a l i z e t h e d o m a i n t o a unit square w i t h a a n d b. T h e values of q, D, a, a n d b are t h e s a m e as in t h e previous example. T a b l e 6 c o n t a i n s t h e values of t h e i n p u t p a r a m e t e r s , m e a n square error, a n d m a x i m u m a b s o l u t e error for q = 20 a n d q = 1 w i t h different values of N. T h e results p r e s e n t e d in T a b l e 7 for q = 20 w i t h N = 35 are v e r y good, considering t h e p o i n t s are s c a t t e r e d . In m o s t of t h e runs ( T a b l e 6), t h e relative p e r c e n t a g e error,

(100 max[ f j -

]j[/)j)

is less t h a n 1%. T a b l e 8 shows t h a t t h e c o m p u t e d values are closed t o t h e a n a l y t i c a l values for q = 1 w i t h N = 30. T h e m a x i m u m error is a p p r o x i m a t e l y 0.02% of its a n a l y t i c a l value. Thus, t h e M Q m e t h o d p r o v i d e s an a c c u r a t e solution to (27) w i t h t h e b o u n d a r y c o n d i t i o n s (24).

292

M. S H A R A N E T AL.

TABLE 7 COMPUTED AND EXACT SOLUTION OF THE POISSON EQUATION fxx + fyy PARAMETERS: q = 20, N = 35, R m = 2.8, R~ = 0.1

+ g( x, y).

Index

X

Y

Computed

Exact

Error

*1 *2 *3 *4 *5 *6 *7 *8 *9 "10 "11 "12 "13 "14 "15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

0.000 0.000 0.000 0.000 0.000 0.250 0.500 0.810 0.180 0.500 1.000 1.000 1.000 1.000 0.800 0.130 0.310 0.770 0.220 0.930 0.860 0.420 0.510 0.680 0.720 0.640 0.840 0.970 0.270 0.530 0.150 0.690 0.480 0.210 0.590

0.000 0.220 0.500 0.700 1.000 0.000 0.000 0.000 1.000 1.000 0.000 1.000 0.500 0.700 1.000 0.260 0.420 0.580 0.730 0.180 0.780 0.910 0.570 0.820 0.620 0.330 0.370 0.680 0.930 0.110 0.210 0.950 0.440 0.570 0.940

9.46105 8.94240 6.83381 4.45365 0.00000 8.79259 6.83966 2.91293 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.57076 6.78130 2.22804 3.83535 1.05650 0.78539 1.13562 4.35235 1.40578 2.46047 4.61329 2.10201 0.23959 1.00667 6.43731 8.75920 0.38783 5.53565 5.78024 0.58868

9.46366 8.94669 6.84478 4.46840 0.00000 8.79692 6.84478 2.91812 0:00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.57616 6.78721 2.22966 3.84322 1.05684 0.78637 1.13662 4.35595 1.40709 2.46234 4.61689 2.10327 0.23987 1.00844 6.44232 8.76403 0.38836 5.53995 5.78855 0.58942

-0.00261 - 0.00428 -0.01097 -0.01476 0.00000 -0.00434 - 0.00512 -0.00519 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.00540 - 0.00591 -0.00162 -0.00787 -0.00035 - 0.00098 - 0.00101 -0.00360 -0.00131 - 0.00187 -0.00361 -0.00126 - 0.00028 -0.00177 -0.00501 - 0.00482 - 0.00052 -0.00429 -0.00831 -0.00074

*Denotes b o u n d a r y point.

Solution of Elliptic PDEs

293 TABLE 8

COMPUTED AND EXACT SOLUTION OF THE POISSON EQUATION, fzx + fyy = g( X, y);

PARAMETERS: q = 1, N = 30, Rm = 2.8, R n = 0.1 Index

X

Y

Computed

Exact

Error

"1 *2 *3 *4 *5 *6 *7 *8 *9 "10 *11 "12 "13 "14 "15 "16 "17 "18 19 20 21 22 23 24 25 26 27 28 29 30

0.000 0.000 0.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 0.200 0.200 0.400 0.400 0.600 0.600 0.800 0.800 0.050 0.130 0.460 0.310 0.070 0.120 0.420 0.510 0.680 0.840 0.970 0.170

0.000 0.250 0.500 0.750 1.000 0.000 0.250 0.500 0.750 1.000 0.000 1.000 0.000 1.000 0.000 1.000 0.000 1.000 0.050 0.260 0.160 0.420 0.580 0.730 0.910 0.570 0.820 0.370 0.680 0.930

0.47259 0.43932 0.34184 0.18864 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.45114 0.00000 0.38796 0.00000 0.28628 0.00000 0.15261 0.00000 0.46991 0.42831 0.35106 0.33905 0.29663 0.19946 0.05679 0.21766 0.07033 0.10504 0.01199 0.05292

0.47318 0.43985 0.34224 0.18907 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.45177 0.00000 0.38862 0.00000 0.28695 0.00000 0.15323 0.00000 0.47051 0.42881 0.35154 0.33936 0.29698 0.19977 0.05683 0.21780 0.07035 0.10516 0.01199 0.05305

- 0.00059 - 0.00052 - 0.00040 -0.00043 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.00063 0.00000 -0.00065 0.00000 - 0.00066 0.00000 -0.00062 0.00000 - 0.00059 -0.00050 -0.00048 - 0.00031 - 0.00035 -0.00031 -0.00004 - 0.00014 - 0.00003 - 0.00012 0.00000 -0.0001

*Denotes b o u n d a r y point.

3.5.

Biha~monic Equation

T h e t r a n s v e r s e d e f l e c t i o n ( f ) of a t h i n r e c t a n g u l a r p l a n e w i t h s i d e s 2 a a n d 2b, of f l e x u r a l r i g i d i t y D, s u b j e c t t o a n o r m a l l o a d q p e r u n i t a r e a is governed by, V4f=

Lx:~ + 2f~yy

+ fyuuy = q / D ,

( a < x < a, - b < y < b). (29)

294

M. SHARAN ET AL.

Equation (29) is subject to the b o u n d a r y conditions f = 0 and fun = 0 along its four edges, where ~? denotes the normal to the boundary. By introducing the variable u = V2f, (29) a m o u n t s to solving Poisson's equation twice in succession

(i) u x ~ + u y y = q/D; - a < x < a, - b < y < b,

(30)

with u = 0 along the four edges, and (ii) fxx + f y y =

u ( x , y); - a <

x < a, - b <

y<

b

(31)

with f = 0 along the four edges x = + a and y = + b.

Notice t h a t (30) and (31) are similar, but not identical, to (23) and (27). In order to obtain f from (31), the solution (30) will be used as the input. The relation (28) is the analytical solution of the biharmonic equation, (29). The values of q, D, a and b are the same as given in Examples 3.3 and 3.4. As in the previous cases, the domain was normalized to a unit square with reference to a and b. T h e values of the input parameters, mean square error and the m a x i m u m absolute error for different values of N and q are given in Table 9. T h e solution of the biharmonic equation, (29) for q = 20 is presented in Table 10 for 35 scattered d a t a points. There is a good agreement between the c o m p u t e d and analytical values (Table 9). The c o m p u t e d values are found to be close, see Table 11, to the exact values for q = 1 with 25 uniform knots. Notice t h a t the accuracy of the solution of biharmonic equation c o m p u t e d in two stages using (30 and 31) and in a single stage from (27) is found to be almost same.

TABLE 9 VALUES OF THE INPUT PARAMETERS, MEAN SQUARE ERROR ( n ) , AND MAXIMUM ABSOLUTE ERROR (MAE) FOR DIFFERENT N FOR V 4 f = Q

S.N

q

N

Rm

Rn

Error

MAE

1 2 3 4

20 20 1 1

25 35 25 35

2.2 2.8 1.5 1.5

0.15 0.1 0.3 0.3

3.89 10 -4 1.58 10 -4 1.20 10 -6 1.72 10 -7

5.22 10 -2 3.49 10 -~ 2.85 10 -3 1.26 10 -3

Solution of Elliptic PDEs

295

TABLE 10 COMPUTED AND EXACT SOLUTION OF THE BIHARMONIC EQUATION

q = 20, N = 3 5 , R , ~ = 2 . 8 , R n = 0 . 1 Index

X

Y

Computed

Exact

Error

1" 2* 3* 4* 5* 6* 7* 8* 9* 10" 11" 12" 13" 14" 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

0.000 0.000 0.000 0200 0.000 0.250 0.500 0.810 0.180 0.500 1.000 1.000 1.000 1.000 0.800 0.130 0.310 0.770 0.220 0.930 0.860 0.420 0.510 0.680 0.720 0.640 0.840 0.970 0.270 0.530 0.150 0.690 0.480 0.210 0.590

0.000 0.220 0.500 0.700 1.000 0.000 0.000 0.000 1.000 1.000 0.000 1.000 0.500 0.700 1.000 0.260 0.420 0.580 0.730 0.180 0.780 0.910 0.570 0.820 0.620 0.330 0.370 0.680 0.930 0.110 0.210 0.950 0.440 0.570 0.940

9.42875 8.91823 6.82939 4.45867 0.00000 8.76810 6.82994 2.91433 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.55057 6.77241 2.22803 3.83735 1.05840 0.78568 1.13516 4.34914 1.40524 2.45991 4.61198 2.10377 0.23966 1.00727 6.43005 8.7366 0.38784 5.53038 5.77787 0.58848

9.46366 8.94669 6.84478 4.46840 0.00000 8.79692 6.84478 2.91812 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.57616 6.78721 2.22966 3.84322 1.05684 0.78637 1.13662 4.35595 1.40709 2.46234 4.61689 2.10327 0.23987 1.00844 6.44232 8.76403 0.38836 5.53995 5.78855 0.58942

- 0.03491 -0.02845 -0.01539 - 0.00973 0.00000 -0.02882 -0.01484 -0.00379 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.02559 -0.01480 -0.00163 - 0.00587 0.00156 -0.00070 - 0.00147 -0.00680 -0.00186 -0.00243 -0.00492 0.00050 -0.00021 -0.00117 - 0.01227 -0.02736 -0.00051 -0.00956 - 0.01068 -0.00093

*Denotes boundary point.

M. SHARAN ET AL.

296

TABLE 11 COMPUTEDAND EXACTSOLUTIONOF THE BIHARMONICEQUATION. PARAMETERS: q = 1, N = 25, R m = 1.5, R , = 0.3 Index

X

Y

Computed

Exact

Error

"1 *2 *3 *4 *5 6* 7 8 9 *10 "11 12 13 14 "15 "16 17 18 19 *20 "21 *22 *23 *24 *25

0.000 0.250 0.500 0.750 1.000 0.000 0.250 0.500 0.750 1.000 0.000 0.250 0.500 0.750 1.000 0.000 0.250 0.500 0.750 1.000 0.000 0.250 0.500 0.750 1.000

0.000 0.000 0.000 0.000 0.000 0.250 0.250 0.250 0.250 0.250 0.500 0.500 0.500 0.500 0.500 0.750 0.750 0.750 0.750 0.750 1.000 1.000 1.000 1.000 1.000

0.47075 0.43843 0.34288 0.19193 0.00000 0.43843 0.40735 0.31876 0.17707 0.00000 0.34288 0.31876 0.24853 0.13783 0.00000 0.19193 0.17707 0.13783 0.07627 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

0.47318 0.43985 0.34224 0.18907 0.00000 0.43985 0.40895 0.31839 0.17606 0.00000 0.34224 0.31839 0.24836 0.13772 0.00000 0.18907 0.17606 0.13772 0.07675 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

-0.00243 - 0.00142 0.00064 0.00285 0.00000 -0.00142 - 0.00060 0.00037 0.00101 0.00000 0.00064 0.00037 0.00018 0.00011 0.00000 0.00285 0.00101 0.00011 -0.00048 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

*Denotes boundary point.

4.

DISCUSSION

W e have successfully applied the M Q m e t h o d to c o m p u t e the solution of elliptic partial differential equations with Dirichlet a n d / o r N e u m a n n boundary conditions. T h e d a t a points in a r b i t r a r y locations with an a r b i t r a r y ordering allows us to o b t a i n a n accurate solution. T h e m e t h o d is shown to work well for the Laplace, Poisson a n d b i h a r m o n i c equations. W h e n comparing the MQ c o m p u t e d values with the exact solution, the agreement is found to be very good. F o r the s u m m a t i o n of the series in (18), (26), a n d (28), sufficiently large n u m b e r of terms ( n ) are considered such t h a t the absolute value of n t h t e r m is less t h a n or equal to 10 -7.

Solution of Elliptic PDEs

297

In the first example, the error was relatively higher on the knots near the point of intersection of the sides with different temperature because of the discontinuity at the boundary. The MQ method performs well for the problems with a curved boundary. The curved boundary is accounted here in a natural way through the location of the data points, unlike the finite difference method in which an interpolation technique is required to express a boundary point in terms of knots. Dubal [15, 16] also pointed out the advantages of using MQ to solve a highly nonlinear Poisson equation representing the three-dimensional initial states of black hole collisions since MQ permits boundary points to be placed exactly on the boundary of the event horizons with great ease. 5.

CONCLUSION

Because the papers published dealing with the applications of MQ to the solutions of PDEs have been primarily exploratory studies, the MQ method lacks the technical maturity of the well-studied finite difference, finite element, finite volume, or pseudo-spectral methods that have been studied for many years by many people. This lack of technical maturity is a hinderance to immediate production code applications. Like so many of the mature schemes, MQ or other radial basis functions require studies to exploit its many advantages. We will highlight the areas of exploration that will enable a powerful tool like MQ radial basis functions to be applied in production codes. Some of the areas in which further research is required are 1. 2. 3. 4.

5.1.

Minimization of the computational effort and potential ill-conditioning; Selection of the shape-parameters and knots; Selection of distance metrics; Solution of PDE problems for which no analytic solution exists.

Minimization of the Computational Effort and Potential Ill-conditioning We have shown that MQ is an excellent approximation scheme for solving the biharmonic equation. The solution of this equation obtained in two stages is approximately the same as obtained in a single stage. Dubai et al. [16] have suggested the use of the singular value decomposition [26] for solving the resulting system of algebraic equations. However, we have observed that in the examples considered here, Gaussian elimination with total pivoting gives better results than the singular value decomposition. Because the coefficient matrix is full, Kansa [19] has suggested domain decomposition and blending. Foley [27] later presented a general formulation

298

M. SHARAN ET AL.

of the map and blend technique and applied it over the spherical earth. Dubal [28] showed that domain decomposition is straightforward to implement, the computational effort is considerably decreased, and problems of ill-conditioning are alleviated.

5.2. Selection of the Shape-parameters and Knots Quite often, a slight variation in R~ and R m may significantly change the mean square errors. Further, there is no systematic way to determine the choice of R n and R m like the finite difference methods where for consistent and stable approximation scheme, a decrease in the grid spacing yields to a more accurate results. More research is required to determine whether this sensitivity is due to ill-conditioning, which is probably the most likely cause of the observed sensitivity. The drawback of finite difference methods is the expensive mesh refinement required in steep gradient regions which can become prohibitive in three dimensions. Thus, it is necessary to develop a suitable procedure for estimating Rn and R m in terms of function values, N, and the location of knots. The study by Carlson and Natarajan [29] showed that sparse approximates of MQ interpolation can be obtained if one permits a small arbitrary error in the interpolation using the greedy algorithm. Franke et al. [30, 31] have applied the greedy algorithm in selecting a minimum number of knots in optimizing accuracy. Relatively few knots are required using MQ interpolation if the locations of the knots can be optimally chosen. We believe that the problem of minimum knots and shape-parameter selection are somewhat connected to the curvature of a function, the loci of the extrema, and inflection points. Franke [7] showed in his example with the Franke function 1, over a 25-point data set, that if the region near an extremum was undersampled, neither MQ, or any other interpolation scheme could properly predict the behavior of the interpolated function; when a data set such as the 33 point or 100 point data set was used which properly sampled the function, the results with MQ were significantly improved. The optimal location of knots was further explored by Kansa and Carlson [32] who showed that very rapidly varying univariate functions could be accurately interpolated with just seven knots if these knots in the interior of the domain were optimally located. Kansa [8, 9] has shown that the power law seems to work very well. Hagen and Kansa [9] use the Brent algorithm to find the optimum constant r or the shape parameters R~ and R m. We have also made an attempt to obtain these shape parameters using the nonlinear optimization technique suggested by Marquardt [25] by choosing these values by experimentation. Quite often, the solution is very sensitive to and a slight variation in R~ and R m (ill-conditioning) may change significantly the mean square errors.

Solution of Elliptic PDEs

299

However, Golberg et al. [21] have developed systematic way to determine the choice of the shape parameter by cross-validation that guarantees very good predictive capabilities.

5.3. Selection of Distance Metrics Another area of exploration is a generalization of the form that the MQ basis function regarding the distance metric. In most applications, the distance in the MQ basis function is assumed to be the Euclidean distance, but other metrics are also possible. A possible form of a generalized MQ basis function is g(x-xj)

=ae{(x-xj)

T

.i-(x-xj)

2~ 1/2

+ rj]

(32)

where t is a rank two-tensor. If, for certain values of x, the basis function becomes imaginary, g is set to zero. Such transformation was explored by Kansa and Carlson [32] and Girosi [33] to obtain better MQ interpolants and approximations. By using only the real projection of the hyperboloids, the MQ basis function, centered at different knots, can have either global or local support. Note, that form of (32) has much in common with wavelets; this connection to wavelets has been partially explored theoretically by Buhmann [34]. It quite trivial to obtain an orthonormalized basis set using the Gram-Schmidt procedure. The connection among the optimization of the knot locations, shape parameters, and metric tensor will be discussed in the next paragraph.

5.4. Solution of PDE Problems for which No Analytic Solution Exists In many previous exploratory studies of MQ approximation in the numerical solution of PDEs, the determination of the parameters such as shape parameters, knots, metric, and amount of rotation relied heavily upon an a priori knowledge of the exact solution. These studies relied upon collocation methods that obtained solutions at a specified set of knots. Recently, Galperin and Zheng [35] and Galperin [36] have recommended a method based on global optimization techniques which may be used in case the exact solution is not known. Their technique is very general and can be applied to finite difference, finite volume, finite element, spectral, or radial basis functions to solve a general set of PDEs. They do not rely upon collocation, but rather globally minimizing errors over the domain and boundaries. Any linear or nonlinear PDE problem can be cast in the following form L(f) = 0

(33)

300

M. SHARAN ET AL.

subject to boundary condition constraints, a(f)

= O,

(34)

and if, appropriate, initial condition constraints, • ( f ) = O.

(35)

One choice of cost-functions, recommended by Galperin and Zheng [35] is the minimization of the integral of the PDE over each subdomain having the following form for elliptic PDEs

C( f ) = qe Qm{ni f L( f ) dv + vfr(r ( f ) dS I

(36)

where C is the cost function, q is a vector of parameters (shape function, knot location, tensor parameters, and expansion coefficient) from the space of admissible parameters, Q; dv and dS represent the differential volume and surface elements, respectively. Galperin and Zheng [35] showed that the minimization of (36) can yield very accurate results with very few MQ basis functions. Hence, they believe that the excellent convergence properties of MQ make it very competitive with familiar low order finite difference and finite element schemes, especially in higher spatial dimension PDE problems. REFERENCES 1 A.H. Nayfeh, Perturbation Method~ Wiley Interscience, New York, 1973. 2 M. Van Dyke, Perturbation Methods in Fluid Mechanic~ Parabolic Press, Palo Alto, 1975. 3 J. Kevorkian and J. D. Cole, Perturbation methods in Applied Mathematica Springer-Verlag, New York, 1981. 4 J.Y. Parlange, Theory of water movements in soils I: One-dimensional absorption, Soil Sci., 111, 134-170 (1971). 5 T. Tsang, Transient state heat transfer and diffusion problems, Ind. Eng. Chem. 52:707 (1960). 6 R.L. Hardy, Theory and applications of the multiquadric biharmonic method: 20 Years of discovery 1968-1988, CompuL Matk AppL 19 (8/9):163-208 (1990). 7 R. Franke, Scattered data interpolation: tests of some methods, Matk Comp. 38:181-200 (1982).

Solution of Elliptic PDEs

301

8 E.J. Kansa, Multiquadrics: a scattered data approximation scheme with application to computational fluid dynamics: 1. surface approximations and partial derivative estimates, CompuL Math. Applic. 19, (8/9):127-154 (1990). 9 E.J. Kansa, Multiquadrics: a scattered data approximation scheme with application to computational fluid dynamics: II parabolic, hyperbolic, and elliptic partial differential equations, CompuL Math, Applic., 19, (8/9):146-161 (1990). 10 C.A. Micchelli, Interpolation of scattered data: distance matrices and conditionally positive definite functions, Const Approx., 2:11-22 (1986). 11 B.C. Baxter, The asymptotic cardinal function of the multiquadric, ](r) = ( r 2 + c2) as c2/A~ • Comput. Math, Applic., 24 (12):1-6, (1992). 12 W.R. Madych, Miscellaneous error bounds for multiquadric and related interpolators, CompuL Math, Applic. 24 (12):121-138 (1992). 13 W . R . Madych and S. A. Nelson, Multivariate interpolation and conditionally positive definite functions Approx. Theo. App£ 4:77-89 (1989). 14 W.R. Madych and S. A. Nelson, Multivariate interpolation and conditionally positive definite functions II, Math, Comput 54:211-230 (1990). 15 M. R. Dubal, Construction of three-dimensional black-hole initial data via multiquadrics, Phys. Rev. D 45:117-118 (1992). 16 M.R. Dubal, S. R. Oliverir£, and R. A. Matzner, Solution of elliptic equations in numerical relativity using multiquadrics, In Approaches to Numerical Relativity R. D. Inverno, ed.), Cambridge University Press, Cambridge, England, 1992. 17 G. Moridis and E. J. Kansa, The Laplace transform multiquadric method highly accurate scheme for numerical solution of partial differential equations, J. Appl Sc~ and Comput 1(2):375-407 (1994). 18 R.E. Carlson and T. A. Foley, The parameter R 2 in multiquadric interpolation, Comput Math, Applic. 21:29-42 (1991). 19 R.E. Hagen and E. J. Kansa, Studies of the R parameter in the multiquadric function applied to ground water pumping, J. Appl Sci. and CompuL 1(2):266-281 (1994). 20 D.M. Marquadt, An algorithm for least squares estimation of nonlinear parameters, J. Soc. Ind. AppL Math, 11:431-441 (1963). 21 M.A. Golberg, C. S. Chen, and S. R. Karur, Improved multiquadric approximation for partial differential equations, preprint 22 A. Makroglou, "Radial basis functions in the numerical solution of Fredholm integral and intergrodifferential equations," In Advances in Computer Methods for Partial Differential equations, VII, 478-484, Editors: R. Vichnevetsky, D. Knight, and G. Richter, Proc. of the 7th IMACS conference on Computer Methods for PDEs, New Brunswick, N J, USA, June 22-24 (1992). 23 H.S. Carslaw and J. C. Jaeger, Heat Conduction in Solids, Oxford University Press, Oxford, 1976. 24 W . F . Hughes and E. W. Gaylord, Basic equations of Engineering Science~ Schanm Publishing Company, New York, 1964. 25 Maithili Sharan, M. P. Singh, and I. Sud, Modelling of oxygen transport in systemic circulation in a hyperbaric environment. Frontiers Meal BioL Eng., 3:2.-3' (1991).

302

M. SHARAN ET AL.

26 W.H. Press, B. P. Flanncry, S. A. Teukolsy, and W. T. Vetterling, Numerical Recipe~ Cambridge University Press, Cambridge, 1986. 27 T. A. Foley, The map and blend scattered data interpolation on a sphere, CompuL Matl~ Applic. 24(12):41-60 (1992). 28 M. R. Dubal, Domain decomposition and local refinement for multiquadric approximations. I: Second-order equations in one dimension, J. AppL Sci~ and CompuL, 1(1):146-171 (1994). 29 R.E. Carlson and B. K. Natarajan, Sparse approximate multiquadric interpolation, CompuL Matl~ Applic. 27(6):99-108 (1994). 30 R. Franke, H. Hagen, and G. M. Nielson, Least squares surface approximation in scattered data using multiquadric function, Adv. CompuL Math. 2:81-99 (1994). 31 R. Franke, H. Hagen, and G. M. Nielson, Repeated knots in least squares multiquadric functions, Computing SuppL 10:(1995). 32 E.J. Kansa and R. E. Carlson, Improved accuracy of multiquadric interpolation using variable shape parameters, CompuL Math. AppIic. 24(12):99-120 (1992). 33 F. Girosi, Some extensions of radial basis functions and the applications to artificial intelligence, CompuL Math. Applic. 24(12):61-80 (1992). 34 M. D. Buhmann, Discrete least squares approximation and prewavelets for radial function spaces, DAMTP 1993/NA4, University of Cambridge, Cambridge, England (1993). 35 E.A. Galperin and Q. Zheng, Solution and control of PDE via Global optimization methods, Compu~ Math Applic. 25(5):103-11 (1993). 36 E. A. Galperin, The cubic algorithm for optimization and control, ND Research'e Pubt Montreal 1990.