Barycentric Series

Interpolating polynomials can be expressed in a barycentric form:

where x_i are the points (nodes or support points) where the function is being interpolated, w_i are the corresponding weights, and f_i are the function values at the support points.

The converse is not true: not every barycentric form is equivalent to a polynomial. In general, a barycentric form represents a rational function: the quotient of two polynomials.

The barycentric form has better numerical properties than the commonly used polynomial form. This is an important reason why its use is preferred over the polynomial forms.

A barycentric basis is defined by a set of support points or nodes and corresponding weights. Each of the basis functions has the value 1 at one of the support points and 0 at all the other support points. This makes it super easy to construct the interpolating curve through a set of points.

In the Extreme Optimization Numerical Libraries for .NET , barycentric basis are represented by the BarycentricBasis class.

The BarycentricBasis class has one constructor. It takes two arguments: a vector containing the support points, and a vector containing the corresponding weights.

In most cases, however, you will want to use one of the static methods of the class to construct a basis. These are listed below.

As mentioned earlier, for a given set of points, specific values for the weights give rise to a basis that consists of polynomials. This is called a Lagrange basis. There are two static methods that can be used to define a polynomial basis.

The Polynomial method constructs a polynomial basis for a set of points. It has two overloads. The first overload takes one argument: a VectorT that contains the support points. The method returns a polynomial basis that interpolates through those points.

The second overload creates a basis for a set of points that are equally spaced on an interval. It takes three arguments. the lower bound and upper bound of the interval, and the number of support points.

The example below creates two barycentric bases. The first is for five points over the interval [0, 2]. The second is for five specific X-values:

// Polynomial basis of equidistant points:
BarycentricBasis pb1 = BarycentricBasis.Polynomial(0.0, 2.0, 5);
// Polynomial basis for specific interpolation points:
var x1 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0);
BarycentricBasis pb2 = BarycentricBasis.Polynomial(x1);

' Polynomial basis of equidistant points:
Dim pb1 As BarycentricBasis = BarycentricBasis.Polynomial(0.0, 2.0, 5)
' Polynomial basis for specific interpolation points:
Dim x1 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0)
Dim pb2 As BarycentricBasis = BarycentricBasis.Polynomial(x1)

// Polynomial basis of equidistant points:
let pb1 = BarycentricBasis.Polynomial(0.0, 2.0, 5);
// Polynomial basis for specific interpolation points:
let x1 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0);
let pb2 = BarycentricBasis.Polynomial(x1);

It should be noted that polynomial interpolation is ill-conditioned in general. This is especially true when the interpolation points are roughly evenly spaced. It is therefore only recommended for relatively low degree.

As mentioned elsehwere, Chebyshev polynomials form a numerically much more desirable basis for polynomials than regular polynomials, when handled properly. Interpolation through so-called Chebyshev points can be done in a numerically stable way.

Standard interpolation using Chebyshev polynomials is done over the interval [-1,1]. In that case, the roots of a Chebyshev polynomial of degree n form a set of n Chebyshev points. For different intervals, a simple linear transformation preserves the desirable properties of a Chebyshev basis.

In barycentric form, the weights for a Chebyshev basis take on a very simple form. The weights are all relatively the same size, which is an indication of numerical stability.

Two methods are available to construct a Chebyshev basis in barycentric form. The Chebyshev method constructs a Chebyshev basis for Chebyshev polynomials of the first kind. The method has two overloads. The first overload simply takes the number of support points and produces a basis over the interval [-1,1]. The second overload takes three arguments: the lower and upper bound of the interpolation interval and the number of support points.

The Chebyshev2 method constructs a Chebyshev basis for Chebyshev polynomials of the second kind. The method again has two overloads with the same meaning.

In the example below, we first create a Chebyshev basis of degree 5 over the standard interval [-1,1], followed by a Chebyshev basis of degree 50 over the interval [0,10]. Finally, we create a Chebyshev basis of the second kind of degree 20 over [-1,1].

// Chebyshev basis of degree 5 over [-1,1]:
BarycentricBasis cb1 = BarycentricBasis.Chebyshev(5);
// Chebyshev basis of degree 50 over [0,10]:
BarycentricBasis cb2 = BarycentricBasis.Chebyshev(0.0, 10.0, 50);
// Chebyshev basis (2nd kind) of degree 20 over [-1,1]:
BarycentricBasis cb3 = BarycentricBasis.Chebyshev2(20);

' Chebyshev basis of degree 5 over [-1,1]:
Dim cb1 As BarycentricBasis = BarycentricBasis.Chebyshev(5)
' Chebyshev basis of degree 50 over [0,10]:
Dim cb2 As BarycentricBasis = BarycentricBasis.Chebyshev(0.0, 10.0, 50)
' Chebyshev basis (2nd kind) of degree 20 over [-1,1]:
Dim cb3 As BarycentricBasis = BarycentricBasis.Chebyshev2(20)

