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


disp(sprintf('**** Example 7.1 Monomial Basis *'));
echo on
%%%% Use polynomial interpolation with monomial basis p(t)=x_1+x_2t+x_3t^2.
		% Find coefficients x_i.
t=[-2 0 1]';	% independent variable
y=[-27 -1 0]';	% measured data values
A=[1 t(1) t(1)^2;1 t(2) t(2)^2;1 t(3) t(3)^2; ]	% Vandermonde matrix

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

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


disp(sprintf('**** Example 7.2 Lagrange Interpolation *'));
echo on
%%%% Use polynomial interpolation with Lagrange basis p(t)=y_1*l_1(t)+y_2*l_2(t)+y_3*l_3(t).
		% where l_j(t)=PI(t-t_k)/PI(t_j-t_k) for k,j=1..n and k not j.
t=[-2 0 1]';	% independant variable
y=[-27 -1 0]';	% measured data values

%%%% There is no need for the solution of the system because A is I so x_i=y=i
%%%% The interpolating polynomial is:
tt=sym('tt');	% symbolic variable tt
p3=simplify(y(1)*((tt-t(2))*(tt-t(3)))/((t(1)-t(2))*(t(1)-t(3))) + ...
   y(2)*((tt-t(1))*(tt-t(3)))/((t(2)-t(1))*(t(2)-t(3)))+ ...
   y(3)*((tt-t(1))*(tt-t(2)))/((t(3)-t(1))*(t(3)-t(2))))
% The evaluation is more complex.


disp(sprintf('**** Example 7.3 Newton Interpolation *'));
echo on
%%%% Use polynomial interpolation with Newton basis p(t)=x_1*n_1(t)+x_2*n_2(t)+x_3*n_3(t).
		% where n_j(t)=PI(t-t_k) for j=1..n and k=1..j-1. Find coefficients x_i.
t=[-2 0 1]';	% independant variable
y=[-27 -1 0]';	% measured data values

A=[1 0 0;1 t(2)-t(1) 0;1 t(3)-t(1) t(3)-t(1)*t(3)-t(2); ]	

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

tt=sym('tt');	% symbolic variable tt
p3=simplify (x(1)+x(2)*(tt-t(1))+x(3)*(tt-t(1))*(tt-t(2))) %Interpolating polinomial 
% It can be evaluated for tt efficiently by Horner's algorithm.



disp(sprintf('**** Example 7.x Plot Monomial, Lagrange and Newton basis functions *'));
echo on
%%%% PLot basis functions for 5 equidistant interpolation points.
t=[0 .25 .5 .75 1];

echo off		% Plot monomial basis functions (independent from data points)
f=zeros(101,5);	% Space for function values
xx=0:0.01:1;		% Plot 5 monomial basis functions
for k=1:5
   for i=1:101
      f(i,k)=xx(i)^(k-1);
   end
end
figure(2);
plot(xx,f)
title('Monomial basis functions t=[0 .25 .5 .75 1].')
echo on

echo off		% Plot Lagrange basis functions (depend on t)
f=zeros(101,5);	% Space for function values
xx=0:0.01:1;		% Plot 5 Lagrange basis functions for t
for i=1:101
   f(i,1)=((xx(i)-t(2))*(xx(i)-t(3))*(xx(i)-t(4))*(xx(i)-t(5)))/...
     		((t(1)-t(2))*(t(1)-t(3))*(t(1)-t(4))*(t(1)-t(5)));
   f(i,2)=((xx(i)-t(1))*(xx(i)-t(3))*(xx(i)-t(4))*(xx(i)-t(5)))/...
     		((t(2)-t(1))*(t(2)-t(3))*(t(2)-t(4))*(t(2)-t(5)));
   f(i,3)=((xx(i)-t(1))*(xx(i)-t(2))*(xx(i)-t(4))*(xx(i)-t(5)))/...
     		((t(3)-t(1))*(t(3)-t(2))*(t(3)-t(4))*(t(3)-t(5)));
   f(i,4)=((xx(i)-t(1))*(xx(i)-t(2))*(xx(i)-t(3))*(xx(i)-t(5)))/...
     		((t(4)-t(1))*(t(4)-t(2))*(t(4)-t(3))*(t(4)-t(5)));
   f(i,5)=((xx(i)-t(1))*(xx(i)-t(2))*(xx(i)-t(3))*(xx(i)-t(4)))/...
     		((t(5)-t(1))*(t(5)-t(2))*(t(5)-t(3))*(t(5)-t(4)));
end
figure(3);
plot(xx,f)
title('Lagrange basis functions for t=[0 .25 .5 .75 1].')
echo on

echo off		% Plot Newton basis functions (depend on t)
f=zeros(101,5);	% Space for function values
xx=0:0.01:1;		% Plot 5 Newton basis functions for t
for i=1:101
   f(i,1)=1;
   f(i,2)=xx(i)-t(1);
   f(i,3)=(xx(i)-t(1))*(xx(i)-t(2));
   f(i,4)=(xx(i)-t(1))*(xx(i)-t(2))*(xx(i)-t(3));
   f(i,5)=(xx(i)-t(1))*(xx(i)-t(2))*(xx(i)-t(3))*(xx(i)-t(4));
end
figure(4);
plot(xx,f)
title('Newton basis functions for t=[0 .25 .5 .75 1].')
echo on


disp(sprintf('**** Example 7.4 Incremental Newton Interpolation *'));
echo on
%%%% Start with the first data point and build the polynomials of higher order 
		% by adding further data points.
t=[-2 0 1]';	% independent variable
y=[-27 -1 0]';	% measured data values

