%THIS SCRIPT CALCULATES 3 DIFFERENCE APPROXIMATION METHODS AND COMPARES
%THEM TO THE ANALYTIC SOLUTION, INVESTIGATING THE ROOT MEAN SQUARE ERROR OF
%THE DIFFERENT APPROXIMATIONS.
% ALIANJECON OTERO
%MAY 29 2025
%%
clc, clear, close all;
% format strings
fs = 'ROOT MEAN SQUARE = %d';
fs2 = 'dx = %.2f';

%parameters
dx = 0.1;
x = 0:dx:2*pi;
N = size(x,2);
y = sin(x);

ypanalyticfda = cos(x(1:N-1));
ypanalyticbda = cos(x(2:N));
ypanalyticcda = cos(x(2:N-1));

%% FORWARD DIFFERENCE APPROXIMATION

%excludes the last point as there is no forward
ypfda = (y(2:N)-y(1:N-1))/dx;

figure(1);
plot(x,y,'k');
hold on;
plot(x(1:N-1),ypanalyticfda,'r');
plot(x(1:N-1),ypfda,'b');
axis([0 2*pi -1 1])
legend('sin(x)','Analytic Derivative','Forward Difference Approximation')

%calculate rmse
fdarmse = rmse(ypfda,ypanalyticfda);

text(0,-1.1,0,num2str(fdarmse,fs))
title(["Forward Difference Approximation",num2str(dx,fs2)])

%% BACKWARD DIFFERENCE APPROXIMATION

%excludes the first point as there is no backward
ypbda = (y(2:N)-y(1:N-1))/dx; 
%this is just the same as fda, just that the derivative of the ith point in fda
%would be the derivative of the ith+1 point in bda. you could absolutely
%remove this.

figure(2);
plot(x,y,'k');
hold on;
plot(x(2:N),ypanalyticbda,'r');
plot(x(2:N),ypbda,'b');
axis([0 2*pi -1 1])
legend('sin(x)','Analytic Derivative','Backward Difference Approximation')

%calculate rmse
bdarmse = rmse(ypbda,ypanalyticbda);

text(0,-1.1,0,num2str(bdarmse,fs))

title(["bda",num2str(dx,fs2)])

%% CENTRAL DIFFERENCE APPROXIMATION

%excludes the first and last point as there is no backward and forward
ypcda = (y(3:N)-y(1:N-2))/(2*dx);
%this method typically yields better results as it is able to sample and
%average two different directions. cda has an error proportional to delta x^2 compared
%to bda and fdas proportionality to delta x.

figure(3)
plot(x,y,'k');
hold on;
plot(x(2:N-1),ypanalyticcda,'r');
plot(x(2:N-1),ypcda,'b');
axis([0 2*pi -1 1])
legend('sin(x)','Analytic Derivative','Central Difference Approximation')

%calculate rmse
cdarmse = rmse(ypcda,ypanalyticcda);

text(0,-1.1,0,num2str(cdarmse,fs))

title(["Central Difference Approximation",num2str(dx,fs2)])

print(figure(1), 'fda.eps', '-depsc');
print(figure(2), 'bda.eps', '-depsc');
print(figure(3), 'cda.eps', '-depsc');