% Analyzing the accuracy of the finite-difference solutions
% to a boundary value problem
% Problem solved:
% y''(x) = g(x)    for  0 < x < 2*pi
% y(0) = y(2*pi) = 1  (Dirichlet problem)
% y'(0) = y'(2*pi) = 0 (Neumann problem)

function ode_01
  
echo off, clc
clear all
 

x0 = pi;
Nmin = 1e1 + 1;
Nmax = 1e4;

N  = Nmin;
g0 = cos(x0);
i = 1;

% Node indexing convention:
% i=0       - left endpoint
% i=1,...,N - interior nodes
% i=N+1     - right endpoint

while ( N <= Nmax );  
  h = 2.0 * pi / (N + 1);
  g1 = cos( h );
% Dirichlet problems  
  [D1a D2a] = BWP1( N, h );
  [D1b D2b] = BWP1c( N, h );  
% Neumann problems  
%  [N1a N2a] = BWP2( N, h );
  [N1a N2a] = BWP2a( N, h );
  [N1b N2b] = BWP2b( N, h );
  Err1(i,1) = h;
  Err2(i,1) = h;
  Err1(i,2) = abs(g0 - D2a) / abs( g0 );
  Err1(i,3) = abs(g0 - D2b) / abs( g0 );
  Err1(i,4) = abs(g0 - N2a) / abs( g0 );
  Err1(i,5) = abs(g0 - N2b) / abs( g0 );
%  Err1(i,4) = abs(g0 - D2c) / abs( g0 );
  Err2(i,2) = abs(g1 - D1a) / abs( g1 );
  Err2(i,3) = abs(g1 - D1b) / abs( g1 );
  Err2(i,4) = abs(g1 - N1a) / abs( g1 );
  Err2(i,5) = abs(g1 - N1b) / abs( g1 );
%  Err2(i,4) = abs(g1 - D1c) / abs( g1 );
  i = i + 1;
  N = N * 2 + 1;
end;

%figure(1); loglog(Err1(:,1),Err1(:,2),'-x',Err1(:,1),Err1(:,3),'-o');
figure(1); loglog(Err1(:,1),Err1(:,2),'-x',Err1(:,1),Err1(:,3),'-o',Err1(:,1),Err1(:,4),'-^',Err1(:,1),Err1(:,5),'-s');
legend('Dirichlet','Dirichlet - compact','Neumann - 1st ord','Neumann - 2nd ord'); 
title('Central Node');
xlabel('h'); ylabel('error');

%figure(2); loglog(Err2(:,1),Err2(:,2),'-x',Err2(:,1),Err2(:,3),'-o');
figure(2); loglog(Err2(:,1),Err2(:,2),'-x',Err2(:,1),Err2(:,3),'-o',Err2(:,1),Err2(:,4),'-^',Err2(:,1),Err2(:,5),'-s');
legend('Dirichlet','Dirichlet - compact','Neumann - 1st ord','Neumann - 2nd ord'); 
title('Endpoint');
xlabel('h'); ylabel('error');


%-------------------------------------------------------------
% Solves the two-point Dirichlet boundary value problem using
% a second-order central finite-difference scheme

function [x1, x2] = BWP1( N, h )  
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
% the main diagonal;    
    b(i) = -2.0;
    
% the upper diagonal;    
    if ( i < N )
      c(i) = 1.0;
    end;
       
% the lower diagonal;    
    if ( i < N )
      a(i) = 1.0;
    end;
    
    if ( i > 1 & i < N )
% RHS terms in the interior;     
      rhs(i) = - h^2 * cos( i * h );
    else
% endpoints
      rhs(i) = - h^2 * cos( i * h ) - 1.0;
    end;
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  x2 = X( ceil(N/2) );
  x1 = X( 1 );
  
  
%-------------------------------------------------------------
% Solves the two-point Dirichlet boundary value problem using
% a second-order COMPACT central finite-difference scheme
  
function [x1, x2] = BWP1c( N, h )  
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
% the main diagonal;    
    b(i) = -2.0;
    
% the upper diagonal;    
    if ( i < N )
      c(i) = 1.0;
    end;
       
% the lower diagonal;    
    if ( i < N )
      a(i) = 1.0;
    end;
    
    if ( i > 1 & i < N )
% RHS terms in the interior;     
      rhs(i) = - h^2 * (cos( (i + 1) * h )+10.0*cos( i * h )+cos( (i - 1) * h ))/12.0;
    else
% endpoints
      rhs(i) = - h^2 * (cos( (i + 1) * h )+10.0*cos( i * h )+cos( (i - 1) * h ))/12.0 - 1.0;
    end;
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  x2 = X( ceil(N/2) );
  x1 = X( 1 );

  
%-------------------------------------------------------------
% Solves the two-point Neumann boundary value problem using
% a second-order central finite-difference scheme with
% first-order approximation at endpoints;
% WARNING - The system matrix is singular!
function [x1, x2] = BWP2( N, h )  
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
% the main diagonal;    
    if ( i > 1 & i < N )
      b(i) = -2.0;
    else
      b(i) = -1.0;
    end;
    
% the upper diagonal;    
    if ( i < N )
      c(i) = 1.0;
    end;
       
% the lower diagonal;    
    if ( i < N )
      a(i) = 1.0;
    end;
    
    rhs(i) = - h^2 * cos( (i - 1) * h );    
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  x2 = X( ceil(N/2) );
  x1 = X( 1 );
    
  

%-------------------------------------------------------------
% Solves the two-point Neumann boundary value problem using
% a second-order central finite-difference scheme with
% first-order approximation at endpoints;
% System matrix is desingularized;
function [x1, x2] = BWP2a( N, h )  
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
% the main diagonal;    
    if ( i > 1 & i < N )
      b(i) = -2.0;
    else
      b(i) = -1.0;
    end;
    
% the upper diagonal;    
    if ( i < N )
      c(i) = 1.0;
    end;
       
% the lower diagonal;    
    if ( i < N - 1 )
      a(i) = 1.0;
    elseif ( i == N - 1 )
      a(i) = 0.0;
    end;
    
    if ( i < N )
% RHS terms in the interior and left endpoint;     
      rhs(i) = - h^2 * cos( (i - 1) * h );
    else
% Right endpoint - need to desingularize the system;
      rhs(i) = - cos( h );
    end;
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  x2 = X( ceil(N/2) );
  x1 = X( 1 );


%-------------------------------------------------------------
% Solves the two-point Neumann boundary value problem using
% a second-order central finite-difference scheme with
% forward/backward second-order approximation at endpoints;
% System matrix is desingularized;
function [x1, x2] = BWP2b( N, h )  
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
% the main diagonal;    
    if ( i == 1 )
      b(i) = -2.0/3.0;
    elseif ( i < N )
      b(i) = -2.0;
    else
      b(i) =  1.0;
    end;
    
% the upper diagonal;    
    if ( i == 1 )
      c(i) = 2.0 / 3.0;
    elseif ( i < N )
      c(i) = 1.0;
    end;
    
% the lower diagonal;    
    if ( i < N - 1 )
      a(i) = 1.0;
    elseif ( i == N - 1 )
      a(i) = 0.0;
    end;
    
    if ( i < N )
% RHS terms in the interior and left endpoint;     
      rhs(i) = - h^2 * cos( i * h );
    else
% Right endpoint - need to desingularize the system;
      rhs(i) = cos( h );
    end;
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  x2 = X( ceil(N/2) );
  x1 = X( 1 );
  

  
%      x = tridiag1( a, b, c, y )
%
%  Solve the  N x N  tridiagonal system for x:
%
%  [ b(1)  c(1)                                  ] [  x(1)  ]   [  y(1)  ] 
%  [ a(1)  b(2)  c(2)                            ] [  x(2)  ]   [  y(2)  ] 
%  [       a(2)  b(3)  c(3)                      ] [        ]   [        ] 
%  [            ...   ...   ...                  ] [  ...   ] = [  ...   ]
%  [                    ...    ...    ...        ] [        ]   [        ]
%  [                        a(N-2) b(N-1) c(N-1) ] [ x(N-1) ]   [ y(N-1) ]
%  [                               a(N-1)  b(N)  ] [  x(N)  ]   [  y(N)  ]
%
%  y must be a vector (row or column); N is determined from its length.
%  a, b, c must be vectors of lengths at least N-1, N, and N-1 respectively,
function x = tridiag1( a, b, c, y )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Check that the input arrays have acceptable sizes

if min(size(y))~=1 | min(size(a))~=1 | min(size(b))~=1 | min(size(c))~=1
   error('a, b, c, y must be vectors');
end
N = length(y);
if length(a)<N-1 | length(b)<N | length(c)<N-1
   error('a, b, c must be vectors of length at least N-1, N, N-1');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solve the problem by back substitution

gam = zeros(1,N);   %  hold the LU decomposition
x = zeros(size(y));

%  phase 1:  LU decomposition
beta = b(1);
if beta==0
   error('beta = 0 at j=1:  matrix is singular');
end

x(1) = y(1) / beta;
for j=2:N
   gam(j) = c(j-1) / beta;
   beta = b(j) - a(j-1)*gam(j);
   if beta==0
      error( [ 'beta = 0 at j=' num2str(j) ':  matrix is singular' ]);
   end
   x(j) = ( y(j) - a(j-1)*x(j-1) )/beta;
end

%  phase 2:  back-substitution
for j=N-1:-1:1
   x(j) = x(j) - gam(j+1)*x(j+1);
end
  
  
