%% MATH 745 (Fall 2023)
%% Solution to Homework Assignment #1

function hw1
  
echo off, clc
clear all


%% QUESTION 1 

x0 = 0.0;
hmax = 0.1;
hmin = 1.0e-10;

x0 = 0.0;
h = hmax;
df0 = df( x0 );
i = 1;
while ( h >= hmin );  
  Err(i,1) = h;
  Err(i,2) = abs( (df0 - FD( x0, h )) );
  Err(i,3) = abs( (df0 - FD2( x0, h )));
  Err(i,4) = abs( (df0 - CD( x0, h )) );
  Err(i,5) = abs( (df0 - CD4( x0, h )) );
  Err(i,6) = abs( (df0 - CSD( x0, h )) );
  i = i + 1;
  h = h / 2.0;
end;



figure(1); loglog(Err(:,1),Err(:,2),'-x',Err(:,1),Err(:,3),'-+',Err(:,1),Err(:,4),'-o',Err(:,1),Err(:,5),'-d',Err(:,1),Err(:,6),'-s');
legend('1st ord Forward Difference','2nd ord Forward Difference', ...
       '2nd ord Central Difference', '4th ord Central Difference',...
       '2nd ord Complex Step');
title('Question 1(a), x=0.0')


x0 = 1.0;
h = hmax;
df0 = df( x0 );
i = 1;
while ( h >= hmin );  
  Err(i,1) = h;
  Err(i,2) = abs( (df0 - FD( x0, h )) / df0 );
  Err(i,3) = abs( (df0 - FD2( x0, h )) / df0 );
  Err(i,4) = abs( (df0 - CD( x0, h )) / df0 );
  Err(i,5) = abs( (df0 - CD4( x0, h )) / df0 );
  Err(i,6) = abs( (df0 - CSD( x0, h )) / df0 );
  i = i + 1;
  h = h / 2.0;
end;


figure(2); loglog(Err(:,1),Err(:,2),'-x',Err(:,1),Err(:,3),'-+',Err(:,1),Err(:,4),'-o',Err(:,1),Err(:,5),'-d',Err(:,1),Err(:,6),'-s');
legend('1st ord Forward Difference','2nd ord Forward Difference', ...
       '2nd ord Central Difference', '4th ord Central Difference',...
       '2nd ord Complex Step');
title('Question 1(b), x=1.0')






%% QUESTION 2

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

N  = Nmin;
i = 1;

global CASE

% 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;
  xx = linspace(0.0, 2*pi,  N + 1);

% solution for the first RHS function  
  CASE = 1;  
  g = soln( xx );
% Dirichlet problems  
  y1 = Solve_Dirichlet_2ord( N-1, h );
  y2 = Solve_Dirichlet_compact( N-1, h );
  Err1(i,1) = h;
  Err1(i,2) = norm(y1 - g(2:N),inf);
  Err1(i,3) = norm(y2 - g(2:N),inf);

% solution for the first RHS function  
  CASE = 2; 
  g = soln( xx );
% Dirichlet problems  
  y1 = Solve_Dirichlet_2ord( N-1, h );
  y2 = Solve_Dirichlet_compact( N-1, h );
  Err1(i,4) = norm(y1 - g(2:N),inf);
  Err1(i,5) = norm(y2 - g(2:N),inf);
 
%  figure(100)
%  plot(xx,[-1 y1 1],'-x', xx,[-1 y2 1],':o',xx,soln(xx),'-r')
   
  i = i + 1;
  N = N * 2 + 1;
end;

figure(3); 
loglog(Err1(:,1),Err1(:,2),'-xr',Err1(:,1),Err1(:,3),'-ob');
legend('Dirichlet','Dirichlet - compact'); 
title('Central Node');
xlabel('h'); ylabel('error');
title('Question  2 - CASE 1')


figure(4); 
loglog(Err1(:,1),Err1(:,4),'-xr',Err1(:,1),Err1(:,5),'-ob');
legend('Dirichlet','Dirichlet - compact'); 
title('Central Node');
xlabel('h'); ylabel('error');
title('Question  2 - CASE 2')


%% QUESTION 3

disp('  Question #3')
disp('   ')

T = 1.0;

h = T / 100;
dh = 0.05*h;

Y = zeros(100,2);

tt = linspace(0,T,100);
figure(5);
plot(tt,10.0 ./(1.0+200.0*tt).^0.5,'-r')
hold on 

% We integrate the ODE repeatedly with increasing step sizes h
% and check when the numerical solution is no longer decreasing 
% (which is an indication of a numerical instability)
while ( (h < T) & (Y(1,2) >= Y(end,2)) ) 
  
  Y = RK4( h, T );

  fprintf('Step size h = %g \n',h)

  plot(Y(:,1),Y(:,2),':xb')

  h = h + dh;
  
