Improved conditioning of the Floater--Hormann interpolants

TL;DR

This paper introduces a modification to Floater--Hormann rational interpolants by adding local polynomial interpolants at the interval ends, which improves their conditioning and enables higher approximation orders in practical applications.

Contribution

The paper proposes a new variant of Floater--Hormann interpolants with enhanced conditioning through end-based local polynomial interpolants, facilitating higher approximation accuracy.

Findings

Improved conditioning of interpolants with the proposed modification.

Enabling higher approximation orders in practical computations.

Reduction in rounding error amplification.

Abstract

The Floater--Hormann family of rational interpolants do not have spurious poles or unattainable points, are efficient to calculate, and have arbitrarily high approximation orders. One concern when using them is that the amplification of rounding errors increases with approximation order, and can make balancing the interpolation error and rounding error difficult. This article proposes to modify the Floater--Hormann interpolants by including additional local polynomial interpolants at the ends of the interval. This appears to improve the conditioning of the interpolants and allow higher approximation orders to be used in practice.

Figures

Tables

Table 1. Table 1: L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} errors of r ( d ) ​ ( x ) superscript 𝑟 𝑑 𝑥 r^{(d)}(x) and r ( d , e ) ​ ( x ) superscript 𝑟 𝑑 𝑒 𝑥 r^{(d,e)}(x) on the interval [ − 5 , 5 ] 5 5 [-5,5] with n 𝑛 n equispaced nodes for f ​ ( x ) = 1 / ( 1 + x 2 ) 𝑓 𝑥 1 1 superscript 𝑥 2 f(x)=1/(1+x^{2}) . For r ( d ) ​ ( x ) superscript 𝑟 𝑑 𝑥 r^{(d)}(x) , d 𝑑 d is the optimal value reported by Floater and Hormann 2007floater . For r ( d , e ) ​ ( x ) superscript 𝑟 𝑑 𝑒 𝑥 r^{(d,e)}(x) , d , e = 14 , 4 formulae-sequence 𝑑 𝑒 14 4 d,e=14,4 except when d 𝑑 d is constrained by n 𝑛 n .

$n$	$d$	$L^{\infty }$, $r^{\left(d\right)}\left(x\right)$	$L^1$, $r^{\left(d\right)}\left(x\right)$	$d,e$	$L^{\infty }$, $r^{(d,e)}\left(x\right)$	$L^1$, $r^{(d,e)}\left(x\right)$
10	0	$3.606\times {10}^{-2}$	$1.601\times {10}^{-1}$	10, 4	$3.005\times {10}^{-2}$	$1.243\times {10}^{-1}$
20	1	$1.536\times {10}^{-3}$	$6.656\times {10}^{-3}$	14, 4	$1.674\times {10}^{-3}$	$4.519\times {10}^{-3}$
40	3	$4.307\times {10}^{-6}$	$1.306\times {10}^{-5}$	14, 4	$3.463\times {10}^{-6}$	$1.220\times {10}^{-5}$
80	7	$2.038\times {10}^{-10}$	$8.003\times {10}^{-11}$	14, 4	$1.214\times {10}^{-11}$	$4.684\times {10}^{-11}$
160	10	$1.887\times {10}^{-15}$	$9.230\times {10}^{-16}$	14, 4	$1.887\times {10}^{-15}$	$9.226\times {10}^{-16}$

Equations

$${\chi }_{i,j}(x)=(-1{)}^i\prod _{k=i}^j\frac{1}{x-x_k}.$$

$$r^{(d)}(x)=\frac{{\sum }_{i=0}^{n-d}{\chi }_{i,i+d}(x)p_{i,i+d}(x)}{{\sum }_{i=0}^{n-d}{\chi }_{i,i+d}(x)}.$$

$$p_{i,j}(x)=\sum _{k=i}^j\prod _{\begin{array}{c}l=i\\ l\ne k\end{array}}^j\frac{x-x_l}{x_k-x_l}{y}_k.$$

$$t^{(d)}(x)=\sum _{i=0}^{n-d}\sum _{j=i}^{i+d}\frac{(-1{)}^i}{x-x_j}\prod _{\begin{array}{c}k=i\\ k\ne j\end{array}}^{i+d}\frac{1}{x_j-x_k}{y}_j.$$

