Barycentric Rational Interpolation

Synopsis

#include <boost/math/interpolators/barycentric_rational.hpp>

namespace boost{ namespace math{
    template<class Real>
    class barycentric_rational
    {
    public:
        template <class InputIterator1, class InputIterator2>
        barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);

        barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t order = 3);

        barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);

        Real operator()(Real x) const;

        Real prime(Real x) const;

        std::vector<Real>&& return_x();

        std::vector<Real>&& return_y();
    };

}}

Barycentric rational interpolation is a high-accuracy interpolation method for non-uniformly spaced samples. It requires 𝑶(N) time for construction, and 𝑶(N) time for each evaluation. Linear time evaluation is not optimal; for instance the cubic B-spline can be evaluated in constant time. However, using the cubic B-spline requires uniformly-spaced samples, which are not always available.

Use of the class requires a vector of independent variables x[0], x[1], .... x[n-1] where x[i+1] > x[i], and a vector of dependent variables y[0], y[1], ... , y[n-1]. The call is trivial:

std::vector<double> x(500);
std::vector<double> y(500);
// populate x, y, then:
boost::math::barycentric_rational<double> interpolant(std::move(x), std::move(y));

This implicitly calls the constructor with approximation order 3, and hence the accuracy is 𝑶(h⁴). In general, if you require an approximation order d, then the error is 𝑶(h^d+1). A call to the constructor with an explicit approximation order is demonstrated below

boost::math::barycentric_rational<double> interpolant(std::move(x), std::move(y), 5);

To evaluate the interpolant, simply use

double x = 2.3;
double y = interpolant(x);

and to evaluate its derivative use

double y = interpolant.prime(x);

If you no longer require the interpolant, then you can get your data back:

std::vector<double> xs = interpolant.return_x();
std::vector<double> ys = interpolant.return_y();

Be aware that once you return your data, the interpolant is dead.

Caveats

Although this algorithm is robust, it can surprise you. The main way this occurs is if the sample spacing at the endpoints is much larger than the spacing in the center. This is to be expected; all interpolants perform better in the opposite regime, where samples are clustered at the endpoints and somewhat uniformly spaced throughout the center.

A desirable property of any interpolator f is that for all x_min ≤ x ≤ x_max, also y_min ≤ f(x) ≤ y_max.

This property does not hold for the barycentric rational interpolator. However, unless you deliberately try to antagonize this interpolator (by, for instance, placing the final value far from all the rest), you will probably not fall victim to this problem.

The reference used for implementation of this algorithm is Barycentric rational interpolation with no poles and a high rate of interpolation, and the evaluation of the derivative is given by Some New Aspects of Rational Interpolation.

Examples

This example shows how to use barycentric rational interpolation, using Walter Kohn's classic paper "Solution of the Schrodinger Equation in Periodic Lattices with an Application to Metallic Lithium" In this paper, Kohn needs to repeatedly solve an ODE (the radial Schrodinger equation) given a potential which is only known at non-equally samples data.

If he'd only had the barycentric rational interpolant of Boost.Math!

References: Kohn, W., and N. Rostoker. "Solution of the Schrodinger equation in periodic lattices with an application to metallic lithium." Physical Review 94.5 (1954): 1111.

#include <boost/math/interpolators/barycentric_rational.hpp>

int main()
{
    // The lithium potential is given in Kohn's paper, Table I:
    std::vector<double> r(45);
    std::vector<double> mrV(45);

    // We'll skip the code for filling the above vectors with data for now...

    // Now we want to interpolate this potential at any r:
    boost::math::barycentric_rational<double> b(r.data(), mrV.data(), r.size());

    for (size_t i = 1; i < 8; ++i)
    {
        double r = i*0.5;
        std::cout <<  "(r, V) = (" << r << ", " << -b(r)/r << ")\n";
    }
}

This further example shows how to use the iterator based constructor, and then uses the function object in our root finding algorithms to locate the points where the potential achieves a specific value.

