Motivation for using Radial Basis Functions to solve PDEs

Edward J. Kansa, Ph.D.


Lawrence Livermore National Laboratory and Embry-Riddle Aeronatical University
Mail Stop L-200 East Bay Campus, Oakland CA 94621,
Livermore, CA 94551-0808 USA Daytona Beach, FL 32114-3900,
Telephone : 925-423-0151
Fax: 925-423-6907
E-mail : kansa1@llnl.gov




1, Introduction

The numerical solution of partial differential equations (PDEs) has been dominated by either finite difference methods (FDM), finite element methods (FEM), and finite volume methods (FVM). These methods can be derived from the assumptions of the local interpolation schemes. These methods require a mesh to support the localized approximations; the construction of a mesh in three or more dimensions is a non-trivial problem. Typically with these methods only the function is continuous across meshes, but not its partial derivatives.

In practice, only low order approximations are used because of the notorious polynomial snaking problem. While higher order schemes are necessary for more accurate approximations of the spatial derivatives, they are not sufficient without monotonicity constraints. Because of the low order schemes typically employed, the spatial truncation errors can only be controlled by using progressively smaller meshes. The mesh spacing, h, must be sufficiently fine to capture the functions partial derivative behavior and to avoid unnecessarily large amounts of numerical artifacts contaminating the solution. Spectral methods while offering very high order spatial schemes typically depend upon tensor product grids in higher dimensions.


2, Scattered Data Interpolation

Let x1,...,xNÎ<W Ì< Rn be a given set of nodes. Let

gj(x) º< g( || x-xj|| ) Î< R,      j = 1,...,N,

be a set of any RBF basis functions. Here || x-xj|| is the Euclidean distance. Given interpolation data values y1,...,yN Î<R at data locations x1,...,xN Î< W Ì< Rn, the RBF interpolant

F(x) = N
å<
j = 1
 ajgj(x)+aN+1
(1)

is obtained by solving the system of N+1 linear equations

N
å<
j = 1
 ajgj(xi)+aN+1
=
yi,      i = 1,...,N,
(2)
N
å<
j = 1
 aj
=
0,

for N+1 unknown expansion coefficients aj. This formulation of interpolation problem, where is due to R. L. Hardy, Hardy (1972,1990), the primary innovator of the RBF method. He adds a constant to the expansion and constrains the sum of the expansion coefficients to be zero. Introducing the notation

p =




1
:
1





Î< RN,      G =




g1(x1)
...
gN(x1)
:
:
:
g1(xN)
...
gN(xN)





Î< RN×N,      H =

G
p
pT
0


Î< R(N+1)×(N+1),

a = ( a1,...,aN+1) T, y = (y1,...,yN,0) T Î< RN+1, we can rewrite the system (2) in the matrix from as

Ha = y.
(3)

Then the interpolation expansion coefficients are given by

a = H-1y.
(4)

One can also easily find the derivatives of the interpolant at the nodes xi, e.g.

F¢<(xi)
=
N

j = 1
 ajgj¢<(xi),      i = 1,...,N,
(5)
F¢¢(xi)
=
N

j = 1
 ajgj¢¢(xi),      i = 1,...,N.
(6)

There is an infinite class of radial basis functions possible. Some of the more commonly used RBFs are:

Thin-plate splines, || x-xj|| ln( ||x-xj|| ) ,
(7)
Linear splines, || x-xj|| ,
(8)
Cubic splines, || x-xj|| 3,
(9)
Gaussians, exp

- || x-xj||
cj2

,
(10)
Multiquadrics (MQ),

1+ || x-xj||
cj2
 2


.
(11)

Equations (7) - (11) represent some of the commonly used global RBFs that are defined over all W. Splines ( 7) - (9) have no adjustable parameters, whereas the Gaussian and MQ do. Hardy (1990) showed that MQ has a physical foundation; it is related to a consistent solution of the biharmonic potential problem. Buhmann and Micchelli (1992) and Chui et al. (1996) have shown that RBFs are related to wavelets, except for the orthonormalization property. The translational invariance property is accounted by ||x-xj|| . RBFs (7) - (9) have uniform scaling, whereas the Gaussians and MQs can have local dilations via the local shape parameter, cj2. The neural network community which relies upon Gaussian RBFs extensively optimizes the local shape parameter to accelerate the convergence of their RBF expansions. Kansa (1990a, 1990b) and Kansa and Hon (1999) found heuristically that even a crude recipe for the shape parameter distribution accelerates convergence in collocation PDE problems using MQ as the RBF.

