Why does my derivative function produce strange results?

Question

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).

Code here and below:

typedef double(*math_func)(double x);
inline double max ( double a, double b ) { return a > 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<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<=N+4;i++){
        vals[i+N+4]=derivative(f,x+h*i);
    }
    //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    for(j=1;j<N;j++){
        double *vals2=calloc(2*N+9,sizeof(double));
        for(i=2;i<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<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    }
    double result=vals[N+4];
    free(vals);
    return result;
}

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").

• When NDEPS=1.5625e-2 (1/64.0): 5th derivative of sin() at x=pi:-1.0003e+00 (looks fine though).
• When NDEPS=1e-5 (1/100000.0): 5th derivative of sin() at x=pi:2.4302e+11 (I'm calling bullshit here).

Why does this happen? Is it related to the nature of the sin() function? Or is it because of floating-point accuracy issues?

Answer 1

Here's an adaptation of your code. The main modification, apart from using more white space, is to make NDEPS into an argument to the nthDerivative() function, so that it can be called with different values and to add copious printing. I've also had to get inventive on the plain derivative() function; the code compiles correctly (but I'm really not trying to make a philosophical or theological statements with the assertion assert(sin == fun);, but it does mean the code compiles without warnings and it recognizes the limitations of this derivative function).

Code

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define max(a, b)    (((a) > (b)) ? (a) : (b))

#define PRIe_double    "%21.15e"

typedef double(*math_func)(double x);

static double derivative(math_func fun, double x) { assert(sin == fun); return cos(x); }

static double nthDerivative(math_func f, double x, int N, double NDEPS)
{
    if (N < 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);
    printf("h = " PRIe_double "\n", h);

    for (i = -(N+4); i <= N+4; i++)
    {
        vals[i+N+4] = derivative(f, x+h*i);
        printf("%2d: deriv(" PRIe_double ") = " PRIe_double "\n", i, x+h*i, vals[i+N+4]);
    }

    for (j = 1; j < N; j++)
    {
        printf("Iteration %d\n", j);
        double *vals2 = calloc(2*N+9, sizeof(double));
        for (i = 2; i < 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 < 2*N+9; i++)
            printf("%2d: " PRIe_double "\n", i, vals[i]);
    }
    double result = vals[N+4];
    free(vals);
    return result;
}

int main(void)
{
    double val = M_PI;
    double eps;
    double r;

    eps = 1.0 / 64.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    eps = 1.0 / 100000.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    return(0);
}

Output

Using GCC 4.7.1 on Mac OS X 10.7.5, the output is:

h = 1.562500000000000e-02
-9: deriv(3.000967653589793e+00) = -9.901285883701071e-01
-8: deriv(3.016592653589793e+00) = -9.921976672293290e-01
-7: deriv(3.032217653589793e+00) = -9.940245152582091e-01
-6: deriv(3.047842653589793e+00) = -9.956086864580017e-01
-5: deriv(3.063467653589793e+00) = -9.969497940760287e-01
-4: deriv(3.079092653589793e+00) = -9.980475107000991e-01
-3: deriv(3.094717653589793e+00) = -9.989015683384429e-01
-2: deriv(3.110342653589793e+00) = -9.995117584851364e-01
-1: deriv(3.125967653589793e+00) = -9.998779321710066e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.157217653589793e+00) = -9.998779321710066e-01
 2: deriv(3.172842653589793e+00) = -9.995117584851364e-01
 3: deriv(3.188467653589793e+00) = -9.989015683384429e-01
 4: deriv(3.204092653589793e+00) = -9.980475107000991e-01
 5: deriv(3.219717653589793e+00) = -9.969497940760287e-01
 6: deriv(3.235342653589793e+00) = -9.956086864580017e-01
 7: deriv(3.250967653589793e+00) = -9.940245152582091e-01
 8: deriv(3.266592653589793e+00) = -9.921976672293291e-01
 9: deriv(3.282217653589793e+00) = -9.901285883701071e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -1.091570566584531e-01
 3: -9.361273104952932e-02
 4: -7.804555123490846e-02
 5: -6.245931771829009e-02
 6: -4.685783565504131e-02
 7: -3.124491392324735e-02
 8: -1.562436419383969e-02
 9: -1.776356839400250e-15
10: 1.562436419384087e-02
11: 3.124491392324558e-02
12: 4.685783565504250e-02
13: 6.245931771828653e-02
14: 7.804555123490846e-02
15: 9.361273104953050e-02
16: 1.091570566584471e-01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.577900251527073e+00
 3: 1.660540592568783e+00
 4: 9.969497901146779e-01
 5: 9.980475067341610e-01
 6: 9.989015643694564e-01
 7: 9.995117545137316e-01
 8: 9.998779281977730e-01
 9: 9.999999960264082e-01
10: 9.998779281978486e-01
11: 9.995117545137319e-01
12: 9.989015643693869e-01
13: 9.980475067340947e-01
14: 9.969497901149178e-01
15: 1.660540592568512e+00
16: -3.577900251527123e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 6.553266640232313e+01
 3: 1.898706817407991e+02
 4: -5.267598134705870e+01
 5: 3.608762837830827e+00
 6: 4.685783548517186e-02
 7: 3.124491378285654e-02
 8: 1.562436412277357e-02
 9: 3.226456139297321e-12
10: -1.562436412279785e-02
11: -3.124491378869069e-02
12: -4.685783548888563e-02
13: -3.608762837816180e+00
14: 5.267598134704988e+01
15: -1.898706817408119e+02
16: -6.553266640231030e+01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 8.382087654791743e+03
 3: -5.062815705775390e+03
 4: -7.597917560836845e+03
 5: 3.261984801532626e+03
 6: -4.336626618856812e+02
 7: 1.791410702361821e+01
 8: -9.998779233550453e-01
 9: -9.999999914294619e-01
10: -9.998779238596159e-01
11: 1.791410702341709e+01
12: -4.336626618847606e+02
13: 3.261984801532444e+03
14: -7.597917560838102e+03
15: -5.062815705774387e+03
16: 8.382087654792240e+03
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -9.999999914294619e-01 (eps = 0.015625)
h = 1.000000000000000e-05
-9: deriv(3.141502653589793e+00) = -9.999999959500000e-01
-8: deriv(3.141512653589793e+00) = -9.999999968000000e-01
-7: deriv(3.141522653589793e+00) = -9.999999975500000e-01
-6: deriv(3.141532653589793e+00) = -9.999999982000000e-01
-5: deriv(3.141542653589793e+00) = -9.999999987500000e-01
-4: deriv(3.141552653589793e+00) = -9.999999992000000e-01
-3: deriv(3.141562653589793e+00) = -9.999999995500000e-01
-2: deriv(3.141572653589793e+00) = -9.999999998000000e-01
-1: deriv(3.141582653589793e+00) = -9.999999999500000e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.141602653589793e+00) = -9.999999999500000e-01
 2: deriv(3.141612653589793e+00) = -9.999999998000000e-01
 3: deriv(3.141622653589793e+00) = -9.999999995500000e-01
 4: deriv(3.141632653589793e+00) = -9.999999992000000e-01
 5: deriv(3.141642653589793e+00) = -9.999999987500000e-01
 6: deriv(3.141652653589793e+00) = -9.999999982000000e-01
 7: deriv(3.141662653589793e+00) = -9.999999975500000e-01
 8: deriv(3.141672653589793e+00) = -9.999999968000000e-01
 9: deriv(3.141682653589793e+00) = -9.999999959500000e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -7.000000116589669e-05
 3: -5.999999941330713e-05
 4: -5.000000598739025e-05
 5: -3.999999683331386e-05
 6: -2.999999600591015e-05
 7: -2.000000257999327e-05
 8: -1.000000082740371e-05
 9: 0.000000000000000e+00
