17.1.2  Spline interpolation

Spline interpolation is a method of piecewise polynomial interpolation in which smaller, low-degree pieces are glued together at given points under suitable boundary conditions. With Xcas you can do both symbolic and numerical spline interpolation.

Theoretical background.

Let σ_n be a subdivision of a real interval [a,b]:

The function s is a spline function of degree l if s is a function from [a,b] to ℝ such that:

• s has continuous derivatives up to the order l−1,
• on each interval of the subdivision σ_n, s is a polynomial of degree less or equal than l.

Proof. Let s be a spline function of degree l on σ_n. For x∈[a,x₁], s is a polynomial A of degree less or equal to l, hence on [a,x₁], s=A(x)=a₀+a₁x+⋯+a_lx^l and A is a linear combination of 1,x,…,x^l.

For x∈[x₁,x₂], s is a polynomial B of degree less or equal to l, hence on [x₁,x₂], s=B(x)=b₀+b₁x+⋯+b_lx^l. Since s has continuous derivatives up to order l−1, it follows B^(j)(x₁)−A^(j)(x₁)=0 for 0≤ j<l, and therefore B(x)−A(x)=α₁(x−x₁)^l, i.e. B(x)=A(x)+α₁(x−x₁)^l, for some α₁. Define the function:

Now s|_[a,x2]=a₀+a₁x+⋯+a_lx^l+α₁q₁(x).

For x∈[x₂,x₃], s is a polynomial C of degree less or equal than l, hence on [x₂,x₃], s=C(x)=c₀+c₁x+⋯+c_lx^l. Since s has continuous derivatives up to order l−1, it follows C^(j)(x₂)−B^(j)(x₂)=0 for 0≤ j<l, and therefore C(x)−B(x)=α₂(x−x₂)ⁿ or C(x)=B(x)+α₂(x−x₂)ⁿ for some α₂∈ℝ. Define the function:

Now s|_[a,x3]=a₀+a₁x+⋯+a_lx^l+α₁q₁(x)+α₂q₂(x).

Continuing, define the functions q_j for all 1≤ j<n:

You obtain s|_[a,b]=a₀+a₁x+⋯+a_lx^l+α₁q₁(x)+⋯+α_n−1q_n−1(x), meaning that s is a linear combination of n+l independent functions 1,x,…,x^l,q₁,…,q_n−1. It follows that the set of all possible s is a real vector space of dimension n+l.

Types of spline functions.

If you want to interpolate a function f on σ_n by a spline function s of degree l, then s must satisfy s(x_k)=y_k=f(x_k) for all 0≤ k≤ n. This gives n+1 conditions, leaving l−1 degrees of freedom. You can therefore add l−1 conditions, these conditions are on the derivatives of s at a and b.

Hermite interpolation, natural interpolation and periodic interpolation are three kinds of interpolation obtained by specifying three kinds of constraints. The uniqueness of the solution of the interpolation problem can be proved for each kind of constraints.

If l is odd (i.e. l=2m−1), there are 2m−2 degrees of freedom. The constraints are defined as follows.

• Hermite interpolation: s^(j)(a)=f^(j)(a) and s^(j)(b)=f^(j)(b), 1≤ j≤ m−1.
• Natural interpolation: s^(j)(a)=s^(j)(b)=0, m ≤ j ≤ 2m−2.
• Periodic interpolation: s^(j)(a)=s^(j)(b), 1≤ j≤ 2m−2.

If l is even (l=2m), there are 2m−1 degrees of freedom. The constraints are defined as follows.

• Hermite interpolation: s^(j)(a)=f^(j)(a), s^(j)(b)=f^(j)(b) and s^(m)(a)=f^(m)(a), 1≤ j≤ m−1.
• Natural interpolation: s^(j)(a)=s^(j)(b)=0 and s^(2m−1)(a)=0, m ≤ j ≤ 2m−2.
• Periodic interpolation: s^(j)(a)=s^(j)(b), 1≤ j≤ 2m−1.

Natural symbolic interpolation.

The spline command finds the natural spline.