The reader is referred to Wendland (1995) for the polynomial compacted supported (CS-RBF) splines and to Kansa and Hon (1999) for the truncated MQ splines.


3, RBF PDEs by Collocation (Asymmetric Method).

The procedure will be illustrated with an elliptic PDE. The spatial treatment for hyperbolic and parabolic PDEs is formulated similarly. Time dependent PDEs can be solved by the method of lines (MOL) technique. In this case one substitutes the RBF expansions of spatially dependent terms which yields a system (or systems) of coupled linear or nonlinear ordinary differential equations (ODEs).

a) Solution procedure for elliptic PDEs.

For elliptic PDEs, we assume the interior and boundary operators, L and B, respectively, are linear, for simplicity, and that they define a well-posed elliptic boundary value problem:

Lu(x) = s(x),     in W Ì< Rn,
(12)
Bu(x)| ¶W = f(x).
(13)

We first introduce a set Qh of nodes

Qh = { (xi)\mid i = 1,N-M Ì< W,      (xi)\mid i = N-M+1,N Ì< ¶W}

We look for the approximate solution uh to (12), ( 13) in the form

uh(x) = N
å<
j = 1 <
 ajgj(x)+aN+1,       x Î<
W
 <
 = WȶW,
(14)

Substituting uh(x) into (12), (13) and using collocation at the nodes in Qh, we obtain the collocation system

Luh(xi)
=
N
å<
j = 1 <
 ajLgj(xi)+aN+1L1 = s(xi),      i = 1,...,N-M,
(15)
Buh(xi)
=
N
å<
j = 1 <
 ajBgj(xi)+aN+1B1 = f(xi),      i = N-M+1,...,N.
(16)

Here L1 and B1 mean the action of L and B on the number 1. Introducing the notation:

GL
=





Lg1(x1)
...
LgN(x1)
:
:
:
Lg1(xN-M)
...
LgN(xN-M)





,      GB =




Bg1(xN-M+1)
...
BgN(xN-M+1)
:
:
:
Bg1(xN)
...
BgN(xN)





,
W
=





GL
Lp
GB
Bp
pT
0





Î< R(N+1)×(N+1),      b = (s(x1),...,s(xN-M),f(xN-M+1),...,f(xN),0) T Î< RN+1,

we can rewrite the system (15), (16) in the matrix from as

Wa = b,
(17)

whose solution is

a = W-1b.
(18)

Then the solution anywhere in the domain is given by (14). If L is a nonlinear operator, then the solution to (17) will require some iterative procedures.

b) Parametrized nonlinear elliptic PDEs.

Typically, in nonlinear problems it is crucial to understand the qualitative dependence of the solutions on the problem parameters. In Fedoseyev, Friedman, and Kansa (1999) some parametrized nonlinear 1D and 2D elliptic PDEs were discretized by the MQ method. Then numerical continuation and bifurcation techniques were used to study qualitative changes in solutions as a problem parameter varied. A significant improvement of the accuracy was observed when a simple optimization strategy with 2 parameters was used. These parameters were 1) the constant shape parameter c = cj for all j, 2) distance [h\tilde] between the boundary ¶W and the layer of the nodes adjacent to ¶W.

c) Solution procedure for hyperbolic and parabolic PDEs.

Assume that in some inertial frame, the time and space variables are separable, and the expansion coefficients depend only upon time. Then the solution procedure for time dependent hyperbolic and parabolic PDEs is similar for the spatial discretization part. The procedure is somewhat modified.

Firstly, the expansion coefficients at time t = 0 are obtained as a solution to the interpolation problem from the initial conditions.

Secondly, the time evolution of the expansion coefficients is generated by the method of lines yielding a coupled system of linear or nonlinear ordinary differential equations(ODEs). Examples can be found in which explicit 4th order Runge-Kutta, Laplace Transform, or even exact ODE solution methods were employed yielding excellent results.

In time marching schemes, it is important to control truncation errors. For example, in the case of hyperbolic PDEs spatial and truncation errors propagate along characteristics. Temporal truncation errors at a node, xk , can contaminate the spatial solution, and vice versa. The reduction of both spatial and temporal truncation errors was demonstrated to be excellent using a combination of Laplace Transforms and MQ, see Moridis and Kansa (1994). Furthermore, in most instances, it was not found necessary to add stabilizing artificial viscosity to advective problems (i.e., upwind differencing methods common with FDM, FEM, or FVM).


4, Hermite-Birkhoff collocation PDE Method (Symmetric Method)

Wu (1996) viewed the solution of the PDE problem as a special case of the Hermite-Birkhoff interpolation problem. Although extra steps are required, the PDE collocation matrix that is generated from the Hermite-Birkhoff method is positive definite, hence the problem always solvable. The Hermite-Birkhoff formulation of the PDE problem is called the symmetric approach. Franke and Schaback (1998) and Fasshauer (1998) used Wendland (1995) compactly supported radial basis functions (CS-RBF).

Wendland's CS-RBFs permit one to adjust the band width of the collocation matrix. Franke and Schaback observed that if the matrix band width were made more narrow, convergence could be only accelerated by using more nodes. But as the CS-RBFs became more global, hence larger matrix band widths, the convergence accelerated.


5, Optimization of Node Location and MQ shape parameters.

Regardless of the spatial approximation scheme used, the more interesting PDE problems require an intelligent placement of the nodes, {xi} , to avoid the under-sampling of important physical phenomena. For example, Hon (1998), using MQ to solve Burgers viscid equations, found that as the Reynolds number increased, an adaptive choice of a subset of nodes was required to resolve the steepening front while using as few nodes as possible. Kansa (1990 b) and Kansa and Hon (1999) optimized the MQ shape-parameter distribution, { cj2} , and observed a twofold effect: 1) the condition number of the W matrix decreased by several orders of magnitude, and 2) the errors in the PDE solution decreased by a couple orders in magnitude.


6, Collocation-free global optimization formulation.

Galperin and Zheng (1993) and Galperin, Pan, and Zheng (1993) argue that all collocation methods are intrinsically ill-conditioned. They proposed to solve PDEs using global optimization methods, and noted many advantages of the MQ-RBFs. The set of initial conditions, boundary conditions on ¶W, and PDE in W are cast into functionals. In turn, a global functional is formed using a weighted sum of these functionals over the initial conditions, boundary conditions, and interior problem. A vector, q = { xi,cj2,aj} in the space, Q, of feasible parameters is varied, until discrepancies are £< h, the prescribed precision. It is noteworthy to mention that Galperin and Zheng found the global optimized solution to a heat conduction PDE to a precision of 0.0012 with just one optimized MQ basis function.

The global optimization approach circumvents any ill-conditioning problems arising from either the asymmetric or symmetric collocation PDE problem. Ill-posed and badly formulated problems can possess h-equivalent solutions that represent physical reality despite the mathematical nonexistence of an exact solution.

Galperin, Kansa, Makroglou, and Nelson (1999) applied global optimization to the numerical solutions of weakly singular Volterra integral equations. They used a method of successive approximations, appending one MQ basis function at a time (j = 1,2,...,m), while optimizing each {xj,cj2,aj} , (a 3-parameter optimization procedure). They observed that four to seven MQ basis functions were required for convergence with the error £< 5×10-7, depending upon the problem.


7, Convergence Properties of the MQ interpolation scheme.

Madych and Nelson (1990) proved that the MQ interpolant and its partial derivatives have exponential convergence. This rate of convergence has been observed numerically by a number of authors. Furthermore, Buhmann and Micchelli (1992) and Chui et al. (1996) showed that RBFs are prewavelets (wavelets are ortho-normalized).

MQs are especially attractive because they have exponential convergence and depend upon the shape parameter distribution that can be globally optimized as an h-equivalent problem, see Galperin and Zheng. Optimization of shape parameters is commonly accepted in neural network problems in which Gaussian RBFs rely upon optimization of the variance of the basis functions.


8, Lagrangian and Moving Node applications with RBFs.

Lagrangian and other moving node schemes used for computational fluid dynamics problems are quite complicated to implement in two and three dimensions, especially if the underlying mesh undergoes considerable distortion and breaking of connectivity. The natural tendency (due to entropy) of such systems is towards a scattered data arrangement. Richard Franke's paper (1982) concluded that of all the interpolation methods for scattered data problems, Hardy's MQ-RBF splines (1972) performed best based upon several criteria. While there are several RBF splines available, in general, MQ seems to have outperformed others.

