Solving Partial Differential Equations by Collocation with Radial Basis Functions

Gregory E. Fasshauer

Abstract. Motivated by [5] we describe a method related to scattered Hermite interpolation for which the solution of elliptic partial differential equations by collocation is well-posed. We compare the method of [5] with our method.

§1. Introduction

In this paper we discuss the numerical solution of elliptic partial differential equations using a collocation approach based on radial basis functions. To make the discussion transparent we will focus on the case of a time independent linear elliptic partial differential equation in IR^{2}.

In the following we assume we are given a set of nodes Ξ = {ξ⃗_{1}, …, ξ⃗_{N}} ⊂ IR^{d}, along with a continuous function φ : [0, ∞) → IR. We then refer to x⃗ ↦ φ(∥x⃗−ξ⃗_{k}∥_{2}), x⃗ ∈ IR^{d}, k ∈ {1, …, N}, as radial basis functions centered at ξ⃗_{k}. Some of the most commonly used radial basis functions are the (reciprocal) multiquadrics φ(r) = (r^{2} + c^{2})^{±1/2}, the Gaussians φ(r) = e^{−cr^{2}}, and the thin plate splines φ(r) = r^{2} log r in IR^{2}. Here we let r = ∥x⃗ − ξ⃗_{k}∥, and c > 0 is a parameter.

Some of the advantages of radial basis functions are their insensitivity to the spatial dimension d, which makes the implementation of this method in higher dimensions much simpler than, e.g., finite elements. Another useful feature of radial basis functions is their radial symmetry and invariance under Euclidean transformations. Furthermore, in the context of scattered data interpolation it is known that some radial basis functions have spectral convergence orders ( e.g., (reciprocal) multiquadrics, Gaussians). This should also be evident in some form when using them for collocation. There are also some well-known disadvantages. Radial basis functions are generally globally supported and poorly conditioned. One usually tries to remedy these problems by localization, preconditioning, and fine tuning of the parameter c.

Proceedings of Chamonix 1996 1 A. Le Méhauté, C. Rabut, and L. L. Schumaker (eds.), pp. 1–8. Copyright o^{c} 1997 by Vanderbilt University Press, Nashville, TN. ISBN 1-xxxxx-xxx-x. All rights of reproduction in any form reserved.

We will not concentrate on any of these issues below, our main point being to establish a basic framework in which a collocation method similar to that proposed by Kansa in [5] is well-posed. An indication of the success of Kansa’s method (which has not yet been shown to be well-posed) are [1, 2, 3, 6, and 9].

§3. Kansa’s Approach

In order to be able to clearly point out the differences between Kansa’s method and our approach we start with a review of scattered data interpolation with radial basis functions in IR^{d}.

In this context we are given data {x⃗_{j}, f_{j}}, j = 1, …, N, x⃗_{j} ∈ IR^{d}, where we can think of the values f_{j} being sampled from a function f : IR^{d} → IR. The goal is to find an interpolant of the form

s(x⃗) = ∑_{k=1}^{N} c_{k}φ(∥x⃗ − x⃗_{k}∥),  x⃗ ∈ IR^{d},  (1)

such that

s(x⃗_{j}) = f_{j}, j = 1, …, N.

The solution of this problem leads to a linear system Ac⃗ = f⃗ with the entries of A given by

A_{jk} = φ(∥x⃗_{j} − x⃗_{k}∥),  j, k = 1, …, N.  (2)

Clearly, there exists a unique solution if and only if A is non-singular. The challenge now is to find the largest possible class of functions φ for which this is true. This is still an unsolved problem. However, e.g., in [7, 8] the non-singularity of A for certain classes of φ was shown by using either a variational approach, or making use of the fact that φ is either strictly positive definite on IR^{d} (reciprocal multiquadrics, Gaussians), or strictly conditionally positive definite of order one (multiquadrics). If one wants to use thin plate splines, one has to append at least a linear polynomial to (1) and enforce some well-known additional conditions. The fact that in this case φ is strictly conditionally positive definite of order two guarantees the non-singularity of the corresponding system matrix.

We are now ready to discuss Kansa’s collocation method. Assume we are given a domain Ω ⊂ IR^{d}, and a linear elliptic partial differential equation of the form

L[u](x⃗) = f(x⃗),  x⃗ in Ω,  (3)

with (for simplicity of description) Dirichlet boundary conditions

u(x⃗) = g(x⃗),  x⃗ on ∂Ω.  (4)

We then choose to represent u by radial basis functions, i.e.,

u(x⃗) = ∑_{k=1}^{N} c_{k}φ(∥x⃗ − ξ⃗_{k}∥),  (5)

