echo on
%%%%% ch07_exa.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;
format short;
disp(sprintf('**** Examples from Chapter 07 with 5 data points ****'));

t=[-2 0 1 2 2.5 ]';	% independent variable
y=[-27 -1 0 -4 10]';	% measured data values

disp(sprintf('**** Example 7.1 Monomial Basis *'));
echo on
%%%% Use polynomial interpolation with monomial basis 
%		p(t)=x_1+x_2*t+x_3*t^2+x_4*t^3+x_5*t^4.
		% Find coefficients x_i.
A=[1 t(1) t(1)^2 t(1)^3 t(1)^4;...
	1 t(2) t(2)^2 t(2)^3 t(2)^4;...
	1 t(3) t(3)^2 t(3)^3 t(3)^4;...
	1 t(4) t(4)^2 t(4)^3 t(4)^4;...
	1 t(5) t(5)^2 t(5)^3 t(5)^4;...
 ]	% Vandermonde matrix

x=A\y		% Solve dense system to get coefficients.

tt=sym('tt');	% symbolic variable tt
p5=simplify (x(1)+x(2)*tt+x(3)*tt^2+x(4)*tt^3+x(5)*tt^4)	% Interpolating polynomial of second order
% It can be evaluated for t efficiently by Horner's algorithm.
echo off		% Plot data and interpolant.
xx=t(1):0.01:t(5);
for i=1:451
   tt=xx(i);
   yy(i)=eval(p5);
end
echo on
figure(1);
plot(xx,yy,'b', t,y,'ro')
title('Data points (-2 -27),(0 -1),(1 0),(2 -4),(2.5 10) and their Interpolants.')
text(-0.5,-15,'Monomial Interpolation for n=5','Color','b')


disp(sprintf('**** Example 7.6 Cubic Spline Interpolation *'));
echo on
%%%% Use a piecewise interpolation with polynomials of degree 3, one for 
% each pair of points.
%(n-1) polynomials with 4 paramethers have to be found so A has dimension
%of 4*(n-1)=16
	%knote values equal to the measured values
A=[1 t(1) t(1)^2 t(1)^3    0 0 0 0   0 0 0 0   0 0 0 0;...
   1 t(2) t(2)^2 t(2)^3    0 0 0 0   0 0 0 0   0 0 0 0;...
	0 0 0 0   1 t(2) t(2)^2 t(2)^3 	 0 0 0 0   0 0 0 0;...
 	0 0 0 0   1 t(3) t(3)^2 t(3)^3 	 0 0 0 0   0 0 0 0;...
 	0 0 0 0   0 0 0 0   1 t(3) t(3)^2 t(3)^3 	  0 0 0 0;...
 	0 0 0 0   0 0 0 0   1 t(4) t(4)^2 t(4)^3 	  0 0 0 0;...
 	0 0 0 0   0 0 0 0   0 0 0 0    1 t(4) t(4)^2 t(4)^3;...
	0 0 0 0   0 0 0 0   0 0 0 0    1 t(5) t(5)^2 t(5)^3;...
   %first derivatives in inner knots are equal
   0 1 2*t(2) 3*t(2)^2   0 -1 -2*t(2) -3*t(2)^2   0 0 0 0  0 0 0 0;...
   0 0 0 0  0 1 2*t(3) 3*t(3)^2   0 -1 -2*t(3) -3*t(3)^2   0 0 0 0;...
   0 0 0 0  0 0 0 0   0 1 2*t(4) 3*t(4)^2   0 -1 -2*t(4) -3*t(4)^2;...
   %second derivatives in inner knots are equal
	0 0 2 6*t(2)  0 0 -2 -6*t(2) 0 0 0 0  0 0 0 0;...
 	0 0 0 0  0 0 2 6*t(3)  0 0 -2 -6*t(3) 0 0 0 0;...
    0 0 0 0  0 0 0 0  0 0 2 6*t(4)  0 0 -2 -6*t(4);...
   %second derivatives in boundary knots = 0 (normal condition)
	0 0 2 6*t(1) 0 0 0 0  0 0 0 0  0 0 0 0;...
   0 0 0 0  0 0 0 0  0 0 0 0   0 0 2 6*t(5);...
]
b=[y(1) y(2) y(2) y(3)  y(3) y(4) y(4) y(5)  0 0 0   0 0 0   0 0]'
x=A\b
%x=[4 coeficients for 4 polynomials] 
tt=sym('tt');		% symbolic variable tt
C1=simplify(x(1)+x(2)*tt+x(3)*tt^2+x(4)*tt^3)
C2=simplify(x(5)+x(6)*tt+x(7)*tt^2+x(8)*tt^3)
C3=simplify(x(9)+x(10)*tt+x(11)*tt^2+x(12)*tt^3)
C4=simplify(x(13)+x(14)*tt+x(15)*tt^2+x(16)*tt^3)

echo off		% Plot Cubic Spline interpolants
C=zeros(451,1);	% Space for function values
xx=t(1):0.01:t(5);	
for i=1:201			% First cubic polynomial 
   tt=xx(i);
   C(i,1)=eval(C1);
end
for i=201:301		% Second cubic polynomial
   tt=xx(i);
   C(i,1)=eval(C2);
end
for i=301:401		% Third cubic polynomial
   tt=xx(i);
   C(i,1)=eval(C3);
end
for i=401:451		% Fourth cubic polynomial
   tt=xx(i);
   C(i,1)=eval(C4);
end
figure(1)
hold on;
plot(xx,C,'g') % Plot monomial again - some other solution?
text(-1.5,-7,'Cubic Spline','Color','g')
echo on