10: 1.000000082740371e-05
11: 2.000000165480742e-05
12: 2.999999508072429e-05
13: 3.999999590812801e-05
14: 5.000000413701854e-05
15: 5.999999663774957e-05
16: 6.999999746515327e-05
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.583333244325556e+00
 3: 1.666666318844711e+00
 4: 1.000000128999663e+00
 5: 1.000000691821058e+00
 6: 9.999995738881511e-01
 7: 9.999997049561470e-01
 8: 1.000000198388602e+00
 9: 1.000000075030488e+00
10: 1.000000144419428e+00
11: 9.999996509869722e-01
12: 9.999995893079155e-01
13: 1.000000645561765e+00
14: 1.000000028771196e+00
15: 1.666666187776715e+00
16: -3.583333074708150e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 1.027777535146502e+05
 3: 2.972222191231725e+05
 4: -8.263881528669108e+04
 5: 5.555518108303881e+03
 6: -6.636923522614542e-02
 7: 4.677328483600658e-02
 8: 1.991719546327412e-02
 9: -3.148201866790915e-03
10: -2.319389535987426e-02
11: -4.176186143567406e-02
12: 6.726872146719150e-02
13: -5.555525175695851e+03
14: 8.263880834779718e+04
15: -2.972222015189417e+05
16: -1.027777456120210e+05
17: 0.000000000000000e+00
18: 0.000000000000000e+00