$${\omega }_{i,j,k}=\prod _{\begin{array}{c}l=i\\ l\ne j\end{array}}^k\frac{1}{x_j-x_l}.$$

$${\xi }_j^{(d)}=\sum _{i=\max (0,j-d)}^{\min (j,n-d)}(-1{)}^i{\omega }_{i,j,i+d}$$

$$t^{(d)}(x)=\sum _{j=0}^n\frac{{\xi }_j^{(d)}}{x-x_j}{y}_j.$$

$$r^{(d)}(x)$$

$$=\sum _{j=0}^n{\beta }_j^{(d)}(x)y_j$$

$${\varphi }_i^{(d)}(x)$$

$${\psi }_i^{(d)}(x)$$

$$r^{(d,e)}(x)=\frac{{\sum }_{i\in I_1}{\varphi }_i^{(d)}(x)p_{0,i}(x)+{\sum }_{i\in I_2}{\chi }_{i,i+d}(x)p_{i,i+d}(x)+{\sum }_{i\in I_3}{\psi }_i^{(d)}(x)p_{i,n}(x)}{{\sum }_{i\in I_1}{\varphi }_i^{(d)}(x)+{\sum }_{i\in I_2}{\chi }_{i,i+d}(x)+{\sum }_{i\in I_3}{\psi }_i^{(d)}(x)}.$$

$$t^{(d,e)}(x)=$$

$$+\sum _{i=n-d+1}^{n-d+e}\frac{(-1{)}^i}{(x-x_n{)}^{i-n+d}}\sum _{j=i}^n\frac{{\omega }_{i,j,n}}{x-x_j}{y}_j$$

$${\zeta }_j^{(d,e)}(x)$$

$$\displaystyle=-\frac{\omega_{0,j,u}}{x-x_{0}}\bigg{(}1-\frac{x_{j}-x_{u}}{x-x_{0}}\bigg{(}\dots\bigg{(}1-\frac{x_{j}-x_{l+1}}{x-x_{0}}\bigg{)}\bigg{)}\bigg{)}$$

$${\eta }_j^{(d,e)}(x)$$

$$\displaystyle=(-1)^{l}\frac{\omega_{l,j,n}}{x-x_{0}}\bigg{(}1-\frac{x_{j}-x_{l}}{x-x_{0}}\bigg{(}\dots\bigg{(}1-\frac{x_{j}-x_{u-1}}{x-x_{0}}\bigg{)}\bigg{)}\bigg{)}$$

$$t^{(d,e)}(x)$$

$$=\sum _{j=0}^n\frac{{\zeta }_j^{(d,e)}(x)+{\xi }_j^{(d)}+{\eta }_j^{(d,e)}(x)}{x-x_j}{y}_j$$

$$r^{(d,e)}(x)$$

$$=\sum _{j=0}^n{\beta }_j^{(d,e)}(x)y_j$$

$$\ell (x)=(-1{)}^{n-d+e}(x-x_0{)}^{e+1}(x-x_1)\dots (x-x_{n-1})(x-x_n{)}^{e+1}$$

$${\mu }_i(x)=\prod _{j=-e}^{i-1}(x-x_j)\prod _{k=i+d+1}^{n+e}(x_k-x).$$

$$s^{(d,e)}(x)=\sum _{i\in I}{\mu }_i(x).$$

$$s^{(d,e)}(x_0)=\sum _{i\in I}{\mu }_i(x_0)={\mu }_{-e}(x_0)>0.$$

$$r^{(d,e)}(x_{\alpha })=\frac{{\sum }_{i\in J_{\alpha }}{\mu }_i(x_{\alpha })q_i(x_{\alpha })}{{\sum }_{i\in J_{\alpha }}^{}{\mu }_i(x_{\alpha })}=y_{\alpha }\frac{{\sum }_{i\in J_{\alpha }}{\mu }_i(x_{\alpha })}{{\sum }_{i\in J_{\alpha }}^{}{\mu }_i(x_{\alpha })}=y_{\alpha }$$