// Chebyshev basis of degree 5 over [-1,1]:
let cb1 = BarycentricBasis.Chebyshev(5);
// Chebyshev basis of degree 50 over [0,10]:
let cb2 = BarycentricBasis.Chebyshev(0.0, 10.0, 50);
// Chebyshev basis (2nd kind) of degree 20 over [-1,1]:
let cb3 = BarycentricBasis.Chebyshev2(20);

Up to this point, we have only considered bases that produce polynomials. A barycentric basis can represent a class of rational functions as well. One particularly interesting kind is a Floater-Hormann basis.

A Floater-Hormann basis is characterised by a set of support points and an integer, often called the order or the blending parameter. It works by blending interpolating polynomials over a sliding window of support points. The result is a set of rational functions that has desirable numerical properties as long as the blending parameter is not too large. Values of 3 or 4 seem to work well in many cases.

The FloaterHormann method constructs a Floater-Hormann basis. It has two overloads. The first overload takes two arguments. The first is a VectorT that contains the support points in ascending order. The second argument is the order or bleding parameter. The method returns a basis of rational functions that interpolate through the support points.

The second overload creates a basis for a set of points that are equally spaced on an interval. It takes four arguments. the lower bound and upper bound of the interval, the number of support points, and finally the order.

// Floater-Hormann basis for 11 equally spaced points
// over [-5,5]. The order of the basis is 3.
BarycentricBasis fhb1 = BarycentricBasis.FloaterHormann(-5.0, 5.0, 12, 3);
// Floater-Hormann basis of order 2 for specific
// support points.
var x2 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0);
BarycentricBasis fhb2 = BarycentricBasis.FloaterHormann(x2, 2);

' Floater-Hormann basis for 11 equally spaced points
' over [-5,5]. The order of the basis Is 3.
Dim fhb1 As BarycentricBasis = BarycentricBasis.FloaterHormann(-5.0, 5.0, 12, 3)
' Floater-Hormann basis of order 2 for specific
' support points.
Dim x2 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0)
Dim fhb2 As BarycentricBasis = BarycentricBasis.FloaterHormann(x2, 2)

// Floater-Hormann basis for 11 equally spaced points
// over [-5,5]. The order of the basis is 3.
let fhb1 = BarycentricBasis.FloaterHormann(-5.0, 5.0, 12, 3);
// Floater-Hormann basis of order 2 for specific
// support points.
let x2 = Vector.Create(1.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0);
let fhb2 = BarycentricBasis.FloaterHormann(x2, 2);

A curve that represents a specific curve from the set defined by a barycentric basis is called a barycentric series. Barycentric series are represented by the BarycentricSeries class.

Because barycentric curves are linear combinations of functions from the basis, BarycentricSeries inherits from the LinearCombination class.

The BarycentricSeries class has one constructor that takes two arguments. The first is a BarycentricBasis. The second is a vector containing the values of the curve at the support points of the basis. This creates a rational curve (or polynomial) that interpolates through the support points.

In the example that follows, we first create a polynomial basis. Then we create a vector with function values. Finally, we construct a barycentric curve using the basis and vector of function values.

var basis = BarycentricBasis.Polynomial(0.0, 2.0, 5);
var f = Vector.Create(1.0, -2.0, 3.0, -4.0, 5.0);
var curve = new BarycentricSeries(basis, f);

Dim basis = BarycentricBasis.Polynomial(0.0, 2.0, 5)
Dim f = Vector.Create(1.0, -2.0, 3.0, -4.0, 5.0)
Dim c = New BarycentricSeries(basis, f)

let basis = BarycentricBasis.Polynomial(0.0, 2.0, 5);
let f = Vector.Create(1.0, -2.0, 3.0, -4.0, 5.0);
let curve = new BarycentricSeries(basis, f);

In addition, the BarycentricSeries class provides a number of static methods that may be used to construct interpolating curves directly.

The GetPolynomialInterpolator constructs the Lagrange interpolating polynomial through a set of points in barycentric form. The method takes two vector arguments that contain the X and Y coordinates of the interpolation points. This method is functionally equivalent to GetInterpolatingPolynomial, but with the benefits in terms of numerical stability of the barycentric form over traditional polynomial form.

The GetChebyshevInterpolator method constructs a curve that interpolates a function through the Chebyshev points over an interval. The method has four arguments. The first is a delegate that evaluates the function that should be interpolated. The second and third arguments are the lower and upper bounds of the interpolation interval. The fourth argument is the degree of the resulting Chebyshev polynomial. The GetChebyshev2Interpolator method is similar but uses the Chebyshev points of the second kind as support points.

The GetFloaterHormannInterpolator method constructs a rational curve that interpolates a function through a set of points. The method has two overloads that mirror the methods that create a Floater-Hormann basis. The first overload takes three arguments. The first two are vectors that contain the X and Y coordinates of the interpolation points. The third argument is the order of the Floater-Hormann basis.

