The Latest RBF QR Algorithm Computer Science Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

RBF approximations have been useful to interpolate data on a sphere and on various other types of domains. Recently, RBF has been used to solve convection type PDEs over a sphere to spectral accuracy. There are two basic decisions to be made in such applications which are:

Which type of RBF to use?

What value should be selected for the shape parameter which is represented by ε and with flat basis function that is stretched out in radial direction with respect ε = 0.

The latest RBF-QR algorithm has made the computations substantially possible even for small values of ε. Results yielded from solving a convective PDE on a sphere are examined here against the various types of RBFs for the entire range of ε values. The methodology used to analyze the results was similar to that of the customary Fourier analysis in equispaced single dimension periodic settings. Particularly, it was found that high precision can be maintained even for very long time integrations. Moreover, it was determined that why RBFs occasionally provide higher accuracy as compared to spherical harmonics. RBF based methods in spherical geometries have the future potential in application areas of weather and climate modeling.

Table of Contents

Table of Contents 2

1.Introduction 3

2. RBF Methodology 7

2.1 RBF Interpolant Basic Form 7

2.2 Radial Basis Functions for PDEs 9

2.3 Time dependent PDE on a sphere 10

2.3.1 Test problem 10

4.Results-Comparisons among Different Types of RBFs 13

5. Discussion-Analysis of the numerical results via properties of the DMs 18

4.1 Single Dimension equi-spaced FD approximations 18

4.2 Dispersive errors for single dimension RBF approximations 20

4.3. Dispersive errors in convection over the sphere 21

6.Conclusions 24

References 25


Solving Partial Differential Equations (PDEs) in spherical geometries can be used in various application areas such as quantum mechanics, astrophysics, geophysics that includes climate and weather modeling, etc. It was lately demonstrated by Flyer and Wright (Flyer and Wright 2007) that a pseudo-spectral (PS) approach based on Radial Basis Function (RBF) are very adequately used to solve wave type partial differential equations, which include shallow water equations and pure convection equations, over spherical surfaces. On comparing the RBF approach to other spectral approaches (for example spherical harmonics methods, double Fourier methods and spectral element methods), many researchers found that RBF approach is very successful due to its following features:

Use of simple algebra, for instance, PDE test cases usually have less than 50 lines of Mat lab for their complete codes,

Quick classification from spherical surfaces to arbitrarily shaped surfaces and

Prospects to combine local refinement and spectral accuracy.

The RBFs mostly rely on shape parameter ε that approaches zero (ε  0) with respect to the limit of increasing flatness. Reducing the value of ε often enhances the subsequent accuracy to some level in case of the following one or two factors either inhibit or reverse this trend:

Occurrence of adverse numerical ill conditioning while employing a direct implementation of the procedure for radial basis function.

Occurrence of Runge phenomenon evocative to that of polynomial interpolation which is most likely to have adverse effects merely in the existence of boundaries or variable node densities besides mainly setting a limit for the reachable accuracy (Fornberg and Zuev 2007).

The above stated two factors are portrayed in the Figure 1, given below, in the three different cases extracted from the literature. The upper most row of the sub-plots demonstrates the variation in error with respect to ε while working out the Poisson equation on a unit disk as explained by Larsson and Fornberg (2003) in their study. The dotted line for comparison depicts the level of accuracy that has been acquired through the previously most precise procedures available which are Chebyshev pseudo-spectral radially and Fourier pseudo-spectral in the angular direction. All cases have a resolution of 16 nodes on the outer side and 48 nodes in the inner side. The Contour- Pade´ was used in the top right sub-plot in order to eliminate computational ill conditioning for small ε (Fornberg and Wright 2004). Such trends are also demonstrated in the middle row of sub-plots in relation to the interpolation of a Gaussian bell over a spherical surface in case of using n = 1849 nodes that are almost uniformly distributed. Fornberg and Piret (2007) in their study of a stable algorithm for flat RBFs on sphere introduced and used RBF-QR algorithm in lieu of the Contour Pade´ algorithm to get rid of the ill conditioning for the right sub-plot. The last sub-plot at the bottom that has been constructed with respect to the data extracted from the study on 'Transport schemes on a sphere using RBFs' by Flyer and Wright (2007) demonstrates the limit for the norm error once the cosine bell is convected over the sphere while solving the time based partial differential equation (Driscoll, 2004) which is explained later in this study.