$$r^{(d,e)}(x)=\frac{{\sum }_{i\in I}{\mu }_i(x)q_i(x)}{{\sum }_{i\in I}^{}{\mu }_i(x)}=f(x)\frac{{\sum }_{i\in I}{\mu }_i(x)}{{\sum }_{i\in I}^{}{\mu }_i(x)}=f(x).$$

Full text

1 Introduction

Let $f:[a,b]\rightarrow\mathbb{R}$ be an unknown function to be interpolated, $x_{0},x_{1}\dots x_{n}$ be strictly increasing values in the interval $[a,b]$, and $y_{0},y_{1}\dots y_{n}$ be measured function values at these points. Depending on the differentiability of $f$ and the distribution of the $x_{i}$, any one of a number of interpolants in the literature could be used. For example, if the distribution of the $x_{i}$ is not fixed, then the interpolating polynomial of minimum degree is accurate and stable when the $x_{i}$ are chosen as one of the various kinds of Chebyshev points 2013trefethen . If some of the $y_{i}$ are outliers and a high degree of differentiability is not required, then spline interpolation is efficient and confines the effect of the outliers to short subintervals 2001deboor .

Rational interpolants (formed by the ratio of two polynomials) offer an intriguing alternative to polynomials and piecewise polynomials like the ones above. Since the space of rational functions contains that of polynomials, rational interpolants can accurately approximate more diverse function behaviors. The theory of rational interpolants is less developed than that of polynomials though, and some of the available constructions suffer from the appearance of spurious poles on the real line and unattainable points. This is inconvenient enough to have contributed to the historically limited use of rational interpolation in practice.

Among the rational interpolants, the Floater–Hormann family 2007floater is notable for a provable absence of unattainable points and poles along the real line, high rates of approximation, and a simple construction. Let $p_{i,j}(x)$ be the unique polynomial of minimum degree that passes through the points $(x_{i},y_{i})\dots(x_{j},y_{j})$ for $i<j$, and $\chi_{i,j}(x)$ be defined as

[TABLE]

Then for any integer $0\leq d\leq n$, the Floater–Hormann interpolant of degree $d$ is a blend of polynomial interpolants through successive sets of $d+1$ points:

[TABLE]

This form is not often used for computations though, since the barycentric form is more computationally efficient and numerically stable 2004berrut . The present derivation of the barycentric form closely follows that of Floater and Hormann 2007floater , but introduces some notation that will be useful in subsequent sections.

The derivation of the barycentric form begins by writing the polynomials $p_{i,j}(x)$ in the Lagrange form:

[TABLE]

Let $t^{(d)}(x)$ be the numerator of $r^{(d)}(x)$ in Equation 2. Substituting the definitions of $\chi_{i,j}(x)$ and $p_{i,j}(x)$ from Equations 1 and 3 and canceling common factors gives

[TABLE]

It will be convenient in the following to introduce a symbol for the barycentric weights, i.e., the constants in the inner summation:

[TABLE]

Exchanging the order of the summations in Equation 4 and defining the Floater–Hormann weights as

[TABLE]

allows the numerator of the interpolant to be written as

[TABLE]

The denominator of the interpolant can be found by requiring that a constant function be interpolated exactly. That is, the denominator is the same as the numerator when all of the $y_{j}$ are equal to one. This gives

[TABLE]

for the barycentric form of the Floater–Hormann interpolant of degree $d$, where $\beta^{(d)}_{j}(x)$ is the $j$th basis function.

The properties of the Floater–Hormann interpolants for equispaced $x_{i}$ are of particular interest since this situation arises often in practice and is difficult for polynomial interpolants. For specificity, consider the Floater–Hormann interpolant with $d=4$ and $n=16$ for equispaced nodes in the interval $[-1,1]$. The basis functions $\beta^{(4)}_{j}(x)$ appear in Figure 1b, and the Lebesgue function $\lambda^{(4)}(x)=\sum_{j}|\beta^{(4)}_{j}(x)|$ in Figure 1c. The Lebesgue function effectively indicates the relative condition number of the interpolant, i.e., the sensitivity to measurement or rounding errors. This reveals one aspect of the Floater–Hormann interpolants that could be improved—the conditioning degrades at the ends of the interval. More specifically, the supremum of the Lebesgue function (the Lebesgue constant) increases exponentially with $d$ 2012bos . Since the rate of approximation also increases with $d$, one of the main concerns when using Floater–Hormann interpolants is finding a value of $d$ that appropriately balances the interpolation error and the rounding error 2012guttel .

