% Using a spectral Fourier collocation method to solve a second
% order elliptic problem *with variable coefficients*
% in a 2pi-periodic domain

function spectr_04
  
echo off, clc
clear all

global nu; global b;  

nu = 1.0;
b  = 5.0;

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

u = Fourier_Collocation( Nmin, 1 ); hold on;
u = Fourier_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 = Fourier_Collocation( N, 0 );
  Err2(i,1) = 2.0 * pi / N;
  Err2(i,2) = abs( u2(N/2+1) - soln(pi) ) / abs( soln(pi) );
  N = N * 2;
  i = i + 1;
end;

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



%-------------------------------------------------------------
% Solver for the problem using Fourier-Collocation method
% N - number of grid-points (power of 2)

function [Ur] = Fourier_Collocation( N, flag )  

global nu; global b;  
  
% Checking of the number of grid points is consistent with FFT;
  if ( ~isposint(log2( N )) )
    disp(' Wrong number of grid-points!')
    return;
  end;
    
% Setting up the grid
  dx = 2.0 * pi / N;
  X = (0.0:dx:(N-1)*dx);
  
  
% Assembling the differentiation matrices;
  for p=1:N
    for q=1:N
      h = (p - q) * dx / 2.0;
      if ( p ~= q )
	D1(p,q) = 0.5 * (-1)^(p+q) * cot( h );
	D2(p,q) = 0.25 * (-1)^(p+q) * N + 0.5 * (-1)^(p+q+1) / (sin(h))^2;
      else
	D1(p,q) = 0.0;
	D2(p,q) = - ( N - 1 ) * ( N - 2 ) / 12.0;
      end;
    end;
% Multiplying (row-wise) the differentiation matrix D1 by 
% the space-dependent coefficient;
    aD1(p,:) = a((p-1)*dx) * D1(p,:);
  end;  
  

% Assembling the system matrix;  
  A = nu * D2 - aD1 + b * eye(N);

% Evaluating the RHS function;
  Fr = f( X );
  
% Solving the system;
  Ur(1:N) = A\transpose(Fr);

% Exact solution (for comparison)  
  Ue = soln( X );
  
% Plotting the figure;  
  if ( flag == 1 )
    figure(2); plot(X,Ur,'-x');
  end;
  if ( flag == 2 )
    figure(2); plot(X,Ur,'-.',X,Ue,'-');
  end;


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

global nu; global b;  
 
  rhs = exp(sin(x)) .* ( nu * ((cos(x)).^2 - sin(x)) - a(x) .* cos(x) + b );
  

%-------------------------------------------------------------
% Function defining the space-dependent coefficient 'a'
function [u] = a( x )  

  u = (sin(x)).^2;


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

  u = exp(sin(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
