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 -1 0 1 2]';	% equidistant independent variable
y=[-27 -10 0 -6 10]';	% measured data values
n=size(t);		% number of points
xx=t(1):0.01:t(n);	% points for ploting	

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.
yy=zeros(401,1);
for i=1:401
   tt=xx(i);
   yy(i)=eval(p5);
end
echo on
figure(1);
plot(xx,yy,'b', t,y,'ro')
title('Data points (-2 -27),(-1 5),(0 0),(1 -6),(2 10) and their Interpolants.')
text(0.3,-10,'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(400,1);	% Space for function values
for i=1:101			% First cubic polynomial 
   tt=xx(i);
   C(i,1)=eval(C1);
end
for i=101:201		% Second cubic polynomial
   tt=xx(i);
   C(i,1)=eval(C2);
end
for i=201:301		% Third cubic polynomial
   tt=xx(i);
   C(i,1)=eval(C3);
end
for i=301:401		% 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.8,0,'Cubic monomial Spline','Color','g')
echo on


disp(sprintf('**** Test - linear B-splines *'));
echo on
%%%% Use polynomial interpolation with cubic B-splines
%		p(t)=x_1*B0(t)+x_2*B1(t).
% Find coefficients x_i.
%Matrix A is banded, main diagonal is 1 and subdiagonals 1/4
A=diag(ones(1,n))

x=A\y		% Solve dense system to get coefficients x, which are equal to y in this case.

%Tabulate cubic B-spline (only ones)
echo off
B=zeros(201,1);	% Space for linear B-spline 
h=t(2)-t(1);
hh=-h:0.01:h; % interval for evaluation +-h
for i=1:101			% left interval 
   tt=hh(i);
   B(i,1)= (1/h)*(tt-t(2));
end
for i=101:201		% right interval (mirror of left part if h constant)
   tt=hh(i);
   B(i,1)= (1/h)*(t(4)-tt);
end
% Evaluate linear B-spline by segments
S=zeros(401,1);	% Space for interpolated function values by linear B-spline
for i=1:101			% calculate by intervals 
   S(i,1)= x(1)*B(i+100,1) + x(2)*B(i,1);
   S(i+100,1)= x(2)*B(i+100,1) + x(3)*B(i,1);
   S(i+200,1)= x(3)*B(i+100,1) + x(4)*B(i,1);
   S(i+300,1)= x(4)*B(i+100,1) + x(5)*B(i,1);
end
figure(1)
hold on;
plot(xx,S,'r') % Plot monomial again - some other solution?
text(-0.5,-7,'Linear B-Spline','Color','r')
echo on


disp(sprintf('**** Test - cubic B-splines *'));
echo on
%%%% Use polynomial interpolation with cubic B-splines
%		p(t)=x_1*B0(t)+x_2*B1(t)+x_3*B2(t)+x_4*B3(t).
% Find coefficients x_i.
%Matrix A is banded, main diagonal is 1 and subdiagonals 1/4
A=diag(ones(1,5))+diag(1/4*ones(1,4),1)+diag(1/4*ones(1,4),-1)

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

%Tabulate cubic B-spline (only ones)
echo off
B=zeros(401,1);	% Space for cubic B-spline 
h=t(2)-t(1);
hh=-2*h:0.01:2*h; % interval for evaluation +-2*h
for i=1:101			% second left interval 
   tt=hh(i);
   B(i,1)= 1/(4*h^3)*(tt-t(1))^3;
end
for i=101:201		% left central interval
   tt=hh(i);
   B(i,1)=1/4*(1 + 3/h*(tt-t(2)) + 3/h^2*(tt-t(2))^2 - 3/h^3*(tt-t(2))^3);
end
for i=201:301		% right central interval
   tt=hh(i);
   B(i,1)=1/4*(1 + 3/h*(t(4)-tt) + 3/h^2*(t(4)-tt)^2 - 3/h^3*(t(4)-tt)^3);
end
for i=301:401		% second right interval (mirror of left part if h constant)
   tt=hh(i);
   B(i,1)= 1/(4*h^3)*(t(5)-tt)^3;
end

% Evaluate cubic B-spline by segments
S=zeros(401,1);	% Space for interpolated function values by cubic B-spline
for i=1:101			% calculate by intervals 
   S(i,1)= x(1)*B(i+200,1) + x(2)*B(i+100,1) + x(3)*B(i,1);
   S(i+100,1)= x(1)*B(i+300,1) + x(2)*B(i+200,1) + x(3)*B(i+100,1) + x(4)*B(i,1);
   S(i+200,1)= x(2)*B(i+300,1) + x(3)*B(i+200,1) + x(4)*B(i+100,1) + x(5)*B(i,1);
   S(i+300,1)= x(3)*B(i+300,1) + x(4)*B(i+200,1) + x(5)*B(i+100,1) ;
end

figure(1)
hold on;
plot(xx,S,'k') % Plot monomial again - some other solution?
text(-1.8,-5,'Cubic B-Spline','Color','k')

