%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function for calculation/visualisation of an
% Hermite interpolation example for the 4th exercise.
%
% The Hermite-polynomial for n+1 nodes is calculated
% for an equidistant mesh on [x0;xn]. Besides the values
% yi, the tangents y'i are calculated (for a function
% defined in the subroutine eval_f).
%
% author: Tobias Neckel
% date: 16.11.05
%
% INPUT: numberOfIntervals number of subintervals n
%
% OUTPUT: -
%
% ATTENTION: matlab starts arrays with index 1 instead of 0,
% so we have to shift indices when using xi, yi!
%
function visHermiteSin(numberOfIntervals)
% set global interval
x0=0;
xn=4*pi;
% get interpolation nodes
xi = (x0 : (xn-x0)/numberOfIntervals : xn)';
[yi, yprime] = eval_f(xi);
% get real function
x = linspace(x0,xn, 50*numberOfIntervals)';
[f,df]=eval_f(x);
% calculate Hermite polynomial
ipp1 = hermite(xi,yi,yprime,x);
% calculate maximum pointwise difference
error_max = max(abs(ipp1-f));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% visualisation: plot f and ipp1
figure_handle=figure;
%%%%%%%%%%%%%%%%%%
% plot for matlab
plot_handle=plot(x,f,'b-.',x,ipp1,'r',xi,yi,'ko');
set(plot_handle,'LineWidth',2);
set(gca,'LineWidth',2,'FontWeight','bold');
titel = sprintf('Hermite-interpolation with %2i nodes',numberOfIntervals+1);
title([titel ' auf [0;4\pi]']);
legend('f(x) = sin(x)','Hermite-polynomial','interpolation nodes');
axis auto;
text(9.5, 0.5, ['max. error = ' num2str(error_max)],'FontWeight','bold');
%%%%%%%%%%%%%%%%%%
% plot for octave
% plot(x,f,'b;f(x) = sin(x);',x,ipp1,'r;Hermite-polynomial;',xi,yi,'co;interpolation nodes;');
% titel = sprintf('Hermite-interpolation with %2i nodes',numberOfIntervals+1);
% title(titel);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subroutine for calculation of the values and
% derivatives of a given function f (f(x) == sin(x))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y,dy]=eval_f(xi)
y = sin(xi);
dy = cos(xi);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subroutine for evaluation of the Hermite interpolant
% at a point x in [x0;xn] for values y0,.., yn, y'0,..., y'n.
%
% The four basis polynomials H0,...,H3 are used to
% evaluate the polynomial on each subinterval [xi;xi+1].
% The transformation factor h for yprime is already
% considered inside the loop and therefore doesn't
% appear in the final evaluation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hermitePoly = hermite(xi,yi,yprime,x);
t = zeros(length(x),1);
y0 = t;
y1 = t;
yprime0 = t;
yprime1 = t;
for i = 1:length(x);
if i == 1
left = xi(1);
right = xi(2);
else
left = xi( max( find(xi < x(i)) ) );
right= xi( min( find(xi>= x(i)) ) );
end
t(i) = (x(i) - left)./(right - left);
y0(i) = yi' * (xi==left );
y1(i) = yi' * (xi==right);
% factor h considered here!
yprime0(i) = (right - left)*yprime' * (xi == left );
yprime1(i) = (right - left)*yprime' * (xi == right);
end
[Hzero, Hone, Htwo, Hthree] = HermiteBasis(t);
hermitePoly = y0.*Hzero + y1.*Hone + yprime0.*Htwo + yprime1.*Hthree;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subroutine for evaluation of the Hermite basis
% polynomials H0,...,H3 (vector values possible)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Hzero, Hone, Htwo, Hthree]=HermiteBasis(t)
Hzero = (1 + 2.*t) .* (1 - t).^2;
Hone = (3 - 2.*t) .* (t.^2);
Htwo = t .* (1 - t).^2;
Hthree = (t - 1) .* (t.^2);
return