% Example 4.7.1 
% Compare the errors in approximations of
%   f(x)=exp(x) on [a,b]
% made by 
% least-squares polynomials % provided by "polyfit"
% Chebyshev polynomials     % provided by "chebyshev_interp_1d"
% Also shows the cubic spline interpolant.

clear all
a = -5; b = 5;

%yi = chebyshev_interp_1d ( nd, xd, yd, ni, xi )
%
%    Input, integer ND, the number of data points, number of nodes.
%       ND must be at least 1.
%    Input, real XD(ND,1), the data locations / nodes.
%    Input, real YD(ND,1), the data values.
%    Input, integer NI, the number of interpolation points.
%    Input, real XI(NI,1), the interpolation points, which
%      must be each be in the interval [ min(XD), max(XD)].
%    Output, real YI(NI,1), the interpolated values. 
%
%yi = chebyshev_value_1d ( nd, c, xmin, xmax, ni, xi )
%
%    Input, integer ND, the number of data points.
%    ND must be at least 1.
%    Input, real C(ND,1), the Chebyshev coefficients.
%    Output, real XMIN, XMAX, the interpolation interval.
%    Input, integer NI, the number of interpolation points.
%    Input, real XI(NI,1), the interpolation points, which
%        must be each be in the interval [XMIN,XMAX].
%    Output, real YI(NI,1), the interpolated values.
%
%[ c, xmin, xmax ] = chebyshev_coef_1d ( nd, xd, yd )
%    Input, integer ND, the number of data points.
%    ND must be at least 1.
%
%    Input, real XD(ND,1), the data locations.
%    Input, real YD(ND,1), the data values.
%    Output, real C(ND,1), the Chebyshev coefficients.
%    Output, real XMIN, XMAX, the interpolation interval.
%

t = linspace(a,b,60)';
exact = cos(t);%exp(t);

nd = 25; % Number of nodes
x_nodes = linspace(a,b,nd)';% interpolating nodes
y_values = cos(x_nodes);%exp(x_nodes);% interpolating values

p=nd-1; %degree of orthogonal "polyfit" polynomial
fprintf('\nThe degree of orthogonal polynomials is %i \n',p)
fprintf('The number of interpolating points is %i \n\n',nd)
%
% Construct the Matlab's least-squares polynomial of degree p
LeastSquares = polyfit(x_nodes,y_values,p);
% Evaluate the interpolating polynomial
yLeastSquares = polyval(LeastSquares,t);

% Evaluate the Chebyshev interpolating polynomial 
 yChebyshev = chebyshev_interp_1d ( nd, x_nodes, y_values, length(t), t );
% Evaluate the Legendre interpolating polynomial 
% Li = p_polynomial_value ( m, n, x )

% Construct the Matlab's not-a-knot cubic splines interpolant
ycubicspline=spline(x_nodes,y_values,t);
% Evaluate the interpolating polynomial
%ycubicspline = polyval(ycubic,t);

format long
LinftyError_Chebyshev=max(exact-yChebyshev);
LinftyError_polyfit=max(exact-yLeastSquares);
LinftyError_cubicspline=max(exact-ycubicspline);
 
table = [ LinftyError_Chebyshev ; LinftyError_polyfit ; LinftyError_cubicspline ];
disp('|Error(Chebyshev)|_L^inf      |Error(polyfit)|_L^inf      |Error(cubic spline)|_L^inf')
format long
fprintf('        %.12f                    %.12f                      %.12f  \n', table)

figure(1)
hold all
plot(t,exact,'-y',x_nodes,y_values,'o')
plot(t,yLeastSquares,':b',t,ycubicspline,'-*g')
plot(t,yChebyshev,'.-.r')
legend('exact','data','polyfit Least-Squares interpolant','cubic spline','Chebyshev')