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.
Let σn be a subdivision of a real interval [a,b]:
a=x0,x1,…,xn=b. |
The function s is a spline function of degree l if s is a function from [a,b] to ℝ such that:
Proof. Let s be a spline function of degree l on σn. For x∈[a,x1], s is a polynomial A of degree less or equal to l, hence on [a,x1], s=A(x)=a0+a1x+⋯+alxl and A is a linear combination of 1,x,…,xl.
For x∈[x1,x2], s is a polynomial B of degree less or equal to l, hence on [x1,x2], s=B(x)=b0+b1x+⋯+blxl. Since s has continuous derivatives up to order l−1, it follows B(j)(x1)−A(j)(x1)=0 for 0≤ j<l, and therefore B(x)−A(x)=α1(x−x1)l, i.e. B(x)=A(x)+α1(x−x1)l, for some α1. Define the function:
q1(x)= |
|
Now s|[a,x2]=a0+a1x+⋯+alxl+α1q1(x).
For x∈[x2,x3], s is a polynomial C of degree less or equal than l, hence on [x2,x3], s=C(x)=c0+c1x+⋯+clxl. Since s has continuous derivatives up to order l−1, it follows C(j)(x2)−B(j)(x2)=0 for 0≤ j<l, and therefore C(x)−B(x)=α2(x−x2)n or C(x)=B(x)+α2(x−x2)n for some α2∈ℝ. Define the function:
q2(x)= |
|
Now s|[a,x3]=a0+a1x+⋯+alxl+α1q1(x)+α2q2(x).
Continuing, define the functions qj for all 1≤ j<n:
qj(x)= |
|
You obtain s|[a,b]=a0+a1x+⋯+alxl+α1q1(x)+⋯+αn−1qn−1(x), meaning that s is a linear combination of n+l independent functions 1,x,…,xl,q1,…,qn−1. It follows that the set of all possible s is a real vector space of dimension n+l.
If you want to interpolate a function f on σn by a spline function s of degree l, then s must satisfy s(xk)=yk=f(xk) 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.
If l is even (l=2m), there are 2m−1 degrees of freedom. The constraints are defined as follows.
The spline command finds the natural spline.
Find the natural spline of degree 3, crossing through the points (0,1), (1,3) and (2,0):
spline([0,1,2],[1,3,0]) |
or:
spline([0,1,2],[1,3,0],x,3) |
|
The first polynomial, −5/4 x3+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)3−15/4 (x−1)2−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:
spline([0,1,2],[1,3,0],x,3,piecewise) |
|
Find the natural spline of degree 4, crossing through the points (0,1), (1,3), (2,0) and (3,−1):
spline([0,1,2,3],[1,3,0,-1],x,4) |
|
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:
spline([0,1,2,3],[1,3,0,-1],[1/2,4/3,9/4],4) |
|
Find the natural spline interpolation of the cosine function with abscissas [0,π/2,3π/2]:
t:=[0,pi/2,3*pi/2]:; spline(t,cos(t)) |
|
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.
The description of spline types below is quoted from GSL documentation.
Note that Steffen interpolation is available only if giac is compiled with GSL version 2.2 or later.
This example is adapted from GSL documentation. To define some synthetic data, enter:
x,y:=makelist(k->k+0.5*sin(k),0,10),makelist(k->k+cos(k^2),0,10):; |
This gives you 11 data points. To interpolate the function at intermediate points, enter:
t:=linspace(min(x),max(x),1000):; z:=interp(x,y,t):; |
Now plot the data and the interpolated values together:
listplot(tran([t,z])); scatterplot(tran([x,y]),display=blue+point_width_2) |
To see the difference between the supported spline types, interpolate in the segment [4,6]. Enter:
t:=linspace(4,6):; c,a,s:=interp(x,y,t,"cubic"),interp(x,y,t,"akima"),interp(x,y,t,"steffen"):; listplot(tran([t,c]),color=blue+quadrant4,legend="cubic"); listplot(tran([t,a]),color=red,legend="akima"); listplot(tran([t,s]),color=magenta,legend="steffen"); |