%
% MATLAB for DUMMIES 11-12
% Exercice 2
%
% Solution detaillee
%  Vincent Legat, Damien Bouvy, Denis-Gabriel Caprace, Gilles De Neyer, 
%  Kevin Degraux, Mohammad Iranmanesh, Christophe Lorant, Remy Moulin, 
%  Matthew Philippe, Karim Slaoui, Olivier Thiry, Laurent Van Eesbeeck,
%  Thomas Vanderstraeten, Antoine Weynants
%
% -------------------------------------------------------------------------
function matlab2()

n = 9; 
X = linspace(0,1,n+1);
U = exp(-2*X).*sin(10*pi*X);
dU = 10*pi*exp(-2*X).*cos(10*pi*X) -2 *U;

x = linspace(0,1,100);
u = exp(-2*x).*sin(10*pi*x);
uh = hermite(n,X,U,dU,x);

close all;
plot(x,u,'b-', x, uh, 'r-', X,U,'b.','Markersize',25);
fprintf('==== uh(22) %14.7e \n',uh(22));
 
end

function [uh] = hermite(n,X,U,dU,x)
%HERMITE Piecewise Cubic Hermite Interpolation
%   UH = HERMITE(n,X,U,dU,x) computes the value
%   of the piecewise cubic interpolation for abscisses defined by x.
%   The interpolation is defined by the values U and the derivatives dU at
%   points X.   The number of intervals is given by n and
%   The vectors X, U and dU must contain at least n+1 elements.
%   It is assumed that the values of X are different and sorted :-)
%
%   As testing input arguments is boring, we do not care about.
%   It is not perfect, but it is the life !
%

%
%   -1- Finding all coefficients of the cubic local polynomials
%
%   On the interval [X(i-1),X(i)], the Hermite interpolation is : 
%
%       uh{i}(x) = a*(x-X(i-1))^3 + b*(x-X(i-1))^2 + c*(x-X(i-1)) + d.
% 
%   The coefficients a,b,c,d are solution of :
%
%       [ 0     0     0   1] [a]   [ U(i-1)]
%       [ 0     0     1   0] [b] = [dU(i-1)]
%       [dx^3  dx^2  dx   1] [c]   [ U(i)  ]
%       [3dx^2 2dx    1   0] [d]   [dU(i)  ]
% 
%   with dx = X(i)-X(i-1).
%   One can analytically obtain the coefficients :
%       
%       a = (dU(i) + dU(i-1) - 2dudx)/dx^2
%       b = (3dudx - 2dU(i-1) - dU(i))/dx   
%       c = dU(i-1)
%       d = U(i-1)
%
%   with dudx = (U(i)-U(i-1))/dx.
%

dx = diff(X);
du = diff(U);
dudx = du ./ dx;

a = (dU(1:n) + dU(2:n+1) - 2* dudx) ./ (dx.*dx);
b = (3*dudx - 2*dU(1:n) - dU(2:n+1)) ./ dx;
c = dU(1:n);
d = U(1:n);

%
%   -2- If you do not believe me :-)
%
%   Of course, it is also possible to solve each small system
%   by the Gaussian elimination...
%   A nice double check of your analytical algebra, but
%   it is not a really very efficient way to work !
%
%   Just uncomment the lines below to check that the analytical
%   and the numerical solutions are the same :-)
%
% for i=1:n
%     A = [ 0         0        0       1 ; 
%           0         0        1       0 ; 
%           dx(i)^3   dx(i)^2  dx(i)   1 ; 
%           3*dx(i)^2 2*dx(i)  1       0];
%     coeff = A \ [U(i) dU(i) U(i+1)  dU(i+1)]';
%     a(i) = coeff(1);
%     b(i) = coeff(2);
%     c(i) = coeff(3);
%     d(i) = coeff(4);
% end
%

%
%   -3- Selecting for each value of x the good polynomial
%

indice = zeros(size(x));
for i=1:n
    indice = indice + (x >= X(i));
end

%
%   -4- Computing the polynomial with the Horner's technique
%       in order to minimize floating point errors :-) 
%

uh = a(indice).*(x - X(indice)) + b(indice);
uh = uh.*(x - X(indice)) + c(indice);
uh = uh.*(x - X(indice)) + d(indice);

end