6.2 Divided Difference Representation of Polynomials

The functions described here manipulate polynomials stored in Newton’s divided-difference representation. The use of divided-differences is described in Abramowitz & Stegun sections 25.1.4 and 25.2.26, and Burden and Faires, chapter 3, and discussed briefly below.

Given a function f(x), an nth degree interpolating polynomial P_{n}(x) can be constructed which agrees with f at n+1 distinct points x_0,x_1,...,x_{n}. This polynomial can be written in a form known as Newton’s divided-difference representation:

P_n(x) = f(x_0) + \sum_(k=1)^n [x_0,x_1,...,x_k] (x-x_0)(x-x_1)...(x-x_(k-1))

where the divided differences [x_0,x_1,...,x_k] are defined in section 25.1.4 of Abramowitz and Stegun. Additionally, it is possible to construct an interpolating polynomial of degree 2n+1 which also matches the first derivatives of f at the points x_0,x_1,...,x_n. This is called the Hermite interpolating polynomial and is defined as

H_(2n+1)(x) = f(z_0) + \sum_(k=1)^(2n+1) [z_0,z_1,...,z_k] (x-z_0)(x-z_1)...(x-z_(k-1))

where the elements of z = \{x_0,x_0,x_1,x_1,...,x_n,x_n\} are defined by z_{2k} = z_{2k+1} = x_k. The divided-differences [z_0,z_1,...,z_k] are discussed in Burden and Faires, section 3.4.

Function: int gsl_poly_dd_init (double dd[], const double xa[], const double ya[], size_t size)

This function computes a divided-difference representation of the interpolating polynomial for the points (x, y) stored in the arrays xa and ya of length size. On output the divided-differences of (xa,ya) are stored in the array dd, also of length size. Using the notation above, dd[k] = [x_0,x_1,...,x_k].

Function: double gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x)

This function evaluates the polynomial stored in divided-difference form in the arrays dd and xa of length size at the point x. An inline version of this function is used when HAVE_INLINE is defined.

Function: int gsl_poly_dd_taylor (double c[], double xp, const double dd[], const double xa[], size_t size, double w[])

This function converts the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation is supplied in the arrays dd and xa of length size. On output the Taylor coefficients of the polynomial expanded about the point xp are stored in the array c also of length size. A workspace of length size must be provided in the array w.

Function: int gsl_poly_dd_hermite_init (double dd[], double za[], const double xa[], const double ya[], const double dya[], const size_t size)

This function computes a divided-difference representation of the interpolating Hermite polynomial for the points (x, y) stored in the arrays xa and ya of length size. Hermite interpolation constructs polynomials which also match first derivatives dy/dx which are provided in the array dya also of length size. The first derivatives can be incorported into the usual divided-difference algorithm by forming a new dataset z = \{x_0,x_0,x_1,x_1,...\}, which is stored in the array za of length 2*size on output. On output the divided-differences of the Hermite representation are stored in the array dd, also of length 2*size. Using the notation above, dd[k] = [z_0,z_1,...,z_k]. The resulting Hermite polynomial can be evaluated by calling gsl_poly_dd_eval and using za for the input argument xa.