It is worthwhile to consider the source of this ill-conditioning. The definition of the Lebesgue function and Figure 1b reveal that it is caused by alternating oscillations in the basis functions at the ends of the interval. The reason that this occurs only at the ends of the interval and not on the interior is most easily seen from Equation 2, where the functions \chi_{i,i+d}(x)\big{/}\sum_{j=0}^{n-d}\chi_{j,j+d}(x) blend the local polynomial interpolants. Figure 1a shows that while the blending functions on the interior decay rapidly enough to damp the oscillations of the local polynomial interpolants, the blending functions at the ends of the interval do not decay at all. This suggests that the source of the ill-conditioning is a deficit of local interpolants at the ends of the interval.

One proposal to improve the conditioning of the Floater–Hormann interpolants extrapolates to points outside of the interval, and uses a Floater–Hormann interpolant on the extended set of points 2013klein . This does resolve the source of the ill-conditioning identified above, but at the cost of introducing instability by the extrapolation process 2017decamargo . There is evidence that even when the extrapolated points are computed in multiple precision arithmetic, the effect of measurement and rounding errors in the initial $y_{i}$ can negate any advantage of this approach 2017decamargo . Moreover, the extrapolation process obscures the dependence of the interpolants on the initial points, makes explicit basis functions difficult to construct, and complicates the study of the Lebesgue function. All of this means that some other procedure to improve the conditioning of the Floater–Hormann interpolants is greatly desired.

2 An Alternative Interpolant

A possible approach to improve the conditioning of the Floater–Hormann interpolants would be to constrain the derivatives at the endpoints, effectively replacing several of the Lagrange interpolants in Equation 2 with Hermite interpolants. Since the original problem does not include any information about the derivatives of $f$, they would need to be estimated from, e.g., one-sided finite difference formulas. Such formulas are themselves based on polynomial interpolants though, and suffer from ill-conditioning when the nodes are equispaced and the degree of the polynomial is high. This could be mitigated by using finite difference formulas derived from polynomial interpolants of degree less than $d$, but the issue remains that any finite difference formulas would complicate the construction of explicit basis functions. Perhaps then the conditioning of the Floater–Hormann interpolants could be improved by directly including some polynomial interpolants of degree less than $d$, with suitably modified blending functions to confine their influence to the ends of the interval.

Let $\phi^{(d)}_{i}(x)$ and $\psi^{(d)}_{i}(x)$ be the modifications of $\chi_{i,j}(x)$ that will be used to blend polynomial interpolants through fewer than $d+1$ points at the lower and upper ends of the interval:

[TABLE]

Further define the three index sets $I_{1}:=\left\{d-e\dots d-1\right\}$, $I_{2}:=\left\{0\dots n-d\right\}$, and $I_{3}:=\left\{n-d+1\dots n-d+e\right\}$. Then for any integers $0\leq d\leq n$ and $0\leq e\leq d$, the proposed interpolant $r^{(d,e)}(x)$ is a Floater–Hormann interpolant of degree $d$ with $e$ additional polynomial interpolants through fewer points at the lower and upper ends of the interval:

[TABLE]

As before, let $t^{(d,e)}(x)$ be the numerator of $r^{(d,e)}(x)$. Using Equation 5 for the numerator of the Floater–Hormann interpolant, substituting the definitions of $\phi^{(d)}_{i}(x)$, $\psi^{(d)}_{i}(x)$ and $p_{i,j}(x)$ from Equations 8, 9 and 3, and canceling common factors gives

[TABLE]

Analogous to the Floater–Hormann weights, define the pair of functions

[TABLE]

[TABLE]

