Code covered by the BSD License  

Highlights from
Chebfun V4

image thumbnail

Chebfun V4

by

 

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

Contact us