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").
Why does this happen? Is it related to the nature of the sin() function? Or is it because of floating-point accuracy issues?
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).
#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);
}
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.
derivative() functionPlugging 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.