%  Code to demonstrate derivative accuracy.
%
%  Using a periodic made up function 
%       u = -sin(x) + sin(2.x)^2*cos(x);
%  comparing its first derivative with numerical derivatives computed
%  using 5 difference methods:
% a) 1st order Backward                   du/dx = B(-1,1,0)/dx u
% b) 2nd order Central                    du/dx = B(-1,0,1)/(2 dx) u
% c) 2nd order Backward                   du/dx = B(1,-4,3,0,0)/(2 dx) u
% d) 4th order Central                    du/dx = B(1,-8,0,8,-1)/(12 dx) u
% e) 4th order Compact Pade   B(1,4,1)/6  du/dx = B(-1,0,1)/(2 dx) u
%
%   Code runs through a sequence of mesh refinements, halving dx with each.
%   Error is defined as the max norm of the difference between 
%   the exact derivative and numerical difference.

% Assignment:  Derive the Taylor Table Analysis for the Order of Accuracy 
%              for each of the methods used.  Compare the analysis results 
%              with the numerical results obtained below using up to 12 levels.
%              (If your machine is not powerful enough, use the maximum 
%               number of levels you can).
%              Write a short (1-2 paragraph) explaination of what you observed.
%
clear;

levels = input('Enter number of levels (>= 2) of grid refinement = ');

DX = zeros(100,1);
ER = zeros(100,6);

for l = 1:levels-2

jmax = 2^(l+2) + 1;
fprintf(' jmax = %5.0f \n',jmax);

% initialize array sizes
u = zeros(jmax,2);
d = zeros(jmax,6);
x = zeros(jmax,1);
e = zeros(jmax,6);

% Grid spacing
dx = 2.*pi/(jmax);
x = dx*[0:jmax-1]';
fprintf(' dx =  %16.8e \n',dx);
DX(l,1) = dx;

% Function and derivative
for j=1:jmax
  arg = x(j,1);
% Function  
  u(j,1) = -sin(arg) + sin(2.*arg)^2*cos(arg);
% 1st Derivative
  u(j,2) = -cos(arg) + 4.0*sin(2.*arg)*cos(2*arg)*cos(arg) - ...
      sin(2*arg)^2*sin(arg);
end

figure(1);
clf;
subplot(3,1,1);
plot(x,u(1:jmax,1));
title('Function');
subplot(3,1,2);
plot(x,u(1:jmax,2),'r');
title('Derivatives');

% a) 1st order Backward   du/dx = B(-1,1,0)/dx u
d(2:jmax,1) = (u(2:jmax,1)-u(1:jmax-1,1))/dx;
% Boundary
d(1,1) = (u(1,1) - u(jmax,1))/dx;

% b) 2nd order Central    du/dx = B(-1,0,1)/(2 dx) u
d(2:jmax-1,2) = 0.5*(u(3:jmax,1)-u(1:jmax-2,1))/dx;
% Boundaries
d(1,2) = 0.5*(u(2,1) - u(jmax,1))/dx;
d(jmax,2) = 0.5*(u(1,1) - u(jmax-1,1))/dx;

% c) 2nd order Backward   du/dx = B(3,-4,1)/(2 dx) u
d(3:jmax,3) = 0.5*(3.*u(3:jmax,1)-4.*u(2:jmax-1,1) + u(1:jmax-2,1))/dx;
d(1,3) = 0.5*(3.*u(1,1)-4.*u(jmax,1) + u(jmax-1,1))/dx;
d(2,3) = 0.5*(3.*u(2,1)-4.*u(1,1) + u(jmax,1))/dx;

% d) 4th order Central    du/dx = B(-1,8,0,-8,1)/(12 dx) u
d(3:jmax-2,4) = (-u(5:jmax,1) + 8.*u(4:jmax-1,1) ...
    -8.*u(2:jmax-3,1)+u(1:jmax-4,1) )/dx/12.;
% Boundaries
d(1,4) = (-u(3,1) + 8.*u(2,1) ...
    -8.*u(jmax,1)+u(jmax-1,1) )/dx/12.;
d(2,4) = (-u(4,1) + 8.*u(3,1) ...
    -8.*u(1,1)+u(jmax,1) )/dx/12.;
d(jmax-1,4) = (-u(1,1) + 8.*u(jmax,1) ...
    -8.*u(jmax-2,1)+u(jmax-3,1) )/dx/12.;
d(jmax,4) = (-u(2,1) + 8.*u(1,1) ...
    -8.*u(jmax-1,1)+u(jmax-2,1) )/dx/12.;

% e) 4th order Compact Pade   B(1,4,1)/6  du/dx = B(-1,0,1)/(2 dx) u
a = ones(jmax,1)/6.0;
b = 4.*a;
c = a;
q = zeros(jmax,1);
s = zeros(jmax,1);
d(1:jmax,5) = d(1:jmax,2);
% Periodic tridiagaonal sover used.
d(1:jmax,5) = trip(jmax,a,b,c,d(1:jmax,5),q,s,1,jmax);

hold on;
% Plot derivatives
plot(x,d(1:jmax,1),'r-.');
plot(x,d(1:jmax,2),'g:');
plot(x,d(1:jmax,3),'c--');
plot(x,d(1:jmax,4),'y-');
plot(x,d(1:jmax,5),'m-');

%  Compute local error
e(1:jmax,1) = abs(d(1:jmax,1) - u(1:jmax,2));
e(1:jmax,2) = abs(d(1:jmax,2) - u(1:jmax,2));
e(1:jmax,3) = abs(d(1:jmax,3) - u(1:jmax,2));
e(1:jmax,4) = abs(d(1:jmax,4) - u(1:jmax,2));
e(1:jmax,5) = abs(d(1:jmax,5) - u(1:jmax,2));

% Plot local errors.
subplot(3,1,3);
semilogy(x,e(1:jmax,1),'r-.');
hold on;
semilogy(x,e(1:jmax,2),'g:');
semilogy(x,e(1:jmax,3),'c--');
semilogy(x,e(1:jmax,4),'y-.');
semilogy(x,e(1:jmax,5),'m--');
title('Errors');

%  Max Norm of error
ER(l,1) = norm(e(1:jmax,1),inf);
ER(l,2) = norm(e(1:jmax,2),inf);
ER(l,3) = norm(e(1:jmax,3),inf);
ER(l,4) = norm(e(1:jmax,4),inf);
ER(l,5) = norm(e(1:jmax,5),inf);

%%pause;
end

% Plot total Error in log-log.  Slope of curves is error order measure.
figure(2);
clf;
loglog(DX(1:levels-2,1),ER(1:levels-2,1),'r+');
loglog(DX(1:levels-2,1),ER(1:levels-2,1),'r--');
hold on;
loglog(DX(1:levels-2,1),ER(1:levels-2,2),'gx');
loglog(DX(1:levels-2,1),ER(1:levels-2,2),'g--');
loglog(DX(1:levels-2,1),ER(1:levels-2,3),'c*');
loglog(DX(1:levels-2,1),ER(1:levels-2,3),'c--');
loglog(DX(1:levels-2,1),ER(1:levels-2,4),'yo');
loglog(DX(1:levels-2,1),ER(1:levels-2,4),'y--');
loglog(DX(1:levels-2,1),ER(1:levels-2,5),'m+');
loglog(DX(1:levels-2,1),ER(1:levels-2,5),'m--');
loglog(DX(1:levels-2,1),DX(1:levels-2,1),'w');
loglog(DX(1:levels-2,1),DX(1:levels-2,1).^2,'w');
loglog(DX(1:levels-2,1),DX(1:levels-2,1).^4,'w');
xlabel('dx');
ylabel('Er');