problem 9.3.5.1

1. Write a program to compute the piecewise Hermite cubic interpolant to f(x) = 1 / (x^2 + 25)

at equidistant nodes in the interval [-1, 1]. Plot f and the interpolant for 2 and 3 nodes. By

sampling the error at 100 points between each node, determine a (computer generated) estimate for

the interpolation error, and plot the log of the error versus the log of the mesh width.

Solution:

The idea is to input n+1 nodes x0, x1, ..., xn between [-1, 1]. Evaluate function values and derivatives at the nodes f0, f1, ... ,fn and f0', f1', ...fn'. There are n intervals. In each interval [xi, xi+1], find a polynomial of order 3

p(x) = a0 + a1*x + a2*x^2 +a3^3 s.t.

p(xi) = f(xi), p(xi+1) = f(xi+1), p'(xi) = f'(xi), p'(xi+1) = f'(xi+1).

There are four unknown coefficients and four conditions. The solution is (page 312 of the text book)

p(x) = f(xi) + f'(xi)(x - xi) + f[xi, xi, xi+1] (x - xi)^2 + f[xi, xi, xi+1, xi+1] (x - xi)^2*(x - xi+1).

The results of f(x), p(x) and error=f(x)-p(x) as functions of x at N = 2, 3, 4 are show below

rungehermite.c can be compiled by "icc -o run rungehermite.c", which use the node number "N" as the input parameter. For example, "./run 2" generates Hermite interpolant using N = 2 nodes.

plotpx plot f(x), p(x) as a function of x at certain N.

ploterror plot error = f(x) - p(x) as a function of x at certain N.

plotmaxerror plot max(f(x) - p(x)) vs mesh width on double log scale using data from the following table.

Maximum error for increasing number of notes

We draw column 2 and 3 on double log scale and log(max error) decreases linearly with decreasing log(mesh width).