• spline takes two mandatory arguments and up to three optional arguments:
  • L_x, a list of abscissas (in increasing order).
  • L_y, a list of ordinates (the same length as L_x).
  • Optionally, x, either a symbolic variable, real number, or list of reals (by default x).
  • Optionally, l, a positive integer for the degree (by default, l=3).
  • Optionally, piecewise, the symbol.
• If x is a variable, then spline(L_x,L_y ⟨,x,l ⟩ ⟨,piecewise ⟩) returns the natural spline function s of degree l, where s(L_x,j)=L_y,j for j=0,…,length(L_x), as a list of polynomials, each polynomial being valid on the corresponding interval determined by L_x, or, if piecewise is given, a piecewise expression representing the entire spline for minL_x≤ x≤ maxL_x.
• If x is a real number (resp. a list of numbers), then spline(L_x,L_y,x ⟨,l ⟩) returns the spline value at x (resp. the list of spline values at individual components of x).

Examples

Find the natural spline of degree 3, crossing through the points (0,1), (1,3) and (2,0):

or:

The first polynomial, −5/4 x³+13/4 x+1, is defined on the interval [0,1] (the first interval defined by the list [0,1,2]) and the second polynomial 5/4 (x−1)³−15/4 (x−1)²−x−1/2+3 is defined on the interval [1,2], the second interval defined by the list [0,1,2]. To obtain the result in a piecewise form, enter:

Find the natural spline of degree 4, crossing through the points (0,1), (1,3), (2,0) and (3,−1):

The output is a list of three polynomial functions of x, defined on the intervals [0,1], [1,2] and [2,3], respectively. To compute the spline value for points x=1/2,4/3,9/4, enter:

Find the natural spline interpolation of the cosine function with abscissas [0,π/2,3π/2]:

Numerical spline interpolation.

The interp command interpolates a discretized function at a list of given points by using fast spline interpolation routines provided by GSL. You should use this command when you do not need the interpolated function itself but only its values at a sequence of points.

• interp takes two or three mandatory arguments and a sequence of optional arguments:
  • data, a sequence or list of two vectors x and y of the same length.
  • v, a vector of x-values at which the spline is evaluated.
  • Optionally, opts, a sequence of options each of which is one of:
    • "stype", a string specifying spline type. Allowed values for stype are cubic, akima and steffen (by default, stype=cubic). See below for descriptions of these types.
    • periodic, the symbol specifying that interpolation should be periodic (by default, interpolation is not periodic). This argument gets ignored if stype=steffen.
    • d, a nonnegative integer specifying the order of the derivative of interpolated function to be evaluated. Allowed values are 0, 1 and 2 (by default, d=0).
• interp([x,y],v ⟨,opts ⟩) or interp(x,y,v ⟨,opts ⟩) returns the list z of spline (derivative) values at points in v.
• Note that data and v must be entirely numeric and data should also be inexact.
• The elements of x must be strictly ascending. The minimal number of data points is 3.
• The elements of t must be between minx and maxx, otherwise the output list will contain undefined elements (as an effect, plotting would fail).

The description of spline types below is quoted from GSL documentation.

cubic

Cubic spline with natural boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The second derivative is chosen to be zero at the first point and last point.

periodic cubic

Cubic spline with periodic boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The derivatives at the first and last points are also matched. Note that the last point in the data must have the same y-value as the first point, otherwise the resulting periodic interpolation will have a discontinuity at the boundary.

(periodic) akima

Non-rounded Akima spline with natural (periodic) boundary conditions. This method uses the non-rounded corner algorithm of Wodicka.

steffen

Steffen’s method guarantees the monotonicity of the interpolating function between the given data points. Therefore, minima and maxima can only occur exactly at the data points, and there can never be spurious oscillations between data points. The interpolated function is piecewise cubic in each interval. The resulting curve and its first derivative are guaranteed to be continuous, but the second derivative may be discontinuous.

Note that Steffen interpolation is available only if giac is compiled with GSL version 2.2 or later.

Example

This example is adapted from GSL documentation. To define some synthetic data, enter:

This gives you 11 data points. To interpolate the function at intermediate points, enter:

Now plot the data and the interpolated values together:

To see the difference between the supported spline types, interpolate in the segment [4,6]. Enter: