% Using a spectral Chebyshev-collocation method 
% to solve a second order elliptic problem in [-1.1]
% The equation
%            u'' = f   in [-1,1]
% BC
%            u(-1) = u(1) = 1
% solution
%            u(x) =e^(cos( u pi / 2))
% RHS 
%            chosen accordingly (see below)

function spectr_07
  
echo off, clc
clear all


Nmin = 2^3;
Nmax = 2^8;

u = Chebyshev_Collocation( Nmin, 1 ); hold on;
u = Chebyshev_Collocation( Nmax, 2 ); hold off;

% Studying convergence of the error;
N = Nmin;
i = 1;
while ( N <= Nmax );  
% Evaluating the solution at the centerpoint
% (which remains fixed regardless of discretization);
  u2 = Chebyshev_Collocation( N, 0 );
  Err2(i,1) = 2.0 * pi / N;
  Err2(i,2) = abs( u2(N/2+1)  - soln(0.0) ) / abs( soln(0.0) );
%  Err2(i,2) = abs( 0.5 * (u2(N/2+1) + u2(N/2)) - soln(0.0) ) / abs( soln(0.0) );
  N = N * 2;
  i = i + 1;
end;

% Plotting the error;
figure(3); loglog(Err2(:,1),Err2(:,2),'-o');
legend('Fourier-Galerkin','Fourier-Collocation'); 
title('Error at centerpoint');
xlabel('h'); ylabel('error');




%-------------------------------------------------------------
% Solver for the problem using Chebyshev-Collocation method
% N - number of grid-points;

function [Uc] = Chebyshev_Collocation( N, flag )  
  
% Generating the Gauss-Lobatto grid points
  X(1:N+1) = cos( (0:N) * pi / N  );

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

% Adding the boundary conditions;
  Fc(1)   = 1.0;
  Fc(N+1) = 1.0;

% Computing the differentiation matrix
  D = cheb( N );
  D2 = D * D;
  
% Incorporating the boundary conditions
  D2(1  ,:) = 0.0;   D2(  1,  1) = 1.0;
  D2(N+1,:) = 0.0;   D2(N+1,N+1) = 1.0;
  
% Solving the system;
  Uc(1:N+1) = D2\transpose(Fc);
  
% Exact solution (for comparison)  
  Ue = soln( X );
    
% Plotting the figure;  
  if ( flag == 1 )
    figure(2); plot(X,Uc,'-x');
  end;
  if ( flag == 2 )
    figure(2); plot(X,Uc,'-.',X,Ue,'-');
  end;

  
  
%-------------------------------------------------------------
% CHEB  compute D = differentiation matrix, x = Chebyshev grid
% Implemented by N.L. Trefethen

function [D,x] = cheb(N)
  if N==0, D=0; x=1; return, end
  x = cos(pi*(0:N)/N)'; 
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';                  
  D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
  D  = D - diag(sum(D'));                 % diagonal entries


  
  
%-------------------------------------------------------------
% The RHS function f(x), note that it depends on the global
% parameters [nu, a, b]  
function [rhs] = f( x )  

global nu; global a; global b;  
 
  rhs = pi^2 / 4.0 * exp(cos( x * pi / 2.0 )) .* ...
	( (sin( x * pi / 2.0 )).^2 - cos( x * pi / 2.0 ) );
  

%-------------------------------------------------------------
% Exact solution to the problem
function [u] = soln( x )  

  u = exp(cos( x * pi / 2.0 ));
  
  
%-------------------------------------------------------------
%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
