%
% This simple program approximates first- and second-order spatial derivatives using
% finite difference schemes.
% The function considered is exp(x) and we compute derivatives at x=1.0
% Order of convergence of the following schemes are studied:
% 1. Forward differece for first derivative
% 2. Backward differece for first derivative
% 3. 2nd order central differece for first derivative
% 4. 2nd order central differece for second derivative
%
% This is a part of teaching materials provided for the course
% "Computational Heat & Fluid Flow (ME 605)" taught at IIT Goa
% during the Spring semester 2020
%
% Written by: Y Sudhakar, IIT Goa
% email: sudhakar@iitgoa.ac.in
%
% Tested with Matlab-R2018b on 15 January 2020
%
% Assignment for the students
% 1. second order one-sided difference for first derivative
% 2. Fourth order central difference approximation for first derivative
% 3. Fourth order central difference approximation for second derivative
%
clc;
clear all;

num_refine = 10;  % number of refinement level considered to plot order of convergence

h      = ones(num_refine,1);      % mesh spacing for each refinement
bd     = ones(num_refine,1);      % approximation by backward differencing
bderr  = ones(num_refine,1);      % error introduced in backward differencing
fd     = ones(num_refine,1);      % approximation by forward differencing
fderr  = ones(num_refine,1);      % error introduced in forward differencing
cd     = ones(num_refine,1);      % approximation by central differencing
cderr  = ones(num_refine,1);      % error introduced in central differencing
cd2    = ones(num_refine,1);      % approximation of second derivative by central differencing
cd2err = ones(num_refine,1);      % error introduced in central differencing for second derivative

% exact value at x=1.0
exac = exp(1.0)

for i=1:num_refine
    h(i)=0.1/(2^(i-1));
    xp=1.0;              % location and solution value at the considered node
    fp=exp(xp);
    xw=1.0-h(i);       % left side node (denoted with 'w' meaning 'west' as in finite volume context)
    fw=exp(xw);
    xe=1.0+h(i);       % right side node (denoted with 'e' meaning 'east' as in finite volume context)
    fe=exp(xe);
    bd(i)=(fp-fw)/h(i);         % backward difference approximation and error
    bderr(i)=abs(exac-bd(i));
    fd(i)=(fe-fp)/h(i);         % forward difference approximation and error
    fderr(i)=abs(exac-fd(i));
    cd(i)=0.5*(fe-fw)/h(i);     % central difference approximation and error
    cderr(i)=abs(exac-cd(i));
    cd2(i)=(fw+fe-2.0*fp)/h(i)/h(i);  % second order derivatives using central differencing
    cd2err(i)=abs(exac-cd2(i));
end

% print(exac)
% print(fd, bd, cd)

% Lines showing first and second order convergence
fir    = ones(num_refine,1);
sec    = ones(num_refine,1);
fou    = ones(num_refine,1);
fir(1) = 0.05;
sec(1) = 0.05; %0.01
fou(1) = 0.05;
i = 2;
while i <= num_refine
    fir(i) = fir(i-1)*(h(i)/h(i-1));
    sec(i) = sec(i-1)*((h(i)/h(i-1))^2);
    fou(i) = fou(i-1)*((h(i)/h(i-1))^4);
    i = i+1;
end

% plot all quantities
set(gca,'TickLabelInterpreter','latex')
loglog(h,bderr);
hold on;
loglog(h,fderr,'o');
hold on;
loglog(h,cderr);
hold on;
loglog(h,cd2err);
hold on;
loglog(h,fir,'--');
hold on;
loglog(h,sec,'--');
hold on;
loglog(h,fou,'--');
legend('Backward','Forward','Central','Central 2nd der.','1st order','2nd order','4th order');
xlabel('${\Delta x}$','interpreter','latex', 'FontWeight','bold');
ylabel('Truncation Error');