where $l$ and $u$ are the lower and upper indices of the respective summations, and the second forms of $\zeta_{j}^{(d,e)}(x)$ and $\eta_{j}^{(d,e)}(x)$ follow from the application of Horner’s method. Exchanging the order of summations over $i$ and $j$ in Equation 11 allows $t^{(d,e)}(x)$ to be written as

[TABLE]

where the second equality uses the convention that a sum is zero whenever the upper index is less than the lower index. Given this form for the numerator, the denominator of the proposed interpolant is again found by requiring that a constant function be interpolated exactly, i.e., by setting all of the $y_{j}$ equal to one. This gives

[TABLE]

for the form of the proposed interpolant analogous to Equation 6, where $\beta^{(d,e)}_{j}(x)$ is the $j$th basis function.

While Equation 14 mimics the barycentric form of the Floater–Hormann interpolants, multiplying the numerator and denominator of Equation 10 by

[TABLE]

reveals that the proposed interpolant is a rational function with numerator and denominator degrees at most $n+2e$ and $n-d+2e$. This means that the proposed interpolant cannot be written in the barycentric form of Berrut and Mittelmann 1997berrut , where the $x$ dependence occurs only through the factors $(x-x_{j})^{-1}$. That said, the proposed interpolant can be made to resemble the barycentric form of Schneider and Werner 1991schneider by finding the full partial fraction decomposition of every term in Equation 11. This is not difficult to do by means of the residue method, but the result is found to be less numerically stable than Equation 14 and is not discussed further.

One of the outstanding qualities of the Floater–Hormann interpolants is the provable absence of poles on the real line. A slight modification of Floater and Hormann’s theorem 2007floater is enough to show that the proposed interpolants share this property.

Theorem 2.1

For all $0\leq d\leq n$ and $0\leq e\leq d$, the rational interpolant $r^{(d,e)}(x)$ in Equation 10 has no poles in $\mathbb{R}$.

Proof

Multiply the numerator and denominator of $r^{(d,e)}(x)$ in Equation 10 by $\ell(x)$ in Equation 16, and let $s^{(d,e)}(x)$ be the denominator. It is sufficient to show that $s^{(d,e)}(x)>0$ for all $x\in\mathbb{R}$. To that end, define additional nodes $x_{i}=x_{0}$ for $i\in\left\{-e,-e+1,\dots,-1\right\}$ and $x_{i}=x_{n}$ for $i\in\left\{n,n+1,\dots,n+e\right\}$, define the index set $I:=\left\{-e,-e+1,\dots,n+e\right\}$, and define the functions

[TABLE]

Then $s^{(d,e)}(x)$ can be written as

[TABLE]

At this point, the proof that $s^{(d,e)}(x)>0$ for $x\in\mathbb{R}\backslash\left\{x_{0},x_{n}\right\}$ is identical to that of Floater and Hormann 2007floater up to a relabeling of the indices. Their proof does not extend to $x\in\left\{x_{0},x_{n}\right\}$ because of the assumption that all of the $x_{i}$ are distinct.

First consider $x=x_{0}$. Since $e\leq d$, the first factor in $\mu_{-e}(x)$ is $(x_{j}-x)$ for some $j>0$, and $\mu_{-e}(x_{0})>0$. For all other $i>-e$, the first factor in $\mu_{i}(x)$ is $(x-x_{-e})$ and $\mu_{i}(x_{0})=0$. Hence

[TABLE]

The reasoning to show that $s^{(d,e)}(x_{n})>0$ involves the last factors in the $\mu_{i}(x)$ but is otherwise the same. ∎

Corollary 1

For all $0\leq d\leq n$ and $0\leq e\leq d$, the rational interpolant $r^{(d,e)}(x)$ has no unattainable points.

Proof

For any $\alpha\in\left\{0,1,\dots,n\right\}$, define the index set $J_{\alpha}:=\left\{i\in I:\mu_{i}(x_{\alpha})\neq 0\right\}$. Theorem 2.1 requires that $J_{\alpha}$ be nonempty, since otherwise $\sum_{i\in I}\mu_{i}(x_{\alpha})=0$ and $r^{(d,e)}(x)$ would have a pole at $x_{\alpha}$. Let $q_{i}(x)$ be the $i$th local polynomial interpolant in Equation 10, and observe from the definition of $r^{(d,e)}(x)$ and the fact that polynomial interpolants do not have unattainable points that $q_{i}(x_{\alpha})=y_{\alpha}$ whenever $\mu_{i}(x_{\alpha})\neq 0$. Then

