Polynomial basis for Hermite interpolation
Pedro Gonnet, September 2010
(Chebfun example approx/HermiteBasis.m)
This script shows how one can use Chebfun to compute a polynomial basis for Hermite interpolation. By Hermite interpolation here we mean approximating a function f by polynomials that interpolate f and a specified number of its derivatives at the endpoints.
For simplicity, we will assume that everything is computed on the interval [-1,1].
d = domain(-1,1);
To evaluate our function on this domain at the endpoints, we can create a basic linop which, when applied to a chebfun, evaluates it at the nodes -1 and 1.
I = eye(d);
L = [];
L{1} = I([1 end],:);
By differentiating this operator, we can construct operators that compute the first k derivatives at the nodes -1 and 1.
k = 3; D = diff(d); for i=2:k+1 L{i} = L{i-1} * D; end;
If we interpolate k derivatives at -1 and 1, our basis will be of degree 2*k+1. To construct this basis, we evaluate the linops at the 2*k+2 Chebyshev nodes and assemble a (2*k+2)x(2*k+2) matrix
M = []; for i=1:k+1 M = [ M ; L{i}(2*k+2) ]; end;
and solve for the unit matrix
V = M \ eye(2*k+2);
The columns of the matrix V now contain the basis polynomials evaluated at the 2*k+2 Chebyshev points. We can recover the basis as Chebfuns by applying the Chebfun constructor to each column.
B = []; for i=1:2*k+2 B = [ B , chebfun( V(:,i) ) ]; end;
We can plot the individual basis functions as
for i=1:2*k+2 subplot(k+1,2,i); plot( B(:,i) , 'LineWidth' , 2 ); end;
This basis can now be used to construct a Hermite interpolant of any function from the function values and derivatives at the interval edges.
f = chebfun( @(x) exp(x).*sin(pi*(x+1)) ); df = diff(f); ddf = diff(df); dddf = diff(ddf); b = [ f(-1) ; f(1) ; df(-1) ; df(1) ; ddf(-1) ; ddf(1) ; dddf(-1) ; dddf(1) ]; h = B * b figure subplot(1,2,1); plot( [ f , h ] , 'LineWidth' , 2 ); title('approximation'); subplot(1,2,2); plot( [ f - h ] , 'LineWidth' , 2 ); title('error');
h =
chebfun column (1 smooth piece)
interval length endpoint values
[ -1, 1] 8 -3.4e-15 -6.7e-16
vertical scale = 1.8
Conversely, we can find the function values and derivatives that would match the original function best in the least-squares sense by computing
b = B \ f; h = B * b; figure subplot(1,2,1); plot( [ f , h ] , 'LineWidth' , 2 ); title('approximation'); subplot(1,2,2); plot( [ f - h ] , 'LineWidth' , 2 ); title('error');