One of the perennial problems of numerical mathematics is the problem of simulating problems whose solutions are continuous in certain regions and discontinuous in other regions. It is proposed here that discontinuities be tracked, not captured. MQ is a continuously differential RBF. Shocks and contact surfaces are discontinuities of lower dimensionality, and except in the limit, continuous functions do not represent discontinuities well.

If a continuous wave begins to steepen to such an extent that its width in the normal direction of propagation tends towards zero, then replace that wave by a discontinuity of lower dimensionality. Front tracking, wave steepening, and front collisions have been fairly simple to implement in 1D with mesh-based methods, but require complex bookkeeping, transformations, etc. if mesh-based methods are used. A huge set of complications is eliminated if a mesh is abandoned.

A Lagrangian FDM scheme tries to minimize the remapping back onto a well-behaved Eulerian mesh because the interpolation scheme is very diffusive. While the interpolation is performed in such a manner to be strictly mass conservative, conservation of momentum and total energy are typically sacrificed. Recent experiments with MQ-RBFs have shown many of the shortcomings with traditional methods can be avoided.

The method of characteristics (MOC) cannot be rigorously constructed beyond one dimension. However, by simultaneously rotating the fluid velocity components and the coordinate system into the normal flow direction, the resulting PDEs are locally one-dimensional. Hence, the MOC can be used to find the advanced time solution.

Remapping is enforced whenever two or more nodes begin to coalesce or leave regions devoid of nodes. Remapping is performed by a modification of the basic interpolation method. The MQ basis functions are integrated in closed form. Because the same spatial basis function is used for all variables, the total mass, momentum component and total energy integrals over the domain, W, are calculated merely by changing the expansion coefficients. The remapping step is constrained to be strictly conservative in all dependent variables by a least squares procedure that combines interpolation with the conservation constraint.


9, Invertible MQ PDE collocation coefficient matrices.

Micchelli (1986) proved that the interpolation problem related to RBFs such a thin-plate-splines and MQ is solvable for distinct nodes. To date, there is no theory that specifies the conditions under which the asymmetric MQ collocation is solvable. However, from experience, the following conditions should be observed.

  1. Nodes should be distinct.
  2. Nodes from W should not coincide with those of the boundary, ¶W.
  3. The problem should be well-posed.
  4. Certain unusual combinations of nodes and shape parameters may produce W matrices in which two or more rows may be nearly linearly dependent or even singular, see Schaback and Hon (1999). To avoid such a situation, it is recommended that the matrix be regularized by perturbing the shape parameter distribution, and avoiding very large values of shape parameters that would result in such a situation. Note that in the global optimization formulation, such combinations of nodes and MQ shape-parameters would not be permitted.
  5. A variable shape parameter distribution gives more distinct rows of the coefficient matrix that considerably lowers its condition number. A good recipe is: cj2 = k1+k2r( xj) , where r( xj) is an estimate of the radius of curvature of the PDE solution, and k1 and k2 are adjustable parameters that control the magnitude of the shape parameters.
  6. Both MQ-RBFs and FEM can yield very ill-conditioned systems of linear equations. It is common practice to pre-condition the matrices as well as use domain decomposition or matrix partitioning methods that not only render parallelization readily but deal with smaller rank matrices that are better conditioned. Wong et al. (1999) found that a multi-zone decomposition of the complex geometry of Tolu Harbour in Hong Kong yielded very reliable results. Kansa and Hon (1999) found that domain decomposition greatly reduced the condition number of the asymmetric PDE collocation matrix and, consequently, the root mean squared errors were reduced by a couple of orders of magnitude.
  7. The truncated MQ RBF developed by Kansa and Hon (1999) would be recommended for large scale problems. The hybrid combination of truncated-MQ + pre-conditioning+ domain decomposition+ alternating Schwartz method + iterative refinement and/or iterative methods give the techniques that are useful in solving large scale problems. The advantage of using exponentially convergent meshless basis functions is that orders of magnitude fewer discretization nodes are required. Since the FEM community has developed a considerable amount of methods and software, the RBF community should incorporate some of the appropriate methodology and software tools for the RBF collocation method.

10, Recommendations for Future Development.

Galperin and Zheng (1993) have pointed out that collocation methods are inherently ill-conditioned. Ill-posed and badly formulated problems can possess h-equivalent solutions that represent physical reality despite the mathematical nonexistence of an exact solution. In the previous references cited of Kansa and Hon, they performed an outer loop optimization of the shape-parameters and nodes using the MQ asymmetric collocation method. Only Galperin et al. have used global optimization on a few limited problems, with extra-ordinary results.

Although it is clear that the numerical solutions of PDE, ODE, integral, and integro-differential equations would greatly benefit from the global optimization, Galerkin-like formulation, the major implementation impediment is the lack of robust multi-parameter global optimization software. Unfortunately, gradient based methods are ill-conditioned, and converge rapidly only under certain restricted conditions. In addition, gradient methods pose the risk of being trapped in a local minimum, rather than in the global minimum. Ferrari and Galperin (1993) have published a software package of a fast one-dimensional adaptive cubic algorithm. It is hopeful that fast multi-dimensional global optimization software packages would be developed soon.


11, References.

Buhmann, M.D., Multivariate interpolation in odd-dimensional Euclidean spaces using multiquadrics; Const. Approx; 6; 21-34 (1990).

Buhmann, M.D. and C.A. Micchelli, ''Spline prewavelets for nonuniform nodes'', Numer. Math., 61: 455-474 (1992).

Chui, C.K., J. Stoeckler, and J.D. Ward, ''Analytic wavelets generated by radial functions'', Adv. Comput. Math. 5: 95-123 (1996).

Fasshauer, G. and J. Jerome, ''Multistep approximations algorithms: improved convergence rates through postconditioning with smoothing kernals'', preprint (1998).

A.I. Fedoseyev, M.J. Friedman, and E. J. Kansa. ''Continuation for nonlinear elliptic partial differential equations discretized by the multiquadric method'', to appear in Int. J. Bif. and Chaos.

Ferrari, A. and E.A. Galperin,''Numerical experiments with one-dimensional adaptive cubic algorithm'', Comput. Math. Applic. 25 (10/11): 47-56 (1993).

Franke,C. and. R. Schaback, ''Solving partial differential equations by collocation using radial basis functions'', Appl. Math. Comput. 93, 73-83 (1998),

Franke, R., ''Scattered data interpolation: tests of some methods'', Math. Comput. 38: 181-200 (1982)

Galperin,E.A. and Q. Zheng, ''Solution and control of PDE via global optimization methods'', Comput. Math. Applic. 25 (10/11): 103-118 (1993).

Galperin,E.A., Z. Pan and Q. Zheng, ''Application of global optimization to implicit solution of partial differential equations'', Comput. Math. Applic. 25 (10/11): 119-124 (1993).

Galperin, E.A. E.J. Kansa, A. Makroglou, and S.A. Nelson, '' Numerical Solutions of weakly singular Volterra equations using global optimization and radial basis functions'', In preparation.

Golberg, M.A. and C.S. Chen, ''Discrete projection methods for integral equations''. Computational Mechanics Publications, Southampton UK and Boston, MA (1997).

Hardy, R.L. ''Multiquadric equations of topography and other irregular surfaces'', J. Geophysics Res. 176: 1905-1915 (1971).

Hardy, R.L. '' Theory and applications of the multiquadric-biharmonic method: 20 years of discovery'', Comput. Math. Applic. 19 (8/9): 163-208 (1990).

Hon, Y.C. and X.Z. Mao,'' An efficient numerical scheme for Burgers equation'', Appl. Math. Comput., 95: 37-50 (1998).

Kansa, E.J. ''Multiquadrics- A scattered data approximation scheme with applications to computational fluid dynamics: I. Surface approximations and partial derivative estimates'', Comput. Math. Appl., 19, (6-8): 127-145 (1990).

Kansa, E.J. ''Multiquadrics- A scattered data approximation scheme with applications to computational fluid dynamics: II. Solutions to parabolic, hyperbolic, and elliptic partial differential equations'', Comput. Math. Appl., 19, (6-8): 147-161 (1990).

Kansa, E. J. and Y.C. Hon, ''Circumventing the ill-conditioning problem with multiquadric radial basis functions: Applications to elliptic partial differential equations'' submitted to Comp. Math. Applic.

Madych,W.R. and S.A. Nelson, Multivariate interpolation and conditionally positive definite functions, II'', Math. Comput. 54, 211-230 (1990).

