Code covered by the BSD License

# Chebfun V4

### Chebfun Team (view profile)

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

Polynomial basis for Hermite interpolation

# 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');