end

fprintf('\nLargest Possible Step size h = %g \n',h-2*dh)


figure(5)
axis([0.0 T  0.0 12.0 ]);
xlabel('t'), ylabel('y')
title('Question 3')


% Implementation of the functions

% functions used in Question #1

% The function to be differentiated;
function [f] = f( x )

  if ( x < 0.0 )
    f = x.^2;
  else
    f = x.^4;
  end;


% The exact derivative of the function f(x)
function [df] = df( x )

  if ( x <= 0.0 )
    df = 2.0*x;
  else
    df = 4.0*x.^3;
  end;  

  
% First-order forward difference formula
function [fd] = FD( x0, h )  
  
  fd = ( f(x0 + h) - f(x0) ) / h;

  
% Second-order forward difference formula  
function [fd] = FD2( x0, h )  
  
  fd = ( - f(x0 + 2.0*h) + 4.0*f(x0 + h) - 3.0*f(x0) ) / (2.0*h);

  
% Second-order central difference formula
function [fd] = CD( x0, h )  
  
  fd = ( f(x0 + h) - f(x0 - h) ) / (2.0 * h);
  
  
% Fourth-order central difference formula  
function [fd] = CD4( x0, h )  
  
  fd = ( - f(x0 + 2.0*h) + 8.0*f(x0 + h) - 8.0*f(x0 - h) + f(x0 - 2.0*h) ) / (12.0*h);
  
% Complex Step Derivative
function [fd] = CSD( x0, h )  
  
  fd = imag( f(x0 + j*h)  ) / h;

% functions used in Question #2

% Exact solutions
function [y] = soln(x)
global CASE  
  
  switch CASE 
   case 1,  
    y =  (1/60)*(-3*x.^5*pi-20*x.^3*pi-30*x.^2*pi+60*x+48*x*pi^5+80*x*pi^3+60*x*pi^2-60*pi)/pi;
   case 2,
    y =  (1/30)*(-x.^6*pi-10*x.^3*pi-15*x.^2*pi+30*x+32*x*pi^6+40*x*pi^3+30*x*pi^2-30*pi)/pi;
  end;

%  The RHS function in Question 3
function [y] = rhs(x)
global CASE  
  
  switch CASE 
   case 1,
    y = x^3+2*x+1;
   case 2,
    y = x^4+2*x+1;
  end;

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

function [y] = Solve_Dirichlet_2ord( 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 and left endpoint;     
      r(i) = - h^2 * rhs( i * h );
    elseif ( i == 1 )
% left endpoint 
      r(i) = - h^2 * rhs( i * h ) + 1.0;
    elseif ( i == N )
% right endpoint 
      r(i) = - h^2 * rhs( i * h ) - 1.0;
    end;
  end;
  
% Solving the system;
  y = tridiag1( a, b, c, r );
  
  
%-------------------------------------------------------------
% Solves the two-point Dirichlet boundary value problem using
% a second-order COMPACT central finite-difference scheme
  
function [y] = Solve_Dirichlet_compact( 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 and left endpoint;     
      r(i) = - h^2 * (rhs( (i + 1) * h )+10.0*rhs( i * h )+rhs( (i - 1) * h ))/12.0;
    elseif ( i == 1 )
% left endpoint 
      r(i) = - h^2 * (rhs( (i + 1) * h )+10.0*rhs( i * h )+rhs( (i - 1) * h ))/12.0 + 1.0;
    elseif ( i == N )
% right endpoint 
      r(i) = - h^2 * (rhs( (i + 1) * h )+10.0*rhs( i * h )+rhs( (i - 1) * h ))/12.0 - 1.0;
    end;
  end;
  
% Solving the system;
  y = tridiag1( a, b, c, r );
  
  
%-------------------------------------------------------------
% 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;
    
    r(i) = - h^2 * rhs( (i - 1) * h );    
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, r );
  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
    


% functions used in Question #3

% the RK4 integrator
function [Y] = RK4( h, T )
  

i = 1;
t = i * h;
Yrn = 10.0;    % the initial condition
while ( t < T )
  t = i * h;
  k1r = F(Yrn);
  
  k2r = F(Yrn + 0.5 * h * k1r);
  
  k3r = F(Yrn + 0.5 * h * k2r);
  
  k4r = F(Yrn + h * k3r);
  
  Yrn1 = Yrn + h * ( k1r/6.0 + (k2r + k3r)/3.0 + k4r/6.0 );

  Yrn = Yrn1;
  
  Y(i,1) = t;
  Y(i,2) = Yrn1;
  i = i + 1;
end;


% the RHS of the ODE
function [fy] = F(y)
  
  fy = - y.^3;
  
  