[TABLE]

and $r^{(d,e)}(x)$ interpolates the data at $x_{\alpha}$. ∎

Finally, $r^{(d,e)}(x)$ reproduces polynomials of degree at most $d-e$. If $f$ is such a polynomial, then $q_{i}(x)=f(x)$ for all $i\in I$ and

[TABLE]

This suggests that the proposed interpolant be compared with a Floater–Hormann interpolant of degree $d-e$. The counterpart to Figure 1 would then be, e.g., Figure 2 where the behavior of the proposed interpolant with $d=8$ and $e=4$ is shown. Comparing the ends of the interval in the two figures reveals an increase in the number of blending functions $\mu_{i}(x)/\sum_{j\in I}\mu_{j}(x)$ in Figure 2a, the damping of the oscillations of the basis functions $\beta^{(8,4)}_{j}(x)$ in Figure 2b, and a reduction in the Lebesgue function $\lambda^{(8,4)}(x)=\sum_{i}|\beta^{(8,4)}_{i}(x)|$ in Figure 2c. That is, the proposed modification of the Floater–Hormann interpolants had the intended effect of improving the conditioning at the ends of the interval. Moreover, if the interpolation error is $O(h^{\delta+1})$ where $h$ is the node spacing and $\delta$ is the degree of the local polynomial interpolants, then the proposed interpolant should have lower interpolation error on the interior of the interval than the corresponding Floater–Hormann interpolant. This is supported by the numerical results in the following section.

3 Numerical results

The Floater–Hormann interpolants have the property that once the weights are known, the computational complexity to find the value of the interpolant at some $x$ is $O(n)$ 2004berrut . Ideally, any modification of the Floater–Hormann interpolants would have the same property. This is the case for $r^{(d,e)}(x)$, for which the additional computational complexity over a Floater–Hormann interpolant is $O(de)$. This can be seen from the $O(d)$ values of $\zeta_{j}^{(d,e)}(x)$ and $\eta_{j}^{(d,e)}(x)$ in Equation 14, each of which requires $O(e)$ operations to calculate using Equations 12 and 13. Practically, while the additional computational complexity is small only for small values of $e$, this is sufficient for the cases of interest.

This section specifically considers the behavior of $r^{(d,e)}(x)$ as an approximant for three variations of Runge’s function $f(x)=1/(1+x^{2})$ on equispaced nodes. The first uses the interval $[-5,5]$ since this appears often in the literature 2007floater ; 2013klein . The second breaks the bilateral symmetry and uses the interval $[-3,7]$ to better represent general $f$. The third uses the interval $[-5,5]$, but perturbs the $y_{i}$ with Gaussian noise with $\sigma=10^{-8}$ to simulate the measurement error when $r^{(d,e)}(x)$ is used as an interpolant. These examples are not comprehensive and the results of this section certainly do not carry the weight of proof, but they do seem to be representative of the performance of $r^{(d,e)}(x)$ in practice.

Figure 3 considers the behavior of $r^{(d,e)}(x)$ for $n=64$ as a function of $d$ and $e$. The Lebesgue constant is shown in Figure 3a, and increases exponentially with $d$ for any fixed $e$ (this has been proven for $e=0$ 2012bos ) with the exception of the region $e\geq d-5$ where the Lebesgue constant is small and nearly constant. This supports the supposition that the behavior of the Lebesgue constant is dominated by the polynomial interpolants of degree at most $d-e$ at the ends of the interval. The $L^{\infty}$ error for Runge’s function on the interval $[-5,5]$ is shown in Figure 3b, and is small and nearly constant in the region defined by $e\leq d-5$, $e\geq(d-5)/5$ and $e\geq d-28$. While this intersects $e=0$ at a single point, increasing the value of $e$ dramatically expands the useful interval of $d$ and helps to stabilize the behavior of the approximant. The $L^{\infty}$ error for Runge’s function on the interval $[-3,7]$ is shown in Figure 3c, and is small in the region defined by $e\leq d-5$, $e\geq(d-5)/5$ and $e\geq d-15$. This better represents general $f$, and while the useful interval of $d$ is smaller the behavior is essentially the same as for Figure 3b. Finally, the $L^{\infty}$ error for Runge’s function on the interval $[-5,5]$ with Gaussian noise is shown in Figure 3d, and is small and nearly constant in the region defined by $e\leq d-4$ and $e\geq d-8$. The resemblance to the Lebesgue constant in Figure 3a is to be expected, since the Lebegsue constant effectively indicates the sensitivity of the approximant to measurement errors.