Signs of madness...

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 2.050347140226725e+10
 3: -1.240740057099195e+10
 4: -1.858796490195886e+10
 5: 7.986101364079453e+09
 6: -1.059021715700324e+09
 7: 4.630176289959386e+07
 8: -3.687893612405425e+03
 9: -2.136279835945886e+03
10: -2.968840021291520e+03
11: 4.630204773690500e+07
12: -1.059022490436399e+09
13: 7.986100736580715e+09
14: -1.858796331554353e+10
15: -1.240739964045201e+10
16: 2.050347017082775e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -2.136279835945886e+03 (eps = 0.000010)

Notice how the results are going beserk with values in the billions in the last iteration. You have at the least a numerical stability problem, or perhaps you need to review the derivative formula. Note that even the run with the larger epsilon has a tendency towards large values on the later iterations.

---

Using 'official' derivative() function

Plugging the derivative function now present in the question into the code above yields even less stable answers in the 4th iteration:

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -6.248925477935563e+09
 3: 1.729549900845405e+11
 4: -2.600544559219368e+11
 5: -2.755286326619338e+11
 6: 8.100546069731433e+11
 7: -2.961111189495415e+09
 8: -7.936480686806423e+11
 9: 2.430177384467434e+11
10: 2.389084910067162e+11
11: -6.461168564124718e+10
12: -4.574822745530297e+10
13: -8.923883451146609e+10
14: 7.042030613792160e+10
15: 4.988386306820556e+10
16: -1.793262395787471e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = 2.430177384467434e+11 (eps = 0.000010)

I wonder to what extent the appearance of the zeros at array indexes 0, 1, 17, and 18 exacerbate the problem.

Answer 2

Numerical differentiation is a difficult problem. The key issue is that the finite difference approximation

f'(x) =(approx) (f(x+h)-f(x-h)) / (2*h)

is a recipe for cancellation.

This means you have to choose a careful tradeoff between two errors: If h is too large, you will incur a large approximation error in the finite difference. If h is too small, you will incur a large (possibly catastrophic) numerical error. This is illustrated quite nicely in the Wikipedia article on numerical differentiation. The problems are exacerbated as the order of the derivative increases.

This is why automatic differentiation is such a big deal -- in essence, it allows you to compute a precise derivative when all that is given is an algorithm that computes the base function. If the function is simple enough that you can determine the derivative symbolically, then you should probably do that.

If you do need to go with numerical differentiation, there is a lot of numerical finesse you can apply. Numerical Recipes has a good treatment -- take a look at Section 5.7, which begins on page 186.

Answer 3

This is a tricky subject. Numerical differentiation is ill-conditioned, and give rise to cancellation and rounding errors. The problem is significantly magnified in your case, where you do repeated differentiation.

You will get much better high order derivative estimations by determining the coefficients for your specific differentiation order N using m points, instead of doing repeated differentiation.

For example, just as you seem to be using the coefficients for the first order derivative:

{1, -8, 8, -1} / (12*h)

the coefficients for approximating a second order derivative using 5 points are:

{-1, 16, -30, 16, -1} / (12*h²)

You can solve for the coefficients for a general problem by making a Taylor series expansion around the sample points, finding the linear combination of these expansions that give the wanted derivative.

for the m points {x + n[1]*h, x+n[2]*h, x+n[3]*h, …, x+n[m]*h}, the coefficients (k) estimating the Nth derivative can be calculated by the following system of equations:

M*k=b

Where M is a m*m matrix:

M[i,j] = n[j]^(i-1), i,j = 1..m

and

b = [0, 0, 0, …, N!/h^N, …, 0, 0, 0],

where N! denotes the factorial of N.