where the points ξ⃗_{1}, …, ξ⃗_{N} are a set of centers for the radial basis functions which usually are selected to coincide with the collocation points X = {x⃗_{1}, …, x⃗_{N}} ⊂ Ω. We assume the simplest possible setting here, i.e., no polynomial terms are added to the expansion (5). The collocation matrix which arises when matching the differential equation (3) and the boundary conditions (4) at the collocation points X will be of the form

A = ( \begin{matrix}Φ \\ L[Φ] \end{matrix} ) ,  (6)

where the two blocks are generated as follows:

\begin{aligned}Φ_{jk} & = φ(∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j} ∈ B, ξ⃗_{k} ∈ X , \\ L[Φ]_{jk} & = L[φ](∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j} ∈ I, ξ⃗_{k} ∈ X . \end{aligned}

Here we have identified (as we will do throughout the paper) the set of centers with the set of collocation points. The set X is split into a set I of interior points, and B of boundary points. The problem is well-posed if the linear system Ac⃗ = y⃗, with y⃗ a vector consisting of entries g(x⃗_{j}), x⃗_{j} ∈ B, followed by f(x⃗_{j}), x⃗_{j} ∈ I, has a unique solution.

Remark. We note that a change in the boundary conditions (4) is as simple as changing a few rows in the matrix A in (6) as well as on the right-hand side y⃗.

Kansa [5] claims the matrix A in (6) is of the same type as that in (2), and therefore refers to [8] to conclude that his method is well-posed. This is not true, however, since the matrix for the collocation problem is composed of rows which are built from different functions (which – depending on the differential operator L – might not even be radial). The results in [8] on the other hand are based on the fact that A is generated by a single function φ. We point out that in [5] the author describes three sets of experiments using his method and comments on the superior performance of multiquadrics in terms of computational complexity and accuracy when compared to finite difference methods. An interesting open question therefore is whether the well-posedness of Kansa’s method can be established.

Before we describe an alternate approach based on Hermite basis functions which ensures well-posedness we would like to point out that in [9] the authors suggest how Kansa’s method can be applied to other types of partial differential equation problems such as non-linear elliptic PDEs, systems of elliptic PDEs, and time-dependent parabolic or hyperbolic PDEs.

§4. An Hermite-based Approach

Our approach is based on scattered Hermite interpolation (see e.g., [4, 10, 12, 13]), which we now also quickly review. In this context we are given data {x⃗_{j}, L_{j}f}, j = 1, …, N, x⃗_{j} ∈ IR^{d} where L = {L_{1}, …, L_{N}} is a linearly independent set of continuous linear functionals. We try to find an interpolant of the form

s(x⃗) = ∑_{k=1}^{N} c_{k}Lξ⃗ ; kφ(∥x⃗ − ξ⃗∥),  x⃗ ∈ IR^{d},  (7)

satisfying

L_{j}s = L_{j}f, j = 1, …, N.

We have used Lξ⃗ to indicate that the functional L acts on φ viewed as a function of the second argument ξ⃗. The linear system Ac⃗ = Lf which arises in this case has matrix entries

A_{jk} = L_{j}Lξ⃗ ; kφ,  j, k = 1, …, N.  (8)

In the references mentioned at the beginning of this section it is shown that A is non-singular for the same classes of φ as given in [8].

We illustrate the situation with a simple example.

Example. Let data \{x⃗_{j}, f(x⃗_{j})\}_{j=1}^{n} and \{x⃗_{j}, \frac{∂f}{∂x}(x⃗_{j})\}_{j=n+1}^{N} with x⃗ = (x, y) ∈ IR^{2} be given. Then

s(x⃗) = ∑_{k=1}^{n} c_{k}φ(∥x⃗ − x⃗_{k}∥) − ∑_{k=n+1}^{N} c_{k} \frac{∂φ}{∂x}(∥x⃗ − x⃗_{k}∥),

and

A = [ \begin{matrix}Φ & −Φ_{x} \\ Φ_{x} & −Φ_{xx} \end{matrix} ] ,

with

\begin{aligned}Φ_{jk} & = φ(∥x⃗_{j} − x⃗_{k}∥),  j, k = 1, …, n, \\ −Φ_{x,jk} & = −\frac{∂φ}{∂x}(∥x⃗_{j} − x⃗_{k}∥),  j = 1, …, n, k = n + 1, …, N, \\ Φ_{x,jk} & = \frac{∂φ}{∂x}(∥x⃗_{j} − x⃗_{k}∥),  j = n + 1, …, N, k = 1, …, n, \\ Φ_{xx,jk} & = \frac{∂^{2}φ}{∂x^{2}}(∥x⃗_{j} − x⃗_{k}∥),  j, k = n + 1, …, N. \end{aligned}

