% ### EXdiff1.m ###      09.04.14

% [source] http://ef.engr.utk.edu/ef230-2011-01/modules/matlab-integration/
% --> simple example of numerical differentiation compared to analytic solution

clear all;
% --------------
% User parameters
np = 30; % number of points on a curve
Nfact= 0.00; % scale factor for additive noise
fType= 0; % choose function: 0-sinusoid, 1-arctan
% --------------
% tStart= start point for time (ditto tEnd)
if fType==0
    func = @(t)4.*sin(t); % function
    dfunc = @(t)4.*cos(t); % derivative of function
    tStart= 0;  tEnd= 2*pi; loc= 'SouthWest';
elseif fType==1
    func = @(t)atan(t); % function
    dfunc = @(t)1./(1+t.^2); % derivative of function
    tStart= -2*pi;  tEnd= 2*pi; loc= 'NorthWest';
end
% ----
t = linspace(tStart,tEnd,np);  % determine time array
y = func(t)+ Nfact*randn(numel(t),1)';  % determined 'sampled' points (which can be noisy)
dydt = diff(y)./diff(t);    % compute derivative via numerical difference
% ----
% visualize
tt = linspace(tStart,tEnd,np*100);    % plot actual curve sans noise (oversample time here)
yy = func(tt);
figure(1); clf;
plot(tt,yy,'m-','LineWidth',2); hold on; grid on;
yy = dfunc(tt);
plot(tt,yy,'k--','LineWidth',2);
plot(t,y,'bs','LineWidth',2);    % plot points and connecting lines
% central difference - plot numerical derivative at midpoint
ttCD = t(1:end-1) + diff(t)./2;
plot(ttCD,dydt,'r+','LineWidth',2,'MarkerSize',8);
xlabel('t'); ylabel('y or dy/dt'); title('Example of numerical differentiation');
legend('y(t)','dy/dt','measured y(t) (w/ noise)','dy/dt (diff.m)','Location',loc);
% ----
% compare to Matlab's built-in gradient.m (which uses centered differences)
if 1==1
    slope= gradient(y,t);
    plot(t,slope,'o','LineWidth',2,'Color',[0 0.75 0])
    legend('y(t)','dy/dt','measured y(t) (w/ noise)','dy/dt (diff.m)','dy/dt (gradient.m)','Location',loc);
end