% ### NDexample1.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; clf; % -------------- % 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); 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