One advantage of the Chebyshev and spline interpolants is that they have few adjustable parameters—the absence of trade-offs in the parameter values makes using them a straightforward affair. For the Floater–Hormann interpolants, there is the temptation to increase $d$ to reduce the interpolation error, but this carries the risk of amplifying the rounding error. While there has been some progress on finding the optimal $d$ when approximating analytic functions 2012guttel , this requires knowledge of the closest singularity in the complex plane. The proposed interpolant apparently complicates the situation further by introducing a second adjustable parameter $e$. Observe though that the line $e=d-8$ intersects the regions of small $L^{\infty}$ error in Figures 3b, 3c and 3d, and that most of stability gains with increasing $e$ have already been achieved on the line $e=4$. These intersect at the point $d,e=12,4$, which seems to offer a reasonable balance of accuracy, stability and computational expense for general use. When $r^{(d,e)}(x)$ is used as an approximant and the measurement error is negligible, increasing $d$ to $12\leq d\leq 16$ can improve the approximation rate at the cost of some stability.

The performance of the proposed approximant as a function of $n$ is compared with that of Chebyshev, cubic spline, and Floater–Hormann approximants for Runge’s function in Figure 4. Figures 4a and 4b use the interval $[-5,5]$, Figures 4c and 4d use the interval $[-3,7]$, and Figures 4e and 4f use the interval $[-5,5]$ with Gaussian noise. There are a number of observations to be made from these figures. First, the Chebyshev and cubic spline approximants show the expected exponential and power law convergence. The Floater–Hormann and proposed approximants initially exhibit exponential convergence and afterward power law convergence, with the transition occurring at an $n$ that increases with $d$. Although Platte, Trefethen and Kuijlaars 2011platte have proven that there is no approximant on equispaced nodes that is stable and converges exponentially in the limit $n\rightarrow\infty$, this is not as serious an issue as one might believe—the Floater–Hormann and proposed approximants have errors that are often substantially less than and reach the level of machine precision well before the Chebyshev approximant. One conclusion then is that there can be some situations where the Floater–Hormann and proposed approximants on equispaced nodes are preferable to the Chebyshev approximant, despite the proven properties 2013trefethen of the latter.

Second, the proposed approximant $r^{(d,e)}(x)$ seems more stable than the Floater–Hormann approximant $r^{(d-e)}(x)$ in three respects. First, the approximantion rate of a Floater–Hormann approximant depends on the parity of $n$ 2007floater , as is visible from the oscillations of the $L^{\infty}$ and $L^{1}$ errors in Figures 4a and 4b. The same oscillations are strongly damped for the $L^{\infty}$ error and almost completely absent for the $L^{1}$ error of the proposed approximant. Second, the $L^{\infty}$ error of the Floater–Hormann approximants in the interval of exponential convergence increases exponentially with $d$ for $d>4$, as is visible from the vertical offset of the curves in Figures 4a, 4c and 4e. The $L^{\infty}$ error of the proposed approximants in the same interval instead collapses onto a single curve in Figures 4a and 4e, with the exception of $n<24$ for $d,e=16,4$. Third, the amplification of the random errors in Figures 4e and 4f is smaller for $r^{(d,e)}(x)$ than for $r^{(d-e)}(x)$, sometimes by an order of magnitude. The reason for this is not clear, since Figure 3a shows that $r^{(d,e)}(x)$ and $r^{(d-e)}(x)$ have similar Lebesgue constants.