#include <boost/math/interpolators/barycentric_rational.hpp>
#include <boost/range/adaptors.hpp>
#include <boost/math/tools/roots.hpp>

int main()
{
    // The lithium potential is given in Kohn's paper, Table I.
    // (We could equally easily use an unordered_map, a list of tuples or pairs, or a 2-dimensional array).
    std::map<double, double> r;

    r[0.02] = 5.727;
    r[0.04] = 5.544;
    r[0.06] = 5.450;
    r[0.08] = 5.351;
    r[0.10] = 5.253;
    r[0.12] = 5.157;
    r[0.14] = 5.058;
    r[0.16] = 4.960;
    r[0.18] = 4.862;
    r[0.20] = 4.762;
    r[0.24] = 4.563;
    r[0.28] = 4.360;
    r[0.32] = 4.1584;
    r[0.36] = 3.9463;
    r[0.40] = 3.7360;
    r[0.44] = 3.5429;
    r[0.48] = 3.3797;
    r[0.52] = 3.2417;
    r[0.56] = 3.1209;
    r[0.60] = 3.0138;
    r[0.68] = 2.8342;
    r[0.76] = 2.6881;
    r[0.84] = 2.5662;
    r[0.92] = 2.4242;
    r[1.00] = 2.3766;
    r[1.08] = 2.3058;
    r[1.16] = 2.2458;
    r[1.24] = 2.2035;
    r[1.32] = 2.1661;
    r[1.40] = 2.1350;
    r[1.48] = 2.1090;
    r[1.64] = 2.0697;
    r[1.80] = 2.0466;
    r[1.96] = 2.0325;
    r[2.12] = 2.0288;
    r[2.28] = 2.0292;
    r[2.44] = 2.0228;
    r[2.60] = 2.0124;
    r[2.76] = 2.0065;
    r[2.92] = 2.0031;
    r[3.08] = 2.0015;
    r[3.24] = 2.0008;
    r[3.40] = 2.0004;
    r[3.56] = 2.0002;
    r[3.72] = 2.0001;

    // Let's discover the absissa that will generate a potential of exactly 3.0,
    // start by creating 2 ranges for the x and y values:
    auto x_range = boost::adaptors::keys(r);
    auto y_range = boost::adaptors::values(r);
    boost::math::barycentric_rational<double> b(x_range.begin(), x_range.end(), y_range.begin());
    //
    // We'll use a lambda expression to provide the functor to our root finder, since we want
    // the abscissa value that yields 3, not zero.  We pass the functor b by value to the
    // lambda expression since barycentric_rational is trivial to copy.
    // Here we're using simple bisection to find the root:
    boost::uintmax_t iterations = (std::numeric_limits<boost::uintmax_t>::max)();
    double abscissa_3 = boost::math::tools::bisect([=](double x) { return b(x) - 3; }, 0.44, 1.24, boost::math::tools::eps_tolerance<double>(), iterations).first;
    std::cout << "Abscissa value that yields a potential of 3 = " << abscissa_3 << std::endl;
    std::cout << "Root was found in " << iterations << " iterations." << std::endl;
    //
    // However, we have a more efficient root finding algorithm than simple bisection:
    iterations = (std::numeric_limits<boost::uintmax_t>::max)();
    abscissa_3 = boost::math::tools::bracket_and_solve_root([=](double x) { return b(x) - 3; }, 0.6, 1.2, false, boost::math::tools::eps_tolerance<double>(), iterations).first;
    std::cout << "Abscissa value that yields a potential of 3 = " << abscissa_3 << std::endl;
    std::cout << "Root was found in " << iterations << " iterations." << std::endl;
}

Program output is:

Abscissa value that yields a potential of 3 = 0.604728
Root was found in 54 iterations.
Abscissa value that yields a potential of 3 = 0.604728
Root was found in 10 iterations.