Now we describe an alternative collocation method based on the theory just reviewed. Assume we are given the same PDE (3) with boundary conditions (4) as in the description of Kansa’s method. In order to be able to apply the results from scattered Hermite interpolation to ensure the non-singularity of the collocation matrix we propose the following expansion for the unknown function u:

u(x⃗) = ∑_{k=1}^{#B} c_{k}φ(∥x⃗ − ξ⃗_{k}∥) + ∑_{k=#B+1}^{N} c_{k}Lξ⃗[φ](∥x⃗ − ξ⃗_{k}∥),  (9)

where #B denotes the number of nodes on the boundary of Ω, and Lξ⃗ is the differential operator used in (3), but acting on φ viewed as a function of the second argument, i.e., L[φ] is equal to Lξ⃗[φ] up to a possible difference in sign. Note the difference in notation. In (7) L is a linear functional, and in (9) a differential operator.

This expansion for u leads to a collocation matrix A which is of the form

A = ( \begin{matrix}Φ & Lξ⃗[Φ] \\ L[Φ] & L[Lξ⃗[Φ]] \end{matrix} ) ,  (10)

where the four blocks are generated as follows:

\begin{aligned}Φ_{jk} & = φ(∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j}, ξ⃗_{k} ∈ B, \\ Lξ⃗[Φ]_{jk} & = Lξ⃗[φ](∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j}, ∈ B, ξ⃗_{k} ∈ I, \\ L[Φ]_{jk} & = L[φ](∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j} ∈ I, ξ⃗_{k} ∈ B, \\ L[Lξ⃗[Φ]]_{jk} & = L[Lξ⃗[φ]](∥x⃗_{j} − ξ⃗_{k}∥),  x⃗_{j}, ξ⃗_{k} ∈ I. \end{aligned}

The matrix (10) is of the same type as the scattered Hermite interpolation matrices (8), and therefore non-singular as long as φ is chosen appropriately. Thus, viewed using the new expansion (9) for u, the collocation approach is certainly well-posed. Another point in favor of the Hermite based approach is that the matrix (10) is (anti-)symmetric as opposed to the completely unstructured matrix (6) of the same size. This property should be of value when trying to devise an efficient implementation of the collocation method. Also note that although A consists of four blocks now, it still is of the same size, namely N × N, as the collocation matrix (6) obtained for Kansa’s approach.

§5. Numerical Examples

In the following test examples we restrict ourselves to two-dimensional Poisson problems whose analytic solution is readily available and therefore can easily be verified. We will refer to a point in IR^{2} as (x, y). In all of the following tests we used multiquadrics in the expansions (5) and (9) of the unknown function u.

Example 1. Consider the Poisson equation

Δu(x, y) = y(1 − y) sin^{3} x, x ∈ (0, π), y ∈ (0, 1),

with Dirichlet boundary conditions

u(x, 0) = u(x, 1) = u(0, y) = u(π, y) = 0.

For this test problem we selected various uniform grids as listed in Tables 1 and 2 on [0, π] × [0, 1]. Tables 1 and 2 show the values of the multiquadric parameter, the relative maximum errors ρ computed on a fine grid of 60 × 60

c  ρ_{K}  ρ_{H}  cond_{K}(A) cond_{H}(A)
5 × 3  1.0 5.248447e-02 2.004420e-01 2.599606e+03 1.627432e+03
8 × 4  1.0 1.126843e-02 1.124710e-02 2.325758e+05 8.167527e+04
10 × 6 1.0 5.809472e-03 6.481697e-03 4.321740e+07 1.808001e+07
16 × 8 1.0 1.347863e-03 1.720007e-03 8.685785e+10 1.496772e+10
20 × 12 1.0 5.053090e-04 5.973294e-04 5.161540e+15 1.234633e+15

Table1. Error progression for increasingly denser data sets (Ex.1, fixed c).

c_{K} c_{H}  ρ_{K}  ρ_{H}  cond_{K}(A) cond_{H}(A)
5 × 3  1.18 1.39 1.627193e-02 4.180428e-02 5.592238e+03 5.231279e+03
8 × 4  1.04 1.11 1.103747e-02 1.062891e-02 3.175078e+05 1.735482e+05
10 × 6 4.80 3.84 2.739293e-03 3.451799e-03 1.193586e+18 1.414927e+15
16 × 8 3.12 3.12 2.707006e-04 2.082886e-04 1.209487e+19 6.609375e+18
20 × 12 2.00 2.30 3.894511e-05 1.273363e-05 3.739554e+19 6.750955e+18

Table 2. Error progression for increasingly denser data sets (Ex.1, “optimal” c).

points, and the approximate condition numbers of A. The range of u on the evaluation grid is approximately [−0.021023, 0.0]. The “optimal” value for c was determined by trial and error. The subscripts K and H refer to Kansa’s and the Hermite based method, respectively.

Fig. 1 shows the distribution of the errors |u(x) − s(x)| on the evaluation grid for the two methods on the 8 × 4 grid used in Table 2. The scale used for the shading is displayed on the right.

Fig. 1. Error for Kansa’s (top), Hermite (bottom) solution for Ex. 1 on 8 × 4 grid.

Example 2. Consider the Poisson equation

Δu(x, y) = sin x − sin^{3} x, x ∈ (0, \frac{π}{2} ), y ∈ (0, 2),

with mixed Dirichlet and Neumann boundary conditions

u(0, y) = u_{x}(\frac{π}{2} , y) = u_{y}(x, 0) = u_{y}(x, 2) = 0.

For this example we selected uniform grids on [0, π/2] × [0, 2] as listed in Table 3. This time we only list the results for the “optimal” choice of c. The values listed are analogous to those in Ex. 1.

c_{K} c_{H}  ρ_{K}  ρ_{H}  cond_{K}(A) cond_{H}(A)
3 × 3  109.0 2.19 9.628085e-01 1.141043e-01 1.592286e+16 5.560886e+02
5 × 5  1.80 1.73 2.181029e-02 4.327029e-02 2.395293e+06 1.271196e+05
7 × 7  1.58 3.56 6.910084e-03 1.871798e-04 5.762316e+08 1.854850e+12
10 × 10 2.80 3.29 9.265197e-05 5.126676e-05 2.842111e+18 7.070804e+17
14 × 14 2.28 2.62 1.138751e-05 1.725526e-06 6.573143e+19 5.891454e+18
20 × 20 1.53 1.91 5.501057e-06 6.217559e-07 5.889491e+19 7.576112e+19

Table 3. Error progression for increasingly denser data sets (Ex.2, “optimal” c).

Discussion: All in all the Hermite method seems to perform slightly better than Kansa’s method. Especially for the cases in which we used relatively many interior points (which is where the method differ). Also, the matrices for the Hermite method generally have smaller condition numbers. An advantage of the Hermite approach over Kansa’s method is that for the differential operator L used here, the collocation matrices resulting from the Hermite approach are symmetric. Therefore the amount of computation can be reduced considerably, which is important for larger problems. Kansa’s method has the advantage of being simpler to implement (since less derivatives of the basis functions are required).

Remark. In order to improve the performance as well as the accuracy one could use the concept of blending together local approximations to form a smooth global approximation coupled with domain decomposition. This would lead to a drastic reduction in the amount of computation when dealing with large problems, and also lends itself to parallelization. Kansa [6] reports success with this method. Another possible venue of attack is the use of locally supported radial basis functions (see e.g., [11]).

References

1. Dubal, M. R., Construction of three-dimensional black-hole initial data via multiquadrics, Phys. Rev. D 45 (1992), 1178–1187.

2. Dubal, M. R., Domain decomposition and local refinement for multiquadric approximations. I: Second-order equations in one-dimension, J. Appl. Sci. Comp. 1 (1994), 146–171.

3. Golberg, M. A., C. S. Chen, and S. R. Karur, Improved multiquadric approximation for partial differential equations, University of Nevada, 1995.

4. Iske, A., Reconstruction of functions from generalized Hermite-Birkhoff data, in Approximation Theory VIII, Vol. 1: Approximation and Interpolation, C. Chui, and L. L. Schumaker (eds), World Scientific Publishing, Singapore, 1995, 257–264.

5. 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 (1990), 147–161.

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

7. Madych, W. R. and S. A. Nelson, Multivariate interpolation: a variational theory, manuscript, 1983.

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

9. Moridis, G. J. and E. J. Kansa, The Laplace transform multiquadric method: A highly accurate scheme for the numerical solution of linear partial differential equations, J. Appl. Sci. Comp. 1 (1994), 375–407.

10. Narcowich, F. J. and J. D. Ward, Generalized Hermite interpolation via matrix-valued conditionally positive definite functions, Math. Comp. 63 (1994), 661–687.

11. Schaback, R., Creating surfaces from scattered data using radial basis functions, in Mathematical Methods for Curves and Surfaces, M. Dæhlen, T. Lyche, L. L. Schumaker (eds), Vanderbilt University Press, Nashville, 1995, 477–496.

12. Sun, X., Scattered Hermite interpolation using radial basis functions, Linear Algebra Appl. 207 (1994), 135–146.

13. Wu, Z., Hermite-Birkhoff interpolation of scattered data by radial basis functions, Approx. Theory Appl. 8 (1992), 1–10.

Gregory E. Fasshauer Department of Mathematics Northwestern University Evanston, IL 60208

fass@math.nwu.edu