Third, the proposed approximants appear to have an advantage over the Floater–Hormann approximants with regard to error. The $L^{\infty}$ and $L^{1}$ errors of $r^{(d,e)}(x)$ are often just above and just below the respective errors of $r^{(d-e)}(x)$ in the interval of power law convergence, as is visible in Figures 4a, 4b, 4c and 4d. Since the $L^{\infty}$ error bounds the pointwise error from above, the Floater–Hormann approximants are preferable here. That said, the proposed approximants can have smaller approximation errors by several orders of magnitude in the interval of exponential convergence, and apparently amplify the random errors substantially less as well. This arguably makes the proposed approximants preferable for the case of general $f$ and $n$.

From Figure 4, the question of the optimal $d$ for a Floater–Hormann approximant with a fixed $n$ seems to involve (at least for Runge’s function) finding the smallest $d$ where the interval of exponential convergence includes $n$—further increasing $d$ can displace the entire curve vertically and increase the approximation error. The advantage of the proposed approximant then is that for modest $e$, the value of $d$ can be safely increased and the interval of exponential convergence expanded without adversely affecting the approximation error. While this sometimes results in smaller overall errors, the more significant practical advantage is that this reduces the necessity of adjusting $d$ to find the optimal value. For example, Table 1 reproduces and expands upon a table in Floater and Hormann 2007floater that reports the optimal values of $d$ for Runge’s function with $n$ equispaced nodes on the interval $[-5,5]$. Since $r^{(d)}(x)$ and $r^{(d,e)}(x)$ are used as approximants and the measurement error is minimal, $d,e=\min(14,n),4$ is used for the proposed approximant instead of the more conservative $d,e=12,4$ suggested above. Observe that even without adjusting $d$ and $e$, the $L^{1}$ error of the proposed approximant is less than that of the optimal Floater–Hormann approximant for every $n$ in the table, and the $L^{\infty}$ error is less than or equal to that of the optimal Floater–Hormann approximant for every $n$ with the exception of $n=20$. For $n=80$, the $L^{\infty}$ error of the proposed approximant is more than an order of magnitude smaller. At the very least then, the proposed approximants could be useful when finding the optimal value of $d$ for a Floater–Hormann approximant would be difficult or time consuming.

4 Conclusion

If one is presented with data on equispaced nodes and desires to interpolate further values, the Floater–Hormann family of rational interpolants is a good choice. They are infintely smooth, more stable than the polynomial interpolant of minimum degree, and often more accurate than cubic spline interpolants. That said, the Floater–Hormann interpolants are a family, and one immediately encounters the question of which one to use in practice. There are at least three answers available in the literature. First, Floater and Hormann 2007floater seem to suggest using a small value of $d$ (perhaps $d=3$) and increasing $n$ until the desired accuracy is achieved. This is a conservative approach, and appears to be widely used 2007press . Second, Güttel and Klein 2012guttel describe a procedure whereby $d$ is given as a function of $n$ depending on the analyticity of $f$. This approach is elegant, but is also more complicated and relies on knowledge of $f$ that is not always available. Third, Klein 2013klein proposed extrapolating to points outside of the interval and constructing Floater–Hormann interpolants on the extended set of points. Certain published results suggest that these Extended Floater–Hormann interpolants do not suffer the usual drawbacks from high values of $d$, but others 2017decamargo suggest that the extrapolation process is a source of significant instability. Klein explicitly ignored this source of error in his analysis.

The rational interpolants proposed in Equations 10 and 14 are modifications of the Floater–Hormann interpolants that blend additional local polynomial interpolants at the ends of the interval. This two-parameter family initially appears to make the above question more difficult to answer, but various numerical examples suggest that a narrow interval of parameter values offers a good balance of accuracy, stability and computational expense. While the accuracy is comparable to that of the Floater–Hormann interpolant with optimal $d$, the proposed interpolants achieve this for constant $d$ and $e$. This is envisioned as simplifying the use of the rational interpolants in practical contexts, where the user might not have enough knowledge of the function $f$ to use a more sophisticated alternative.