The second overload takes a total of five arguments. The first is a delegate that evaluates the function to interpolate. The second and third arguments are the lower and uppoer bounds of the interpolation interval. The fourth argument is the number of points to use for the construction of the interpolating function. The last argument is again the order of the Floater-Hormann basis.

The example below constructs interpolating barycentric curves using each of the methods listed above:

// x and y values contain the data points:
const int n = 13;
var xValues = Vector.CreateRange(0.0, 2.0, n);
var yValues = Vector.Sin(xValues);
// Polynomial interpolant for equidistant points:
var poly = BarycentricSeries.GetPolynomialInterpolator(xValues, yValues);
// Chebyshev interpolant (1st kind):
var cheb1 = BarycentricSeries.GetChebyshevInterpolator(Math.Sin, 0.0, 2.0, n);
// Chebyshev interpolant (2nd kind):
var cheb2 = BarycentricSeries.GetChebyshev2Interpolator(Math.Sin, 0.0, 2.0, n);
// Two ways to create Floater-Hormann interpolant:
// - using equidistant points (order 2)
var fh2 = BarycentricSeries.GetFloaterHormannInterpolator(Math.Sin, 0.0, 2.0, n, 2);
// - using supplied support points (order 4)
var fh4 = BarycentricSeries.GetFloaterHormannInterpolator(xValues, yValues, 2);

' x And y values contain the data points
Const n As Integer = 13
Dim xValues = Vector.CreateRange(0.0, 2.0, n)
Dim yValues = Vector.Sin(xValues)
' Polynomial interpolant for equidistant points:
Dim poly = BarycentricSeries.GetPolynomialInterpolator(xValues, yValues)
' Chebyshev interpolant (1st kind):
Dim cheb1 = BarycentricSeries.GetChebyshevInterpolator(AddressOf Math.Sin, 0.0, 2.0, n)
' Chebyshev interpolant (2nd kind):
Dim cheb2 = BarycentricSeries.GetChebyshev2Interpolator(AddressOf Math.Sin, 0.0, 2.0, n)
' Two ways to create Floater-Hormann interpolant:
' - using equidistant points (order 2)
Dim fh2 = BarycentricSeries.GetFloaterHormannInterpolator(AddressOf Math.Sin, 0.0, 2.0, n, 2)
' - using supplied support points (order 4)
Dim fh4 = BarycentricSeries.GetFloaterHormannInterpolator(xValues, yValues, 2)

// x and y values contain the data points:
let n = 13;
let xValues = Vector.CreateRange(0.0, 2.0, n);
let yValues = Vector.Sin(xValues);
// Polynomial interpolant for equidistant points:
let poly = BarycentricSeries.GetPolynomialInterpolator(xValues, yValues);
// Chebyshev interpolant (1st kind):
let cheb1 = BarycentricSeries.GetChebyshevInterpolator(func Math.Sin, 0.0, 2.0, n);
// Chebyshev interpolant (2nd kind):
let cheb2 = BarycentricSeries.GetChebyshev2Interpolator(func Math.Sin, 0.0, 2.0, n);
// Two ways to create Floater-Hormann interpolant:
// - using equidistant points (order 2)
let fh2 = BarycentricSeries.GetFloaterHormannInterpolator(func Math.Sin, 0.0, 2.0, n, 2);
// - using supplied support points (order 4)
let fh4 = BarycentricSeries.GetFloaterHormannInterpolator(xValues, yValues, 2);

The BarycentricSeries class implements nearly all methods and properties of the Curve class. The one exception is integration (quadrature). Because there is no closed-form formula for the integral of a rational curve, numerical integration is used instead. Because most barycentric curves used in practice are well-behaved, this is not really an issue.

The ValueAt method returns the value of the polynomial at a specified point. SlopeAt returns the derivative.

var b1 = BarycentricSeries.GetChebyshevInterpolator(Math.Tan, -1, 1, 12);
Console.WriteLine("b1.ValueAt(0.5) = {0}", b1.ValueAt(0.5));
Console.WriteLine("b1.SlopeAt(0.5) = {0}", b1.SlopeAt(0.5));

Dim b1 = BarycentricSeries.GetChebyshevInterpolator(AddressOf Math.Tan, -1, 1, 12)
Console.WriteLine("b1.ValueAt(0.5) = {0}", b1.ValueAt(0.5))
Console.WriteLine("b1.SlopeAt(0.5) = {0}", b1.SlopeAt(0.5))

let b1 = BarycentricSeries.GetChebyshevInterpolator(func Math.Tan, -1, 1, 12);
Console.WriteLine("b1.ValueAt(0.5) = {0}", b1.ValueAt(0.5));
Console.WriteLine("b1.SlopeAt(0.5) = {0}", b1.SlopeAt(0.5));

Integral evaluates the definite integral over a specified interval. As mentioned above, the integral is always evaluated numerically.