%
% this program demonstrates issues with finite difference approximation.
% written by Ming Gu, Sept.2014 for MATH 228A
%

close all
%
% define a function 
%
func  = inline('sin(x)');
dfunc = inline('cos(x)');  % df/dx
ddfunc = inline('-sin(x)'); % d^2f/dx^2
xbar  = 1;

%
% forward, backward, central FD for dfunc.
%

h     = [1e-1 5e-2 1e-2 5e-3 1e-3];
fx    = func(xbar);
fxp   = func(xbar+h);
fxm   = func(xbar-h);
dfx   = dfunc(xbar);

DP    = (fxp-fx)./h;
DM    = (fx-fxm)./h;
D0    = (fxp-fxm)./(2*h);

EP    = abs(DP-dfx);
EM    = abs(DM-dfx);
E0    = abs(D0-dfx);

figure(1);clf;
loglog(h,EP,'bd-');
hold on
loglog(h,EM,'ro-.');
loglog(h,E0,'ks--');
legend('D+','D-','D0','Location','NorthWest');
title('Difference in Finite Difference Errors','FontSize',14);

%
% FD for ddfunc.
%

fx    = func(xbar);
ddfx  = ddfunc(xbar);
h     = [1e-1 5e-2 1e-2 5e-3 1e-3];
%
% hp != hm
%

hm    = h;
hp    = 2 *h;
fxp   = func(xbar+hp);
fxm   = func(xbar-hm);

DD1   = 2 *((fxp-fx)./hp - (fx-fxm)./hm)./(hp+hm);
E1    = abs(DD1-ddfx);

%
% hp = hm
%

hm    = h;
hp    = h;
fxp   = func(xbar+hp);
fxm   = func(xbar-hm);

DD2   = 2 *((fxp-fx)./hp - (fx-fxm)./hm)./(hp+hm);
E2    = abs(DD2-ddfx);

figure(2);clf;
loglog(h,E1,'bd-');
hold on
loglog(h,E2,'ks--');
legend('hp != hm','hp == hm','Location','NorthWest');
title('Difference in Finite Difference Errors','FontSize',14);


%
% forward, backward, central FD for dfunc. 
% same as before, but for much smaller h's
%

h     = [[1e-1 5e-2 1e-2 5e-3 1e-3], [1e-1 5e-2 1e-2 5e-3 1e-3]*1e-4, [1e-1 5e-2 1e-2 5e-3 1e-3]*1e-8];

fx    = func(xbar);
fxp   = func(xbar+h);
fxm   = func(xbar-h);
dfx   = dfunc(xbar);

DP    = (fxp-fx)./h;
DM    = (fx-fxm)./h;
D0    = (fxp-fxm)./(2*h);

EP    = abs(DP-dfx);
EM    = abs(DM-dfx);
E0    = abs(D0-dfx);

figure(3);clf;
loglog(h,EP,'bd-');
hold on
loglog(h,EM,'ro-.');
loglog(h,E0,'ks--');
legend('D+','D-','D0','Location','NorthWest');
title('Difference in Finite Difference Errors','FontSize',14);

%
% FD for ddfunc.
% same as before, but for much smaller h's
%

fx    = func(xbar);
ddfx  = ddfunc(xbar);
h     = [[1e-1 5e-2 1e-2 5e-3 1e-3], [1e-1 5e-2 1e-2 5e-3 1e-3]*1e-4, [1e-1 5e-2 1e-2 5e-3 1e-3]*1e-8];

%
% hp != hm
%

hm    = h;
hp    = 2 *h;
fxp   = func(xbar+hp);
fxm   = func(xbar-hm);

DD1   = 2 *((fxp-fx)./hp - (fx-fxm)./hm)./(hp+hm);
E1    = abs(DD1-ddfx);

%
% hp = hm
%

hm    = h;
hp    = h;
fxp   = func(xbar+hp);
fxm   = func(xbar-hm);

DD2   = 2 *((fxp-fx)./hp - (fx-fxm)./hm)./(hp+hm);
E2    = abs(DD2-ddfx);

figure(4);clf;
loglog(h,E1,'bd-');
hold on
loglog(h,E2,'ks--');
legend('hp != hm','hp == hm','Location','NorthWest');
title('Difference in Finite Difference Errors','FontSize',14);
DIR = '/home/mgu/Matlab/228A/';
for k = 1:4
   figure(k);
   print('-dpsc',[DIR,'FD',num2str(k),'.ps']);
   cmd = cat(2,'ps2pdf ', [DIR,'FD',num2str(k),'.ps'], ' ', [DIR,'FD',num2str(k),'.pdf']);
   [r,s]= system(cmd);
end
