On the numerical solution of differential equations with radial basis functions

Gregory E. Fasshauer

Department of Applied Mathematics, Illinois Institute of Technology, Chicago, IL 60616, U.S.A. E-Mail: fass@amadeus.csam.iit.edu

Abstract

In this paper we report on two different experiments dealing with the numerical solution of differential equations by radial basis functions: 1) the solution of a two-point boundary value problem; 2) the solution of a twodimensional Poisson equation. In the second experiment we contrast a multilevel collocation algorithm based on locally supported basis functions with two different direct solution approaches (one based on locally supported basis functions, the other on globally supported multiquadrics). In both experiments the effects of a smoothing operation are studied.

1. Introduction

Since their introduction in the 1970s radial basis functions (RBFs) have attracted a lot of attention from mathematicians working in the area of approximation theory. More recently, however, a large number of scientists and engineers have discovered their merits when working with multidimensional data. Furthermore, the main focus of the applications seems to have shifted from scattered data approximation to the numerical solution of partial differential equations (PDEs). It is this second class of problems we will address in this paper.

The RBFs being used fall into two general categories: globally supported, and locally supported functions. Globally supported RBFs (i.e. multiquadrics, Gaussians, thin plate splines, etc.) were the ones used early on, and for problems whose size does not exceed 400-500 pieces of data these methods should still be the method of choice. For larger problems the exact solution of the associated (dense) linear systems becomes too expensive. Therefore, researchers interested in continuing to work with globally supported RBFs have turned to the construction of fast (approximate) solvers and evaluation algorithms. Most work in this area is being done by the group in Cambridge, England (see e.g. Powell[11]), as well as by Rick Beatson (see e.g. Beatson & Light[1]).

It is the class of locally supported functions, however, which we will focus our attention on in this paper. Those compactly supported RBFs used most widely were introduced by Wendland[13]. Simple direct formulas for computing these functions are (see Fasshauer[3])

\begin{aligned}ϕ_{ℓ,0}(r) & ≐ (1 − r)_{+}^{ℓ}, \\ ϕ_{ℓ,1}(r) & ≐ (1 − r)_{+}^{ℓ+1} [(ℓ + 1)r + 1] , \\ ϕ_{ℓ,2}(r) & ≐ (1 − r)_{+}^{ℓ+2} [(ℓ^{2} + 4ℓ + 3)r^{2} + (3ℓ + 6)r + 3] , \\ ϕ_{ℓ,3}(r) & ≐ (1 − r)_{+}^{ℓ+3} [(ℓ^{3} + 9ℓ^{2} + 23ℓ + 15)r^{3} \\ & +(6ℓ^{2} + 36ℓ + 45)r^{2} + (15ℓ + 45)r + 15] . \end{aligned}  (1)

Here ≐ denotes equality up to a constant factor, and r = ∥x∥, where x ∈ IR^{d}, so that the composition of the univariate function ϕ_{ℓ,k} with the norm indeed yields a radial function. Note that the cut-off function (·)_{+} ensures the compact support of these functions. Moreover, it should be pointed out that the ϕ_{ℓ,k} are in fact piecewise polynomials. The two indices ℓ and k are related to the space dimension d, and to the smoothness of the basis functions. More precisely, if one intends to work in IR^{d}, then one should take ℓ = ⌊\frac{d}{2}⌋ + k + 1, where k is the value of the second index stating that ϕ_{ℓ,k} ∈ C^{2k}, since the functions are positive definite for this choice, and thus the associated system matrix is assured to be invertible.

The solution of partial differential equations with radial basis functions is a subject which is still in its infant stages, even though some first attempts were reported roughly 10 years ago (see Kansa[10]). The examples we are about to discuss are linear elliptic problems which we solve numerically by a symmetric collocation approach. To this end we assume that the PDE is of the form

Lu = F,

where L = [L_{1}, …, L_{N}]^{T} is a vector of linear operators, and the righthand side F = [f_{1}, …, f_{N}]^{T}. Together they define the (system of) partial differential equations along with boundary conditions. We also assume we are given a set of centers X = {x_{1}, …, x_{n}} at which we attempt to satisfy the PDE pointwise (collocation). The expansion one uses for the solution is of the form

u(x) = ∑_{ν=1}^{N} ∑_{j=1}^{#X_{ν}} c_{νj}L_{ν}^{(2)}ϕ(∥x − x_{j}∥_{2}),  x ∈ IR^{d} .  (2)

Here the X_{ν} are subsets of X with cardinality #X_{ν} corresponding to the components of L. The superscript^{(2)} indicates that L acts on ϕ as a function of the second variable, i.e., the center location. Expansion (2) is motivated by a connection to scattered Hermite interpolation which ensures the invertibility of the related collocation matrix (for more details see Fasshauer[2,4] or Franke & Schaback[7]).

The remainder of the paper is organized as follows. In Section 2 we outline our algorithmic approach which is a Newton-like multilevel method with interlaced smoothing operations, and then discuss our two examples in the following two sections. The last section is devoted to summarizing remarks.

2. Multilevel Collocation with Smoothing

It is not recommended that one use compactly supported functions in the same way as one uses globally supported functions, i.e., in a direct, singlegrid approach. The advantage that compactly supported functions add in terms of computational efficiency is offset by their modest approximation properties (see Example 2 in Section 4). To obtain the best results with compactly supported functions (for large data sets) it is recommended that a multilevel algorithm be employed (a more detailed discussion of this issue can be found in Fasshauer[4] or in Schaback[12]).

The basic multilevel algorithm for the solution of the problem Lu = F on Ω ⊂ IR^{d} by collocation can be described as follows:

Multilevel Collocation Algorithm.

1. Generate a nested sequence of computational grids X_{1} ⊂ ⋯ ⊂ X_{K} ⊂ Ω, and let u_{0} = 0.

2. For k from 1 to K do a. Compute the residual r_{k} = F − Lu_{k−1} on the grid X_{k}. b. Solve Lu_{k} = r_{k} on X_{k} with u_{k} as in Eqn. (2). c. Update the level-k approximation as

u_{k} ← u_{k−1} + u_{k}.

By linearity and a simple telescoping argument it is clear that u_{K} will satisfy the differential equation at all points of the finest grid X_{K}.

In Fasshauer & Jerome[5] it was shown how a multilevel algorithm of this kind can be interpreted as Newton iteration at the operator level. This implies that step 2 above is viewed as

u_{k} − u_{k−1} = −T_{X_{k}}F (u_{k−1}).  (3)

Here the operator F computes the residual, i.e., F (u) = Lu − F, and the operator T_{X} (which encompasses the solution of the differential equation via collocation in our examples) can be interpreted as an approximation to the inverse of the (Frechèt) derivative of F.

It is well known (see e.g. Jerome[9]) that – unlike the classical Newton method – the operator version does in general not converge quadratically, but suffers from a so-called loss of derivatives. The main result of Jerome[9], which was adapted for Besov spaces in Fasshauer & Jerome[5], is that the addition of an appropriate smoothing operation S_{k} (defined via convolution with an appropriate convolution kernel ψ) at each step of the Newton iteration almost restores quadratic convergence. In other words, the iterative scheme

u_{k} − u_{k−1} = −S_{k}T_{X_{k}}F (u_{k−1})  (4)

converges superlinearly. The precise statement of the convergence of the algorithm (in a variety of Besov spaces) is given in Theorem 3.1 of Fasshauer & Jerome[5]. The two most important hypotheses required for the superlinear convergence result to hold are:

• The smoothing operation needs to be an approximate identity satisfying a Bernstein and a Jackson inequality in a very general scale of Besov spaces (assumption (A2) in Fasshauer & Jerome[5]).

• The meshsize of each computational grid has to be chosen adaptively according to the size of the most recent residual ((3.4) in Fasshauer & Jerome[5]).

Thus, the addition of an appropriate smoothing operation in step 2.c of the multilevel collocation algorithm will generate a (superlinearly) convergent scheme, provided that the sequence of computational grids satisfies the meshsize constraint mentioned above. Another way to interpret the results of the smoothing operation is that the smoothing restores the convergence properties of the numerical inversion operator T_{X} (see (3.2) and (3.6) of Fasshauer & Jerome[5]).

3. Smoothing with Spline Kernels

The main properties required in the proof of the Jackson inequality for the smoothing operation in Fasshauer & Jerome[5] were:

• ψ̂ ∈ C_{0}^{∞}(IR^{d}),

• ψ̂ ≡ 0 in a neighborhood B_{ϵ}(0) of the origin.

Here ψ̂ denotes the d-variate Fourier transform of ψ. Smoothing kernels satisfying all of these properties are not very easily found. One example, however, was listed already in Jerome[9]. Since that kernel is not easy to implement, other possibilities (which satisfy only some of the theoretical requirements) were suggested in Fasshauer & Jerome[5]. One such possibility are kernels whose Fourier transforms are compactly supported on IR^{d}, and which also satisfies the condition ψ̂ ≡ 1 in B_{ϵ}(0). However, the regularity of the Fourier transform is assumed to be limited. The simplest kernel of this sort is the sinc function

ψ(x) = \frac{sin(x)}{πx},

whose Fourier transform is the characteristic function of the interval [−1, 1]. For general d, the inverse Fourier transform of the characteristic function of the unit ball in IR^{d} is given by

ψ(x) = \frac{J_{d/2}(∥x∥)}{(2π∥x∥)^{d/2}} ,

where J_{ν} is the Bessel function of the first kind of order ν.

Reverting to the case d = 1, kernels whose Fourier transform is smoother can be constructed using Hermite interpolation conditions at ±1 and ±ϵ. Therefore the kernel whose Fourier transform is given by the continuous function

ψ̂(ω) = \{ \begin{matrix}1 & |ω| ≤ 1/2, \\ 2(1 − ω) & 1/2 < |ω| ≤ 1, \\ 0 & otherwise, \end{matrix}

is given by

ψ(x) = 2\frac{cos(\frac{x}{2}) − cos(x)}{πx^{2}},

or the C^{1} Fourier transform

ψ̂(ω) = \{ \begin{matrix}1 & |ω| ≤ 1/2, \\ 4(4ω^{3} − 9ω^{2} + 6ω − 1) & 1/2 < |ω| ≤ 1, \\ 0 & otherwise, \end{matrix}

leads to

ψ(x) = 24\frac{4 (cos(\frac{x}{2}) − cos(x)) − x (sin(x) + sin(\frac{x}{2}))}{πx^{4}}.  (5)

Of course it is also possible to construct kernels whose transform is smoother.

Example 1. Consider the two point boundary-value problem

\begin{matrix}−u′′(x) + π^{2}u(x) = & 2π^{2} sin πx,  x ∈ (0, 1), \\ u(0) = & u(1) = 0, \end{matrix}  (6)

with solution u(x) = sin πx. As computational grids X_{k} we take 2^{k+3} + 1 uniformly spaced points on [0, 1] as indicated in Table 1. Even though we are dealing with a one-dimensional problem we use the radial basis function ϕ_{5,3} as defined in Eqn. (1), which is six times continuously differentiable. Note that since the differential equation (6) is of second order, the RBF chosen needs to be at least C^{4} for use with the symmetric collocation approach (2). The linear systems arising at each level are solved with the conjugate gradient method and Jacobi preconditioning. In order to have a meaningful initial approximation the support of the basis functions for the coarsest grid is chosen so large that the matrix is a dense one. At subsequent levels the

mesh w/o smoothing  w/ smoothing  % nonzero
ℓ_{2}-error  rate  ℓ_{2}-error  rate
17 5.408189 10^{−6}  5.681042 10^{−6}  100
33 1.160280 10^{−7} 5.54 5.590327 10^{−8} 6.67  78.0
65 1.548268 10^{−8} 2.91 6.070907 10^{−9} 3.20  46.6
129 1.118532 10^{−8} 0.47 5.538264 10^{−9} 0.13  25.3
257 1.093402 10^{−8} 0.03 9.823599 10^{−9} −0.83  13.2
513 1.092675 10^{−8} 0.00 1.087111 10^{−8} −0.15  6.71
1025 1.093092 10^{−8} 0.00 1.261928 10^{−8} −0.22  3.39
2049 1.093351 10^{−8} 0.00 1.357822 10^{−8} −0.11  1.70

Table 1: Errors and convergence rates with and without C^{1}-spline smoothing for BVP (6).

support size for the new basis functions is taken to be half that of the one for the previous level. We do this in order to keep the bandwidth of the system matrices constant (note that the mesh size is also halved in each iteration). The rightmost column of Table 1 shows the percentage of nonzero entries in the system matrix at each level.

Columns 2 and 3 of Table 1 show how the multilevel collocation algorithm performs without the recommended smoothing. Clearly, it ceases to converge after the first 4 steps. The smoothing which is used for the results presented in columns 4 and 5 is based on the C^{1} spline kernel (5). The smoothing operation S_{k} in Eqn. (4) depends on the level k, in fact ψ_{t_{k}}(x) = t_{k}^{d}ψ(t_{k}x). In Fasshauer & Jerome[5] the smoothing speeds t_{k} were defined as t_{k} = ρ^{θβ^{k}}, where the three parameters ρ, θ and β can be chosen by the user (with some dependency on the smoothness of the problem). For the results of Table 1 we have used ρ = 10, θ = 1.1, and β = 1.2. The reader can see at least some benefits of the smoothing in this example since the errors (as well as the convergence rates) are better for the first few steps of smoothing. Later on, however, the performance of the algorithm with smoothing is even worse than without. This may be attributable to the lack of smoothness, or possibly to a stronger influence of the oscillations in the smoothing kernel on these finer grids.

4. A Boundary-Value Problem in 2D

We now give an example of a Poisson equation in two dimensions. The solution by multilevel collocation with compactly supported RBFs is, to our knowledge, the first such attempt. Besides discussing the influence of smoothing on the performance of the algorithm, we will also compare the multilevel results with direct, single-grid results based on multiquadrics and compactly supported RBFs.

Example 2. Consider the Poisson equation

\begin{aligned}∇^{2}u(x, y) = & − 30000e^{−1000((x− \frac{1}{2} )^{6}+(y− \frac{1}{2} )^{6})}× \\ & ((x − \frac{1}{2})^{4} − 1200(x − \frac{1}{2})^{10} + (y − \frac{1}{2})^{4} − 1200(y − \frac{1}{2})^{10}) , \end{aligned}

(7)

inside the unit square with Dirichlet boundary conditions

u(x, y) = e^{−1000((x− \frac{1}{2} )^{6}+(y− \frac{1}{2} )^{6})}, x = 0 or x = 1 or y = 0 or y = 1.

A plot of the solution is given on the left in Figure 1.

For this example the computational grids X_{k} are given by (2^{k+1} + 1) × (2^{k+1} + 1) uniformly spaced points in the unit square. Again, we use the RBF ϕ_{5,3}. The arrangement of the information in Table 2 is the same as in Table 1. However, the smoothing operation is now performed by convolving with the Gauss-Weierstrass kernel

ψ_{t}(x) = \frac{t^{2}}{(4π)}e^{− \frac{t^{2}∥x∥^{2}}{4}} , x ∈ IR^{2} .

The parameters determining the smoothing speeds are ρ = 1.2, θ = 1.5, and β = 1.9. It should be pointed out that the Gauss-Weierstrass kernel does not satisfy all the theoretical requirements for the kernel either. In particular, as a positive kernel it is saturated (i.e., the required Jackson inequality does not hold). However, the results below clearly show that its use improves the performance of the multilevel collocation algorithm. Plots of the two different fits on the 33 × 33 grid are shown in Figure 1.

mesh  w/o smoothing  w/ smoothing  % nonzero
ℓ_{2}-error  rate  ℓ_{2}-error  rate
5 × 5  1.656886 10^{0}  2.584879 10^{−1}  100
9 × 9 4.476316 10^{−1} 1.89 3.722518 10^{−2} 2.80  72.6
17 × 17 4.008961 10^{−1} 0.16 1.873797 10^{−2} 0.99  29.5
33 × 33 4.088999 10^{−1} −0.03 1.269354 10^{−2} 0.56  9.35
65 × 65 4.149843 10^{−1} −0.02 1.255188 10^{−2} 0.02  2.63

Table 2: Errors and convergence rates for multilevel collocation without and with Gauss-Weierstrass smoothing for BVP (7).

Finally, we would like to offer the following interesting comparison with collocation on a single grid. We do this with both the compactly supported RBF ϕ_{5,3} used above, and with globally supported multiquadrics, ϕ(r) = \sqrt{r^{2} + c^{2}}. The values of c used for the multiquadrics were chosen on an Figure 1: Plot of exact solution, and multilevel fits on 1089 points without and with smoothing for BVP (7).

ad-hoc basis, aided by a small amount of trail and error, and can therefore probably be improved upon. The NA entries for the multiquadrics in Table 3 indicate that the resulting dense linear systems (of size 1089 × 1089, and 4425×4225) were too large to solve on a desktop PC. The support size (and therefore the bandwidth of the matrices for the locally supported results were kept constant throughout resulting in the same percentage of nonzero entries as for the multilevel method. Even though the second entry for the compactly supported functions indicates a relatively good fit, the graph of the approximate solution in Figure 2 shows otherwise.

mesh  MQ  CS-RBF
ℓ_{2}-error  c  ℓ_{2}-error  % nonzero
5 × 5 7.571929 10^{−2} 0.2 1.656886 10^{0}  100
9 × 9 3.954529 10^{−3} 0.15 3.322141 10^{−2}  72.6
17 × 17 5.816741 10^{−4} 0.1 1.980783 10^{−1}  29.5
33 × 33  NA  NA 2.492955 10^{−1}  9.35
65 × 65  NA  NA 2.665762 10^{−1}  2.63

Table 3: Errors for single-grid collocation with multiquadrics and compactly supported RBFs for BVP (7).

The comparison with single-grid collocation shows that globally supported functions should be preferred on small data sets, and that compactly Figure 2: CS-RBF fit on 81 and 1089 points, and MQ fit on 289 points for BVP (7).

supported functions should probably be avoided altogether for direct computations on a single grid. Moreover, the solutions obtained with multiquadrics were better than those obtained with the multilevel method (even if up to 4225 points were used instead of 81). Of course, the support size of the compactly supported functions could be increased, but that would contradict the philosophy of wanting sparse matrices for efficient computation.

5. Concluding Remarks

In this paper we have described some numerical experiments with radial basis functions using the multilevel collocation method for partial differential equations. The application of a spline kernel smoothing delivers some limited improvements over the unsmoothed algorithm, but the difficulties caused by the oscillations of these kernels seem to outweigh their benefits. The Gauss-Weierstrass smoothing applied to a 2D problem, on the other hand, performed in a rather stable way, and yielded clear improvements over the unsmoothed method. The best results, however, were obtained using a direct approach with globally supported multiquadrics and a rather small set of collocation points. The use of globally supported functions should therefore not be ignored since they are also much easier to implement than a multilevel algorithm with added smoothing.

References

1. Beatson, R. K. and W. A. Light, Fast evaluation of radial basis functions: Methods for 2-dimensional polyharmonic splines, IMA J. Numer. Anal. 17 (1997), 343–372.

2. Fasshauer, G. E., Solving partial differential equations by collocation with radial basis functions, in Surface Fitting and Multiresolution Methods, A. Le Méhauté, C. Rabut, and L. L. Schumaker (eds), Vanderbilt University Press, Nashville TN, 1997, 131–138.

3. Fasshauer, G. E., On smoothing for multilevel approximation with radial basis functions, in Approximation Theory IX, Charles K. Chui and Larry L. Schumaker (eds), Vanderbilt University Press,, Nashville, TN, 1999, to appear.

4. Fasshauer, G. E., Solving differential equations with radial basis functions: multilevel methods and smoothing, Adv. in Comput. Math., to appear.

5. Fasshauer, G. E. and J. W. Jerome, Multistep approximation algorithms: Improved convergence rates through postconditioning with smoothing kernels, Adv. in Comput. Math., to appear.

6. Floater, M. S. and A. Iske, Multistep scattered data interpolation using compactly supported radial basis functions, J. Comput. Applied Math. 73 (1996), 65–78.

7. Franke, C. and R. Schaback, Solving partial differential equations by collocation using radial basis functions, Appl. Math. Comp. 93 (1998), 73–82.

8. Franke, C. and R. Schaback, Convergence orders of meshless collocation methods using radial basis functions, Adv. in Comput. Math. 8 (1998), 381–399.

9. Jerome, J. W., An adaptive Newton algorithm based on numerical inversion: regularization as postconditioner, Numer. Math. 47 (1985), 123–138.

10. 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.

11. Powell, M. J. D., Recent research at Cambridge on radial basis functions, preprint, 1998.

12. Schaback, R., On the efficiency of interpolation by radial basis functions, in Surface Fitting and Multiresolution Methods, A. Le Méhauté, C. Rabut, and L. L. Schumaker (eds), Vanderbilt University Press, Nashville TN, 1997, 309–318.

13. Wendland, H., Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. in Comput. Math. 4 (1995), 389–396.

14. Wendland, H., Meshless Galerkin approximation using radial basis functions, Math. Comp., to appear.

15. Wendland, H., Numerical solution of variational problems by radial basis functions, in Approximation Theory IX, Charles K. Chui and Larry L. Schumaker (eds), Vanderbilt University Press,, Nashville, TN, 1999, to appear.