Mathematics 171 – Numerical Analysis 1 Lab Day 4 – Cubic Spline Interpolation November 5, 1998

Goals and Background

The goals of today’s lab are to introduce the Maple spline procedure for free cubic splines and another “homegrown” Maple procedure for clamped cubic splines. We will use them to construct cubic spline interpolating functions. Recall that on Tuesday we saw that given x_{0}, …, x_{n} and function values f(x_{0}), …, f(x_{n}) to be interpolated, the coefficients a_{j}, b_{j}, c_{j}, d_{j} in the free cubic spline interpolating function

S(x) = { S_{0}(x) = a_{0} + b_{0}(x − x_{0}) + c_{0}(x − x_{0})^{2} + d_{0}(x − x_{0})^{3}, x_{0} ≤ x ≤ x_{1} ; S_{1}(x) = a_{1} + b_{1}(x − x_{1}) + c_{1}(x − x_{1})^{2} + d_{1}(x − x_{1})^{3}, x_{1} ≤ x ≤ x_{2} ; ⋮ ; S_{n−1}(x) = a_{n−1} + b_{n−1}(x − x_{n−1}) + c_{n−1}(x − x_{n−1})^{2} + d_{n−1}(x − x_{n−1})^{3}, ; x_{n−1} ≤ x ≤ x_{n}

satisfy a_{j} = f(x_{j}) for 0 ≤ j ≤ n − 1 and

\begin{aligned}b_{j} & = \frac{a_{j+1} − a_{j}}{h_{j}} − \frac{h_{j}}{3} (2c_{j} + c_{j+1}) \\ d_{j} & = \frac{c_{j+1} − c_{j}}{3h_{j}} \end{aligned}

for all j = 0, …, n. The c_{j} are found by solving the (n − 1) × (n − 1) system of linear equations:

\begin{matrix}( \begin{matrix}2(h_{0} + h_{1}) & h_{1} & 0 & ⋯ & 0 & 0 \\ h_{1} & 2(h_{1} + h_{2}) & h_{2} & ⋯ & 0 & 0 \\ ⋮ & & & & ⋮ & ⋮ \\ 0 & ⋯ & 0 & 0 & h_{n−2} & 2(h_{n−2} + h_{n−1}) \end{matrix} ) · ( c_{1} ; ⋮ ; c_{n−1} ) \\ = ( \frac{3(a_{2}−a_{1})}{h_{1}} − \frac{3(a_{1}−a_{0})}{h_{0}} ; ⋮ ; \frac{3(a_{n}−a_{n−1})}{h_{n−1}} − \frac{3(a_{n−2}−a_{n−1})}{h_{n−1}} ) \end{matrix}

where h_{j} = x_{j+1} − x_{j} and c_{0} = c_{n} = 0. There is a very similiar system of linear equations whose solutions give the c_{j} coefficients for clamped cubic splines.

A New Maple Command

Maple has a spline interpolation function that can be used to construct free cubic spline interpolating functions. It is invoked by first entering the command:

                         readlib(spline);

(This does something very close to what the with command for reading in external packages does.) Then the format of the spline command is similar to the interp command from the last lab.

                     spline(xlist,ylist,var);

The output from the spline command is the piecewise cubic definition of the natural cubic spline function S(x) – in Maple terms it is an expression, not a function, though(!) You can assign this to a symbolic name, plot, use it to define a function, or whatever you like.

Laboratory Questions

A) In class on Tuesday, we worked out an example where x_{0} = 1, x_{1} = 1, x_{2} = 3, x_{3} = 4 and f(x_{0}) = a_{0} = 5, f(x_{1}) = a_{1} = 4, f(x_{2}) = a_{2} = 0, f(x_{3}) = a_{3} = 1. We found that the cubic spline was given by

\{ \begin{matrix}S_{0}(x) = 5 − \frac{5}{8}x − \frac{3}{8}x^{3}, & 0 ≤ x ≤ 1 \\ S_{1}(x) = 4 − \frac{7}{4}(x − 1) − \frac{9}{8}(x − 1)^{2} + \frac{1}{2}(x − 1)^{3}, & 1 ≤ x ≤ 3 \\ S_{2}(x) = −\frac{1}{4}(x − 3) + \frac{15}{8}(x − 3)^{2} − \frac{5}{8}(x − 3)^{3}, & 3 ≤ x ≤ 4 \end{matrix}

Verify that S satisfies the continuity, differentiability, and second-order differentiability conditions for cubic spline functions. Plot S(x) and the points (0, 5), (1, 4), (3, 0), (4, 1) together (see the Maple online help page for the plot command for some ideas about plotting sets of points and use the option style = point).

B) Use Maple’s spline command to do problems 3d and 4d in section 3.4 of Burden and Faires. Also, plot the spline together with f(x) = x cos x − 2x^{2} + 3x − 1, and discuss the goodness of the “fit” between the spline function and f(x).

C) Use Maple’s spline command to do problem 27 in section 3.4 of Burden and Faires.

D) Maple does not have a built-in procedure or a library package for computing clamped cubic splines. However, you can get the source code for one by entering the command:

cp /home/fac/little/Classes/Numer-98/clampedspine.map clampedspline.map

(The source listing is also available on the course homepage.) Do this and read the code into your Maple session. The format for using the procedure is similar to the one for spline, except you also need to supply the values for the derivative of the spline function at the endpoints (df0, dfn are f′(x_{0}) and f′(x_{n}) respectively):

              clampedspline(xlist,ylist,df0,dfn,var);

Use this procedure to solve problem 28 in section 3.4 of Burden and Faires.

Assignment

Individual lab write-ups due in class Thursday, November 12.