The basic purpose of the current study is to further investigate this last test case particularly if the numerical ill conditioning for low values of ε has been eradicated. In this study, we will make use of both numerical computations and some novel analysis to fulfill the following tasks.

Explore the variations in the accuracy through

Type of Radial Basis Function

Parameter of Shape for entire range of ε > 0

Length of time integration

Explore the case of ε = 0 as this is often identical to using spherical harmonics in lieu of radial basis functions.

Poisson equation on unit disk




Ill-conditioning elliminatedA very concise introduction to radial basis function interpolation is provided in the next section along with quoting some related results from the literature for example the importance of the flat basis function limit. Then the section 3 includes the convective partial differential equation and results making use of radial basis function discretization. In the context of our study, the performance disparity between the types of smooth radial basis function are found to be insignificant however the performances of the non-smooth radial basis function are further lower. Section 4 analyzes and explains some of these results, generalizing them to scattered node cases over the sphere, the Fourier-based arguments are usually used on periodic problems with equal spacing in single dimension (1-D). Section 5 focuses on summaries and conclusions as it discusses the reason behind the superiority of smooth basis functions over their non-smooth counterparts, in the present context of long time integration, irrespective of the smoothness of the solution being convected.

Interpolation on Sphere



Time dependent convection on Sphere


Figure 1: Radial-Basis-Function errors in 3 different applications being illustrated as functions of ε, the shape parameter. The data used to construct the subplots categorized in the above three rows have been reproduced from the studies of Larsson and Fornberg (2005), Fornberg and Piret (2007) and Flyer and Wright (2007), respectively. MQ refers to Multi-Quadratic RBF, IQ refers to Inverse Quadratic RBF and GA refers Gaussian RBF which are further discussed in the following section. The current work will produce the missing right hand side sub-plot of the third row.

2. RBF Methodology

2.1 RBF Interpolant Basic Form

The basic form of an RBF interpolant is represented by


s(x) = Σ λiΦ(â•‘x - xi‌‌‌‌‌â•‘) ___________(1)


where â•‘.â•‘ represents the Euclidean norm. The expansion coefficients λi should satisfy

Aλ = f _______________(2)

so as to take the values of fi at locations xi where i = 1, 2, 3 ,…., n for the Euclidean norm.

Ai,j = Φ(â•‘x - xi‌‌‌‌‌â•‘) are the entries to the matrix A. The numerical use of equation (2) followed by equation (1) is denoted by RBF-Direct. This study will focus on radial functions Φ(r) provided in the Table 1 given below. The parameter ε is called the shape parameter which is included in all but piece-wise smooth global cases TPS and CU. The Wendland functions (Wendland 2005) listed in the table are the functions of lowest degree which assure that the matrix A is positive definite in 2 and 3 dimensions for all distinct node locations. IQ, AMQ and GA hold this property in all dimensions. The matrix A is still symmetric in the cases: CU, MQ and TPS. Non-singularity remains ascertained for MQ however positive definiteness is eliminated. Additional issues are encountered in the TPS and CU cases. Generalized variations of equation (1) are for MQ.

Table 1:

Name of Radial Basis Function



Smooth, Global:

Inverse Quadratic

Inverse Multi-Quadratic











Piecewise Smooth, Global:

Thin Plate Spline






Piecewise Smooth, Compact:

(for 0< r < 1/ε, for r> 1/ε equal to 0)

Wendland type order 2

Wendland type order 4

Wendland type order 6




(1-εr)4(4εr +1)



The shape parameter ε is used to control the flatness of radial basis functions. For Wendland functions, their orders denote their degree of smoothness C2, C4 and C6, respectively.

n n

s(x) = α + ΣλiΦ(â•‘x - xiâ•‘) if Σλi = 0 ___________________(3)

i=1 i=1

and in the case of 3-D, for TPS and CU


s(x) = α + βx + γy + δz + ΣλiΦ(â•‘x - xiâ•‘) _______________(4)


n n n n

Only if, Σλi = Σλixi = Σλiyi = Σλizi = 0

i=1 i=1 i=1 i=1

The most significant aspect of Wendland functions is their compact support if the value of ε is large, that results in sparse A-matrices, and thus, with potentiality for high computational speeds in both solving and evaluating, such as in case of using conjugate gradient type methods. In the context of present case with nodes on the unit sphere, if the value of ε is less than 0.5 then all sparsity is lost.

2.2 Radial Basis Functions for PDEs

Kansa presented collocation with RBFs in 1990 as mechanism to approximate spatial derivatives and therefore to solve PDEs, numerically (Kansa 1990). This approach is, by convention, spectrally accurate for smooth RBFs (Buhmann 1993). Other significant benefit of this RBF approach, in comparison to finite element, finite difference and finite volume methods, is that it supersedes the typically intricate task of developing computational meshes over irregular domains with the scattering computational nodes that are more convenient.

According to Driscoll and Fornberg (2002), the interpolant typically converges to Lagrange's interpolation polynomial for globally smooth RBFs in single dimension in the case of flat basis function limit approaching to zero (ε  0). This leads to the perception that RBF approach for PDEs is a generalization of the pseudo-spectral (PS) method to irregular domains and to scattered nodes (Boyd 2001). Thus, this approach now often categorized as the RBF-PS method. The limit ε  0 will regenerate the Fourier-PS method in single dimension periodic setting (Driscoll and Fornberg 2002). Fornberg and Piret (2007) found the corresponding limit to agree with the spherical harmonics pseudo-spectral (SPH-PS) method for the nodes on a sphere. As mentioned previously in the introduction that Flyer and Wright (2007) were the first to solve purely convective, that is non-dissipative, PDEs on spherical surface through using RBFs. Ill-conditioning inhibited their numerical explorations from being expanded to even arbitrarily small values of ε since they implemented their scheme based on the direct use of equation (2) that is followed by equation (1).

2.3 Time dependent PDE on a sphere

2.3.1 Test problem

The test problem used in this study is described as a solid body that rotates/convects around an axis inclined by the angle 'α' corresponding to the polar axis (Figure 2a). The spherical coordinates are defined below according to the convention in many applications.

The standard spherical coordinates differ from the above definition in the use of latitude (in comparison to co-latitude) as well as in the direction referred by θ and φ. The governing PDE in the present coordinates becomes


One complete revolution will relate to the time t = 2π.

The pole singularities penetrate by means of the factor 'tan θ' where θ = + π/2. These singularities are not tangible, but emerge only as result of the (θ,φ)-system that is singular at these points. These singularities will be lost with formation of RBF differentiation matrix (DM that represents the spatial operator in terms of the node values of u. DM denotes the physical operator but not the coordinate system as it occurred to be expressed in whilst an in-between derivation step).

The figure 3 shows the initial condition used in this study which can be defined in Cartesian coordinates as:

Figure 2b: Spherical latitude-longitude-type coordinate system.

Figure 2a: Flow directions in the 'solid body' convection test problem.

Figure 3b: The illustration of cosine bell in grey scale over the spherical surface from positive x direction.

Figure 3a: The illustration of the cosine bell as a function of x based on the work of Fasshauer (2007).

Figure 3c: The illustration of the cosine bell on an 'unrolled' φ, θ-plane marked with the n = 1849 ME node locations (Fasshauer 2007).

Where R is constrained to the surface of the unit sphere and is equal to 1/3.

Cosine bells are considered as the standard initial conditions for convective test calculations over a sphere (Jakob-Chien, Hack and Williamson 1995) due to various reasons that include:

The dispersive errors after convection are easy to display and interpret due to their compact support.

Testing can be conducted over different spatial scales due to the easy-to-change peak width R.

Very high order methods are prevented, by the discontinuous 2nd derivative at the base of the bell, from generating misleadingly good results as opposed to those expected from cases having more physically relevant data.

Results-Comparisons among Different Types of RBFs

The data points in Figure 3c have been used as the input data and represented in the table below.

Table 2: The Input Data Table.

N = 43 x 43 = 1849




































































































Figure 4: The SPH interpolant to the cosine bell for 1849 ME node points.

Figure 5. The max norm errors if the initial cosine bell is brought to 1849 node RBF representation. The ε values related to 0 percent and 90 percent sparsity of the matrix 'A 'are marked for the W6 cases. For all cases, the direct interpolation is compared with the least square approach. The divergence of RBF-Direct for the declining ε is portrayed in sub-plots a and b (shown as almost vertical solid lines at ε = 1). Similar divergence occurs in above sub-plots c and d and in the figures 6, 8, and 9 given below but it is not marked explicitly.

Figure 6: Errors when times t =10 and t = 10, 000 in the form of functions of ε for (a) IMQ and (b) W6. The error (ε-independent) of the commencing SPH representation of the initial data at

t=0 is given by the thin dashed line.

Figure 7: Progression of the error with time when (a) 0 < t < 10 and (b) 0 < t < 10, 000. The values of ε GA, IMQ, and W6 were close to zero while ε is not included in the case of TPS.

It has been found that some types of RBFs tend to yield more precise results than the rest. MQ is considered far better than GA merely because GA does not support any 'exact polynomial reproduction' properties on infinite lattices where ε is greater than zero (Buhmann 2003). The precision acquired through different RBF methods is highly based over problems such as the interpolation of non-smooth data poses various demands on the methods as opposed to long term solution of convective PDEs. The combination of RBF-Direct and RBF-QR allows us to run the convective test problem for the entire set of RBFs, for the complete range of ε (i.e., > 0), yielding the results that are shown in Figures 8 and 9. As observed previously, smooth RBFs provide exceptional precision even for prolonged integrations.

Figure 8: The errors, for the entire range of RBF types considered in this study, as a function of ε at t = 10. In GA, the ill-conditioning for RBF-Direct appears a bit earlier as compared to other smooth RBF types which leaves a small gap between the ranges of RBF-QR and RBF-Direct as shown in Figure 9. The MQ results are significantly more accurate as compared to those for other RBF types in case of larger values of ε but not as good as the results achieved by all smooth RBFs for small values of ε.

Figure 9: The errors, for the entire range of RBF types considered in this study, as a function of ε at t = 10,000.

Discussion-Analysis of the numerical results via properties of the DMs

Spectral Analysis explains the above obtained numerical results, which is conceptually akin to Fourier analysis of equally spaced finite difference and RBF approximations to the model equation

∂u + ∂u = 0 ________________________(7)

∂t + ∂x

4.1 Single Dimension equi-spaced FD approximations

It is already known that how accuracy of a FD scheme can be explained through its behavior on individual Fourier modes. Since

d eiÏŽx = iÏŽeiÏŽx


The eigen-function to the operator d/dx is given by u(x) = eiÏŽx. For the discrete situation with 'h' as the grid spacing, ÏŽ Ñ” [-ÏŽmax, ÏŽmax] is the possible frequency range, when ÏŽmax = π/h, due to aliasing. Similarly, applying a finite difference 2, that is the 2nd order FD at the center, scheme to u(x) = eiÏŽx provides:

u(x+h) - u(x-h) = eiÏŽ(x+h) - eiÏŽ(x-h) = i sinÏŽh eiÏŽx ___________________(8)

2h 2h h

This means that Fourier mode is also an eigen function however the eigen-value has shifted from iÏŽ to (sinÏŽh)/h. If we just ignore 'i' for now, then the eigen-values would be given by fPS(ÏŽ) = ÏŽ and fFD2(ÏŽ) = (sinÏŽh)/h which are illustrated in the figure 10 combined with similar factors even for some higher order FD methods. Whereas the PD method rightly treats all modes presented on a h-spacing grid, notable errors are found in lower order FD methods in all modes. These errors relate to the errors in the phase speed of the traveling waves in case of time integration of equation (7).

Figure 10: The effect of d/dx and FD approximations of its different orders when enforced to a basic Fourier mode eiÏŽx. ÏŽmax = π/h denotes the highest mode that can be present on h-spacing grid. (Reproduced from Fornberg, 1996).

The number of high modes of phase increases with the length of time integration and thus vanishes while improving accuracy. The severity of this degradation over time rises with the decline in the order of the finite difference approximation.

Figure 11 shows the solutions for narrow Gaussian initial bell through direct numerical evaluation. It is obvious that the accuracy of low order FD2 method is lost very soon while the high order FD methods are able to hold the pulse integrity for some longer period but severe dispersive trailing wave train occurs in all cases.

Figure 11: The narrow Gaussian pulse's integrity (portrayed as dots on a grid with h = 1) being integrated exactly in time by employing FD2, FD4 and FD6 methods, respectively, in space.

4.2 Dispersive errors for single dimension RBF approximations

In case of ε  0, the curves in Figure 10 are observed to swiftly approach the perfect PS straight line case. The equation below was numerically evaluated for the same test cases as defined for FD methods in Figure 11 so as to acquire the results for IQ-RBF shown in Figure 12. Although large value of ε is quite acceptable at short times however smaller values of ε are evidently seen to support high long term accuracy. Despite the RBF scheme is spectrally precise for all values of ε, however, these results for constant h=1 and different ε are implicative of those for increasing order FD schemes which is due to the similarities just depicted in how the eigen-values of the FD approximations differ with the Fourier frequency.

4.3. Dispersive errors in convection over the sphere

SPH modes Yvμ(x), μ < μmax, ν = -ν,….,-1,0,1,…,ν, in this geometry, develop a counterpart to a truncated set of Fourier modes in a periodic single dimension case. Figure 13 represents some low modes up to μmax = 4. All the eigen-values of DM will be purely imaginary.

Figure 12: Results for IQ-RBF with different ε-values for all the test cases that were used in

Figure 11.

Figure 13: The illustration for the SPH basis functions Yvμ for orders μ < μmax = 4. The dashed

lines portray where they change the sign.

Fornberg and Piret (2007) found that as the shape parameter ε approaches to zero, the space spanned through n = μ2max globally smooth Radial Basis Function will agree completely with the spherical harmonics space for μ < μmax. Now in this limit of shape parameter we can add to these observations that the n = n = μ2max eigen values of Radial Basis Function DM will reach:

μmax + 1 eigen-values 0

μmax eigen-values +i and equally many -i

μmax -1 eigen-values +2i and equally many -2i




so on till

eigen-values +μi and also one eigen-value -μi _____________________(16)

This can be concluded from the further 2 observations:

The truncated expansions of spherical harmonics constitute a closed set in relation to any rotation of the coordinate system that yields the result such like this independent on the value of α in equation (5) provided previously in section 2.3.1 of this paper.

An examination of the patterns illustrated in Figure 13 above, for 'α' equal to zero, depicts that the functions μmax + 1 in the center column, ν = 0, are not influenced by any movement around the polar axis that leads towards μmax + 1 eigen-values 0. The functions, μmax for ν = -1 and similarly μmax for ν = +1, iterate themselves after 1 complete revolution whereas the next group of functions repeat themselves after half revolution and the next group of functions repeat themselves after one-third revolution and so on and so forth.

Maintaining α equal to zero and making the graphical representation less scattered while selecting μmax = 23 with n = μ2 max = (23)2 = 529 as opposed to μmax = 43 with n = μ2 max = (43)2 = 1849 that is used else where in this study, the eigen-values - as functions of the parameters, μ and ν, of spherical harmonics in the case of ε equal to zero - become as represented in figure 14. This flat triangular section of a plane relates to the pseudo-spectral straight line shown in Figure 10 and the likewise straight line in case of the shape parameter ε equals to zero. Increasing the value of the shape parameter from zero will result in deviations from the ideal pattern of eigen-value represented in Figure 14. Figure 15a illustrates differently in comparison to Figure 10 that how the eigen-values in the case of single dimension FD change with respect to the order p of the FD schemes. The range of eigen-values along the imaginary axis reduces by a factor of π as we move from p = ∞ to p = 2. In the same manner, Figure 15b illustrates the variation in eigen values of spectral harmonics (ε = 0) as the value of the shape parameter increases from zero.

The largest eigen-value reduces swiftly with the increase in the value of the shape parameter, which leaves the lowest eigen-values unchanged for the longest time.

Figure 15: The eigen-values as computed by using IQ Radial Basis Function in the limit of ε = 0 to the convective operator in the sphere case. Fornberg, Wright and Larsson (2004) associated these eigen-values with values of μ and ν as shown here.

Figure 15: (a) The n = 51 single dimension eigen-values of convection Partial Differential Equation for distinct orders of FD schemes if h = 1 and for a periodic domain from -25 to 25 instead of -∞ to ∞. (b) The eigen-values for n = 529 sphere convection problem for various values of the shape parameter for MQ and in case the results for other smooth Radial Basis Function types are quite alike). The vertical scales not shown as the current problem of concern is to qualitatively demonstrate how the limits in the 2 cases change with the order of FD and with the value of the shape parameter, respectively.


The numerical results conclude that very long time integrations are entirely viable with the RBF method however the shape parameter ε then has to be very quite low, particularly if the integration time is increased as it makes the computations beyond reach for the RBF-Direct method when employing the standard 16-digit arithmetic of double precision. All types of smooth RBFs virtually provide the same precision however the piecewise smooth RBFs are not at all competent in the current scenario. These general results are not affected by the smoothness of the convected solution, following from the properties of the DMs.