% Approximation of functions of the interval [1,1] using 
% the Gauss-Lobatto-Chebyshev approach

 
function spectr_06
   
echo off, clc
clear all

% Creating uniform grid (for comparison);
M = 100;
dx = 2.0 / (M - 1);
X0 = (-1.0:dx:1.0);

N = 17;      % N=65, 33, 17, ...
% Generating the Gauss-Lobatto grid points
X1(1:N) = cos( (0:N-1) * pi / (N - 1) );


u1 = Chebyshev_Coefficients( N, 1 );
I0 = Chebyshev_Interpol ( X0, N, u1 );


E0 = f( X0 );
E1 = f( X1 );

dE = ( I0 - E0 ) ./ abs( E0 );

figure(1); plot(X0,I0,':',X0,E0,'-',X1,E1,'o');
title('Chebyshev Approximation'); legend('Approx','Exact','GLC grid points');

figure(2); plot(X0,dE,'-');
title('Chebyshev Approximation - Error');



%-------------------------------------------------------------
% Determines the Chebyshev expansion coefficients via FFT
% N - number of grid-points (power of 2)

function [Uch] = Chebyshev_Coefficients( N, flag )  

% Checking of the number of grid points is consistent with FFT;
  if ( ~isposint(log2( 2 * (N-1) )) )
    disp(' Wrong number of grid-points!')
    return;
  end;

% Generating the Gauss-Lobatto grid points
  X(1:N) = cos( (0:N-1) * pi / (N - 1) );

% Evaluating the function on the grid-points
  U1(1:N) = f( X' );

% Calculating the cosine transform
  U2 = real(fft([U1'; flipud(U1(2:N-1)')]));

% Extracting the Chebyshev coefficients
  Uch = U2(1:N).*[0.5; ones(N-2,1); 0.5] / (N-1); 


%-------------------------------------------------------------
% Evaluating the Chebyshev interpolant
% x   - independet variable
% N   - order of the interpolant
% Uch - Chebyshev expansion coefficients

function [Ux] = Chebyshev_Interpol ( x, N, Uch )

Ux = 0.0;
for k=1:N
  Ux = Ux + Uch(k) * Tk( x, k-1 );
end;


%-------------------------------------------------------------
% Chebyshev polynomial of order k
% x - independent variable
% k - order of the polynomial

function [u] = Tk( x, k )
  if ( min(x) < -1 | max(x) > 1 )
    disp('Argument out of range! Exiting ...')
    return;
  end;
  
  u = cos( k * acos( x ) );
 
 
%-------------------------------------------------------------
% The exponential function via Taylor series
% x - independent variable
% n - order of expansion;

function [u] = ftaylor( x, n )

u = 1.0;
for k=1:n
  u = u + x.^k / factorial(k);
end;


%-------------------------------------------------------------
% Interpolated function
% x - independent variable

function [u] = f(x)

  u = exp( x );


%-------------------------------------------------------------
%ISPOSINT True for positive integer values. 
%       Mark Beale, 11-31-97
%       Copyright 1992-2001 The MathWorks, Inc.
% $Revision: 1.5 $
function f=isposint(v)
 
f = 1;
if ~isa(v,'double') | any(size(v) ~= [1 1]) | ...
  ~isreal(v) | v<0 | round(v) ~= v
  f = 0;
end
