// numdiff : compare numerical derivatives of cos(x) in different approximations
//    show, how the error decreases with decreasing data distance 
// (corresponding to increasing number of data points)

function [a,b,c]=numdiff(n) //n= number of data points
x=0:%pi/n:%pi-0.00001 ;   // Interval over which we test the approximations error

y=sin(x); // the function itself

exact=cos(x); // its exact derivative
fwd=exact;   // initialise forward difference approx.
mid=exact;   // initialise symmetric 
fifth=exact; // initialise fifth order approx

for i=1:n-1
 fwd(i)=(y(i+1)-y(i))/(x(i+1)-x(i)); // forward difference
end

for i=2:n-1
 mid(i)=(y(i+1)-y(i-1))/(x(i+1)-x(i-1)); // central difference, denominator = 2*h
end

for i=3:n-2
 fifth(i)=(y(i-2)-8*y(i-1)+8*y(i+1)-y(i+2))/(12*(x(i+1)-x(i))); // fifth order Lagr., valid for equidistant data!
end

//plot(x,exact,"g-");  // uncomment these lines to show more
//plot(x,fwd,"ro");
//plot(x,mid,"bo");
//plot(x,fifth,"b+");

a=max(abs(fwd-exact)/exact);
b=max(abs(mid-exact)/exact);
c=max(abs(fifth-exact)/exact);

endfunction


// now use that function for testing different number of points n, decreasing h
a=0;b=0;c=0;
n=5;
for i=1:9;
 [a(i),b(i),c(i)]=numdiff(n);   // call the function to get back numerical derivatives maximal errors
 nn(i)=n;                       // for how many data points was this result?
 n=n*2; // in each cycle of the for-loop, increase n by a factor 2
end

clf();
plot(log10(nn),log10(b),"g+")  // log-log-plot of relative max. errors against number of data points n
plot(log10(nn),log10(c),"rx")

plot(log10(nn),log10(a),"b.")

reglin(log10(nn)',log10(a)')  // determine the 3 slopes in the plots
reglin(log10(nn)',log10(b)')
reglin(log10(nn(2:9))',log10(c(2:9))') // leave out first point, take only points 2...9.

// END OF CODE,   (C) Stephan Frickenhaus