Mathematics 371 – Numerical Analysis Lab Day 3 – Using Cubic Spline Interpolation October 19, 2001

Goals and Background

The goals of today’s lab are to use 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 Wednesday 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 − 1. 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

As we saw in the handout from Monday’s class, 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:

                        with(CurveFitting);

Then the format of the spline command is similar to the PolynomialInterpolation command from the last problem set.

                     Spline(xlist,ylist,var);

where xlist is the list of x-coordinates of the interpolation points (that is, the x_{j}), and ylist is the list of y-coordinates f(x_{j}). These are Maple lists so, they should be lists of floating point numbers, separated by commas in square brackets: [ , ]. 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 make it into a function, assign it to a symbolic name, plot it, or whatever you like.

Laboratory Questions

A) In class on Wednesday, we worked out an example where x_{0} = 0, 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}

Use the Spline command to check this work, 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 (in particular, the option style = point). You will probably want to make the combined plot of S(x) and the points by using the display command from the plots package: Generate the two plots separately and assign them to symbolic names (this is one case where you probably don’t want to see the output, so use a colon at the end of the plot commands), then use display combine them.

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. Note that this is the same collection of data points as you used in problem 25 from section 3.1 on the last problem set, but there you used interpolating polynomials. Compare and contrast your results using interpolating polynomials and interpolating splines.

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 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 Friday, October 26.