- Email: [email protected]
Hermite function interpolation on a finite uniform grid: Defeating the Runge phenomenon and replacing radial basis functions
Applied Mathematics Letters 26 (2013) 995–997
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Research announcement
Hermite function interpolation on a finite uniform grid: Defeating the Runge phenomenon and replacing radial basis functions John P. Boyd ∗ , Luis F. Alfaro Department of Atmospheric, Oceanic and Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109-2143, United States
article
info
Article history: Received 15 March 2013 Received in revised form 3 May 2013 Accepted 3 May 2013 Keywords: Hermite function Radial basis functions Interpolation Runge phenomenon
abstract We show that Hermite functions are successful for interpolation on a finite interval, even on a uniform grid where polynomial interpolation fails. The Runge Phenomenon is not completely abolished, but is greatly diminished. Finite interval Hermite interpolation is ill-conditioned and therefore limited to a maximum of about 250 interpolation points on a uniform univariate grid, but it is still far superior to the approach using Gaussian radial basis functions (RBFs). The Hermite functions are a complete spectral basis for the infinite interval; the motivation for employing them here derives from a careful study that showed, for small RBF shape parameter α and small number of points N, that Gaussian RBF cardinal functions can be accurately approximated by the product of a Gaussian with the usual polynomial cardinal function. Direct comparisons show that when N and α are not both small, Hermite interpolation is greatly superior in accuracy, condition number and efficiency to the RBF methods that inspired it. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Boyd and Alfaro [1] show that on a uniform grid in one dimension, the radial basis function (RBF) cardinal function, poly poly CjRBF , is approximately equal to the product (K(x)/K(xj ))Cj (x) where Cj (x) is the usual Lagrangian polynomial cardinal function and K(x) is a function which is well approximated by a Gaussian. In this work, we jettison the RBF paradigm and instead treat polynomials multiplied by a Gaussian as the fundamental approximation. This basis is equivalent to applying a series of Hermite functions, ψn (y), in approximation on a finite interval, which can be normalized to x ∈ [−1, 1] without loss of generality. One may then apply the usual interpolation conditions fN (xj ) = f (xj ) for all j = 1, 2, . . . , N where f (x) is the ‘‘target’’ function to be approximated and fN is the approximation fN ≡
N −1
an ψn (y[x])
(1)
n =0
where ψn (y) is a normalized Hermite function defined by the recursion
ψ0 (y) ≡ π −1/4 exp(−(1/2)y2 ) √ ψ1 (y) ≡ π −1/4 2 y exp(−(1/2)y2 ) ∗
Corresponding author. E-mail addresses: [email protected], [email protected] (J.P. Boyd).
0893-9659/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.aml.2013.05.011
(2) (3)
996
J.P. Boyd, L.F. Alfaro / Applied Mathematics Letters 26 (2013) 995–997
Fig. 1. Isolines of log 10(L∞ ) error. Comparison of Gaussian RBFs (right) with Hermite pseudo-RBFs (left) for f (x) = (1 + x)/(1 + 25x2 ). Maple high precision arithmetic was used to suppress the distortions of roundoff error.
ψn+1 (y) =
2 n+1
y ψn (y) −
n n+1
ψn−1 (y).
(4)
The argument of the Hermite functions, y, is scaled such that exp(−(1/2)y2 ) = exp(−[(N − 1)/2] log(4)i2 x2 ), implying y=
√
log(4) N − 1 i x
(5)
where i is a user-choosable parameter; the reason for adopting this somewhat odd convention is to ensure that the exponential growth of the polynomial cardinal functions as x → ±1 is exactly canceled when i = 1. (The optimum choice of i is indeed 1 or slightly smaller.) Hermite interpolation is much cheaper than using Gaussian RBFs because factoring the RBF interpolation matrix costs O(N 3 ) operations. For Hermite functions, it is vastly cheaper to use the ‘‘first barycentric’’ form of the interpolant, N wj fN = Ω (x) f (xj ), x − xj j =1
Ω ≡ exp(−[(N − 1)/2] log(4)i x ) 2 2
N j =1
(x − xj ),
wj =
dΩ dx
(xj )
−1
.
(6)
This costs only O(N ) operations per point or O(N 2 ) for a set of N points, an order of magnitude saving over matrix inversion. However, the barycentric option is not possible for RBFs because the RBF cardinal functions are not known in explicit form. The contour plots in Fig. 1 show the isolines of the maximum pointwise errors in approximating Runge’s exemplary f (x) on a uniform grid (where polynomial interpolation diverges). RBFs (right) give an error as small as 10−6 for the largest N and for the optimum choice of the RBF parameter α . However, Hermite functions (left) give errors a million times smaller! The Hermite function interpolant shares two defects with polynomial and RBF interpolation: the Runge phenomenon and ill-conditioning. The latter is not completely eliminated by either Hermite or RBF interpolation, but the size of the Runge Zone (where singularities of f (x) imply divergence) is contracted by an order of magnitude, and may be further manipulated by varying the ‘‘shape’’ parameter (RBFs) or i (Hermite). Fig. 2 shows that ill-conditioning limits the best attainable accuracy to roughly the square root of machine epsilon. An analysis presented in [2] shows that when the Gaussian width of the RBFs is chosen to match that of the Hermite functions, the logarithm of the condition number is about seven times larger for RBFs—a huge advantage for Hermite function interpolation. The dashed curve shows that although the rate of Hermite function convergence is geometric, exponential convergence, Hermite interpolation on a uniform grid is not as accurate as polynomial interpolation on the highly nonuniform Chebyshev grid. But the nonuniformity of the Chebyshev grid creates its own drawbacks [3,1]. Our forthcoming full-length article will provide theory, maps of the Runge Zone, connections with Fourier Extension theory [4–7,2] and other details omitted in this brief announcement.
J.P. Boyd, L.F. Alfaro / Applied Mathematics Letters 26 (2013) 995–997
997
10-6
10-9
10-12
10-15
10-18
10-21 100
200
300
400 N
500
600
700
Fig. 2. Hermite function interpolation (on a uniform grid) of f (x) = (1 + x)/(1 + 25x2 ) as computed in Maple using 16 decimal digits of precision (red, long dashes), 32 digits (blue, short dashes), and 64 digits (solid black). For comparison, polynomial interpolation on the very nonuniform (and very wellconditioned) Chebyshev grid in 16-digit precision is also shown as the thin, dashed green line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Acknowledgments This work was supported by the National Science Foundation through grant OCE 1059703 and by the Undergraduate Research Opportunities Program (UROP) at the University of Michigan. We thank the two reviewers for their comments. References [1] J.P. Boyd, L. Alfaro, Constructing RBF cardinal functions from polynomial and sinc cardinal functions, J. Comput. Appl. Math. (2013) in preparation. [2] J.P. Boyd, L. Alfaro, Hermite function interpolation on a finite interval: defeating the Runge phenomenon and the partial obsolescence of radial basis functions, Comput. Math. Appl. (2013) in preparation. [3] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001, p. 665. [4] B. Adcock, D. Huybrechs, On the resolution power of Fourier extensions for oscillatory functions, J. Comput. Appl. Math. (2012) submitted for publication. [5] B. Adcock, D. Huybrechs, J. Martin-Vaquero, On the numerical stability of Fourier extensions for oscillatory functions, Found. Comput. Math. (2012) submitted for publication. [6] J.P. Boyd, Fourier embedded domain methods: extending a function defined on an irregular region to a rectangle so that the extension is spatially periodic and C ∞ , Appl. Math. Comput. 161 (2005) 591–597. [7] J.P. Boyd, A comparison of numerical algorithms for Fourier extension of the first, second and third kinds, J. Comput. Phys. 178 (2002) 118–160.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Research announcement
Hermite function interpolation on a finite uniform grid: Defeating the Runge phenomenon and replacing radial basis functions John P. Boyd ∗ , Luis F. Alfaro Department of Atmospheric, Oceanic and Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109-2143, United States
article
info
Article history: Received 15 March 2013 Received in revised form 3 May 2013 Accepted 3 May 2013 Keywords: Hermite function Radial basis functions Interpolation Runge phenomenon
abstract We show that Hermite functions are successful for interpolation on a finite interval, even on a uniform grid where polynomial interpolation fails. The Runge Phenomenon is not completely abolished, but is greatly diminished. Finite interval Hermite interpolation is ill-conditioned and therefore limited to a maximum of about 250 interpolation points on a uniform univariate grid, but it is still far superior to the approach using Gaussian radial basis functions (RBFs). The Hermite functions are a complete spectral basis for the infinite interval; the motivation for employing them here derives from a careful study that showed, for small RBF shape parameter α and small number of points N, that Gaussian RBF cardinal functions can be accurately approximated by the product of a Gaussian with the usual polynomial cardinal function. Direct comparisons show that when N and α are not both small, Hermite interpolation is greatly superior in accuracy, condition number and efficiency to the RBF methods that inspired it. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Boyd and Alfaro [1] show that on a uniform grid in one dimension, the radial basis function (RBF) cardinal function, poly poly CjRBF , is approximately equal to the product (K(x)/K(xj ))Cj (x) where Cj (x) is the usual Lagrangian polynomial cardinal function and K(x) is a function which is well approximated by a Gaussian. In this work, we jettison the RBF paradigm and instead treat polynomials multiplied by a Gaussian as the fundamental approximation. This basis is equivalent to applying a series of Hermite functions, ψn (y), in approximation on a finite interval, which can be normalized to x ∈ [−1, 1] without loss of generality. One may then apply the usual interpolation conditions fN (xj ) = f (xj ) for all j = 1, 2, . . . , N where f (x) is the ‘‘target’’ function to be approximated and fN is the approximation fN ≡
N −1
an ψn (y[x])
(1)
n =0
where ψn (y) is a normalized Hermite function defined by the recursion
ψ0 (y) ≡ π −1/4 exp(−(1/2)y2 ) √ ψ1 (y) ≡ π −1/4 2 y exp(−(1/2)y2 ) ∗
Corresponding author. E-mail addresses: [email protected], [email protected] (J.P. Boyd).
0893-9659/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.aml.2013.05.011
(2) (3)
996
J.P. Boyd, L.F. Alfaro / Applied Mathematics Letters 26 (2013) 995–997
Fig. 1. Isolines of log 10(L∞ ) error. Comparison of Gaussian RBFs (right) with Hermite pseudo-RBFs (left) for f (x) = (1 + x)/(1 + 25x2 ). Maple high precision arithmetic was used to suppress the distortions of roundoff error.
ψn+1 (y) =
2 n+1
y ψn (y) −
n n+1
ψn−1 (y).
(4)
The argument of the Hermite functions, y, is scaled such that exp(−(1/2)y2 ) = exp(−[(N − 1)/2] log(4)i2 x2 ), implying y=
√
log(4) N − 1 i x
(5)
where i is a user-choosable parameter; the reason for adopting this somewhat odd convention is to ensure that the exponential growth of the polynomial cardinal functions as x → ±1 is exactly canceled when i = 1. (The optimum choice of i is indeed 1 or slightly smaller.) Hermite interpolation is much cheaper than using Gaussian RBFs because factoring the RBF interpolation matrix costs O(N 3 ) operations. For Hermite functions, it is vastly cheaper to use the ‘‘first barycentric’’ form of the interpolant, N wj fN = Ω (x) f (xj ), x − xj j =1
Ω ≡ exp(−[(N − 1)/2] log(4)i x ) 2 2
N j =1
(x − xj ),
wj =
dΩ dx
(xj )
−1
.
(6)
This costs only O(N ) operations per point or O(N 2 ) for a set of N points, an order of magnitude saving over matrix inversion. However, the barycentric option is not possible for RBFs because the RBF cardinal functions are not known in explicit form. The contour plots in Fig. 1 show the isolines of the maximum pointwise errors in approximating Runge’s exemplary f (x) on a uniform grid (where polynomial interpolation diverges). RBFs (right) give an error as small as 10−6 for the largest N and for the optimum choice of the RBF parameter α . However, Hermite functions (left) give errors a million times smaller! The Hermite function interpolant shares two defects with polynomial and RBF interpolation: the Runge phenomenon and ill-conditioning. The latter is not completely eliminated by either Hermite or RBF interpolation, but the size of the Runge Zone (where singularities of f (x) imply divergence) is contracted by an order of magnitude, and may be further manipulated by varying the ‘‘shape’’ parameter (RBFs) or i (Hermite). Fig. 2 shows that ill-conditioning limits the best attainable accuracy to roughly the square root of machine epsilon. An analysis presented in [2] shows that when the Gaussian width of the RBFs is chosen to match that of the Hermite functions, the logarithm of the condition number is about seven times larger for RBFs—a huge advantage for Hermite function interpolation. The dashed curve shows that although the rate of Hermite function convergence is geometric, exponential convergence, Hermite interpolation on a uniform grid is not as accurate as polynomial interpolation on the highly nonuniform Chebyshev grid. But the nonuniformity of the Chebyshev grid creates its own drawbacks [3,1]. Our forthcoming full-length article will provide theory, maps of the Runge Zone, connections with Fourier Extension theory [4–7,2] and other details omitted in this brief announcement.
J.P. Boyd, L.F. Alfaro / Applied Mathematics Letters 26 (2013) 995–997
997
10-6
10-9
10-12
10-15
10-18
10-21 100
200
300
400 N
500
600
700
Fig. 2. Hermite function interpolation (on a uniform grid) of f (x) = (1 + x)/(1 + 25x2 ) as computed in Maple using 16 decimal digits of precision (red, long dashes), 32 digits (blue, short dashes), and 64 digits (solid black). For comparison, polynomial interpolation on the very nonuniform (and very wellconditioned) Chebyshev grid in 16-digit precision is also shown as the thin, dashed green line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Acknowledgments This work was supported by the National Science Foundation through grant OCE 1059703 and by the Undergraduate Research Opportunities Program (UROP) at the University of Michigan. We thank the two reviewers for their comments. References [1] J.P. Boyd, L. Alfaro, Constructing RBF cardinal functions from polynomial and sinc cardinal functions, J. Comput. Appl. Math. (2013) in preparation. [2] J.P. Boyd, L. Alfaro, Hermite function interpolation on a finite interval: defeating the Runge phenomenon and the partial obsolescence of radial basis functions, Comput. Math. Appl. (2013) in preparation. [3] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001, p. 665. [4] B. Adcock, D. Huybrechs, On the resolution power of Fourier extensions for oscillatory functions, J. Comput. Appl. Math. (2012) submitted for publication. [5] B. Adcock, D. Huybrechs, J. Martin-Vaquero, On the numerical stability of Fourier extensions for oscillatory functions, Found. Comput. Math. (2012) submitted for publication. [6] J.P. Boyd, Fourier embedded domain methods: extending a function defined on an irregular region to a rectangle so that the extension is spatially periodic and C ∞ , Appl. Math. Comput. 161 (2005) 591–597. [7] J.P. Boyd, A comparison of numerical algorithms for Fourier extension of the first, second and third kinds, J. Comput. Phys. 178 (2002) 118–160.