Micchelli, C.A., ''Interpolation of scattered data: Distance matrices and conditionally positive definite functions'', Constr. Approx. 2: 11-22 (1986).

Moridis, G. J., and E. J. Kansa, ''The Laplace Transform Multiquadrics method: A highly accurate scheme for the numerical solution of linear PDEs'', J. Appl. Sci. & Comp. 1, (2): 375-407 (1994).

Schaback, R. and Y.C. Hon, ''On unsymmetric collocation with radial basis functions'', preprint (1999).

Wendland, H. ''Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree'', Adv. Comput. Math. 4: 389-396 (1995).

Wong, S.M., Y.C. Hon, T.S. Li, S.L. Chung, and E.J. Kansa, ''Multi-zone decomposition of time-dependent problems using the mulitquadric scheme'', Comput. Math. Applic.. 37: 23-43 (1999).

Wu, Zongmin , ''Solving PDE with radial basis functions and the error estimation'', Adv. in Comput. Math., Lectures in Pure and Applied Mathematics, 202, eds. Z. Chen, Y.Li, C.A. Micchelli, Y. Xu, and M Dekkon (1998).


12, Multiquadric Publications of E.J. Kansa.

1. E. J. Kansa, ''Application of Hardys multiquadric interpolation to hydrodynamics'', Simulations 4: 111-117 (1986).

2. E. J. Kansa, ''Multiquadrics- A scattered data approximation scheme with applications to computational fluid dynamics: I. Surface approximations and partial derivative estimates'', Comput. Math. Appl., 19, (6-8): 127-145 (1990).

3. E. J. Kansa, ''Multiquadrics- A scattered data approximation scheme with applications to computational fluid dynamics: II. Solutions to parabolic, hyperbolic, and elliptic partial differential equations'', Comput. Math. Appl., 19, (6-8): 147-161 (1990).

4. E. J. Kansa, ''A strictly conservative spatial approximation scheme for the governing engineering and physics equations over irregular regions and inhomogeneously scattered nodes'', Comput. Math. Appl. 24 (6): 169-190 (1992).

5. E. J. Kansa and R. E. Carlson, ''Improved accuracy of multiquadric interpolation using variable shape parameters'', Comput. Math. Appl. 24, (12): 99-120 (1992).

6.G. J. Moridis and E. J. Kansa, ''The Laplace Transform Multiquadrics method: A highly accurate scheme for the numerical solution of linear PDEs'', J. Appl. Sci. & Comp. 1, (2): 375-407 (1994).

7. R. E. Hagan and E. J. Kansa, ''Studies of the R parameter in the multiquadric function applied to groundwater pumping'', J. Appl. Sci. & Comp. 1, (2): 266-282 (1994).

8. E. J. Kansa and R. E. Carlson, ''Radial basis functions: A class of grid-free scattered data approximations'', J. Comput. Fluid. Dynamics. 3, (4), 489-496 (1995).

9. E.A. Galperin, E.J. Kansa, A. Makroglou, and S.A. Nelson, '' Numerical Solutions of weakly singular Volterra equations using global optimization and radial basis functions'', Proc. Hermas96, Ed. E.A. Lipitakis, Hellenic European Research on Mathematics and Informatics 96, LEA, Athens. 3: 550-559 (1997).

10. Y.C. Hon, K.F. Cheung, X. Z Mao, and E. J. Kansa, ''A multiquadric solution for the shallow water equations'', Journal of Hydraulogy 125 (5): 524-533 (1999).

11. S. M. Wong, Y.C.Hon, T.S.Li, S.L.Chung, and E.J.Kansa, ''Multi-zone decomposition for simulation of time-dependent problems using the multiquadric scheme'', Comput. Math. Applic. 37: 23-43 (1999).

12. E.J. Kansa and Y.C. Hon, ''Circumventing the ill-conditioning problem with multiquadric radial basis functions: Applications to elliptic partial differential equations'', submitted to Comput. Math. Applic.

13. A.I. Fedoseyev, M.J. Friedman, and E. J. Kansa. ''Continuation for nonlinear elliptic partial differential equations discretized by the multiquadric method'', to appear in Int. J. Bif. and Chaos.



Guest Editor
Computers and Mathematics with Applications 24 (12) (1992).
J. Applied Science and Computing, 1 (1), (2), and (3) (1994).

File translated from TEX by TTH, version 2.00. On 24 Aug 1999, 17:42.