title: Why does my derivative function produce strange results?
body: <p>I am writing a small calculus library for personal use. In it are the standard calculus tools - taking the first derivative, nth derivative, Riemann sums, etc. One problem I have faced is that the nth derivative function is returning patently bogus results for certain values of h (the finite difference).</p>
  
  <p><a href="http://pastebin.com/6ttezSde" rel="nofollow">Code here</a> and below:</p>
  
  <pre><code>typedef double(*math_func)(double x);
  inline double max ( double a, double b ) { return a &gt; b ? a : b; }
  
  //Uses the five-point-stencil algorithm.
  double derivative(math_func f,double x){
      double h=max(sqrt(DBL_EPSILON)*x,1e-8);
      return ((-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h));
  }
  
  #define NDEPS (value)
  double nthDerivative(math_func f, double x, int N){
      if(N&lt;0) return NAN; //bogus value of N
      if(N==0) return f(x);
      if(N==1) return derivative(f,x);
      double* vals=calloc(2*N+9,sizeof(double)); //buffer region around the real values
      if(vals==NULL){ //oops! no memory
          return NAN;
      }
      int i,j;
      //don't take too small a finite difference
      double h=max(sqrt(DBL_EPSILON)*x,NDEPS);
      for(i=-(N+4);i&lt;=N+4;i++){
          vals[i+N+4]=derivative(f,x+h*i);
      }
      //for(i=0;i&lt;2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
      for(j=1;j&lt;N;j++){
          double *vals2=calloc(2*N+9,sizeof(double));
          for(i=2;i&lt;2*N+7;i++){
              vals2[i]=(-vals[i+2]+8*vals[i+1]-8*vals[i-1]+vals[i-2])/(12*h);
          }
          free(vals);
          vals=vals2;
          //for(i=0;i&lt;2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
      }
      double result=vals[N+4];
      free(vals);
      return result;
  }
  </code></pre>
  
  <p>A sample problem I give to test this function is the 5th derivative of sin(x) when x=pi. As we all know, the answer is -1. The problem comes when I try to vary the value of NDEPS ("Nth derivative epsilon").</p>
  
  <ul>
  <li>When NDEPS=1.5625e-2 (1/64.0): 5th derivative of sin() at x=pi:-1.0003e+00 (looks fine though).</li>
  <li>When NDEPS=1e-5 (1/100000.0): 5th derivative of sin() at x=pi:2.4302e+11 (I'm calling bullshit here).</li>
  </ul>
  
  <p>Why does this happen? Is it related to the nature of the sin() function? Or is it because of floating-point accuracy issues?</p>