tt=sym('tt');	% symbolic variable tt
p1=y(1)		% first data point interpolated by the constant polynomial

x2=(y(2)-y(1))/(t(2)-t(1));
p2=simplify(y(1)+x2*(tt-t(1)))	% second data point gives the new polynomial of higher degree

x3=(y(3)-(y(1)+x2*(t(3)-t(1))))/((t(3)-t(1))*(t(3)-t(2)));
p3=simplify(p2+x3*((tt-t(1))*(tt-t(2))))	% third data point gives the new polinomial of higher degree


disp(sprintf('**** Example 7.5 Divided Differences for Newton Interpolation *'));
echo on
%%%% Start with values in data points, build recursively divided differences
t=[-2 0 1]';	% independent variable
y=[-27 -1 0]';	% measured data values

tt=sym('tt');	% symbolic variable tt
f1=y(1);		% first differences are data point values
f2=y(2);
f3=y(3);

f12=(f2-f1)/(t(2)-t(1));
f23=(f3-f2)/(t(3)-t(2));

f123=(f23-f12)/(t(3)-t(1));

p3=simplify(f1 + f12*(tt-t(1)) + f123*((tt-t(2))*(tt-t(1)))) %Again the interpolation polinomial


disp(sprintf('**** Example 7.x Generate and Plot Legendre and Chebyshev polynomials *'));
echo on
%%%% Generate P_k and T_k orthogonal polynomials by three-term recurrence and
		% plot polynomials for 5 equidistant interpolation points.
t=[0 .25 .5 .75 1];

%%%% Generate Legendre polynomials P_k
tt=sym('tt');		% symbolic variable tt
P0=1
P1=simplify(tt)
k=1;
P2=simplify(((2*k+1)*tt*P1-k*P0)/(k+1))
k=2;
P3=simplify(((2*k+1)*tt*P2-k*P1)/(k+1))
k=3;
P4=simplify(((2*k+1)*tt*P3-k*P2)/(k+1))

echo off		% Plot Legendre polynomials P_k
P=zeros(101,5);	% Space for function values
xx=-1:0.01:1;	
for i=1:201
   tt=xx(i);
	P(i,1)=1;
   P(i,2)=eval(P1);
   P(i,3)=eval(P2);
   P(i,4)=eval(P3);
   P(i,5)=eval(P4);
end
figure(5)
plot(xx,P)
title('Legendre polynomials P_k for t=[0 .25 .5 .75 1].')
echo on

%%%% Generate Chebyshev polynomials T_k
tt=sym('tt');		% symbolic variable tt
T0=1
T1=simplify(tt)
k=1;
T2=simplify(2*tt*T1-T0)
k=2;
T3=simplify(2*tt*T2-T1)
k=3;
T4=simplify(2*tt*T3-T2)

echo off		% Plot Legendre polynomials T_k
T=zeros(101,5);	% Space for function values
xx=-1:0.01:1;	
for i=1:201
   tt=xx(i);

	T(i,1)=1;
   T(i,2)=eval(T1);
   T(i,3)=eval(T2);
   T(i,4)=eval(T3);
   T(i,5)=eval(T4);
end
figure(6)
plot(xx,T)
title('Chebyshev  polynomials T_k for t=[0 .25 .5 .75 1].')
echo on


disp(sprintf('**** Example 7.x Interpolation by Legendre Polynomials *'));
echo on
%%%% Use Legendre orthogonal basis polynomial for interpolation three data 
		% points p(t)=x_1*P_0(t)+x_2*P_1(t)+x_3*P_2(t), where P_j(t) is 
		% Legendre polinomial (using (k+1)P_(k+1)=(2k+1)*t*P_k-kP_(k-1)). 
		% Find coefficients x_i.
t=[-2 0 1]';		% independent variable
y=[-27 -1 0]';		% measured data values

A=[1 t(1) (3*t(1)^2-1)/2; 1 t(2) (3*t(2)^2-1)/2; 1 t(3) (3*t(3)^2-1)/2; ]	

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

tt=sym('tt');		% symbolic variable tt
p3=simplify (x(1) + x(2)*tt + x(3)*(3*tt^2-1)/2) %Interpolating polynomial same as before


disp(sprintf('**** Example 7.6 Cubic Spline Interpolation *'));
echo on
%%%% Use a piecewise interpolation with polynomials of degree 3.
t=[-2 0 1]';		% independent variable
y=[-27 -1 0]';		% measured data values

A=[1 t(1) t(1)^2 t(1)^3 0 0 0 0; 1 t(2) t(2)^2 t(2)^3 0 0 0 0;...
   0 0 0 0 1 t(2) t(2)^2 t(2)^3; 0 0 0 0 1 t(3) t(3)^2 t(3)^3;...
   0 1 2*t(2) 3*t(2)^2 0 -1 -2*t(2) -3*t(2);...
   0 0 2 6*t(2) 0 0 -2 -6*t(3); 0 0 2 6*t(1) 0 0 0 0; 0 0 0 0 0 0 2 6*t(3)]
b=[y(1) y(2) y(2) y(3) 0 0 0 0]'
x=A\b
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)

echo off		% Plot Cubic Spline interpolants
C=zeros(301,1);	% Space for function values
xx=-2:0.01:1;	
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
figure(1)
plot(xx,yy, t,y,'ro',xx,C) % Plot monomial again - some other solution?
title('Data points (-2 -27),(0 -1),(1 0)and their Interpolants.')
text(-1.8,-3,'Monomial Interpolation for n=3')
text(-1,-18,'Cubic Spline')
echo on