Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define NDEPS (value)
- typedef double(*math_func)(double x);
- double nthDerivative(math_func f, double x, int N){
- //perfectly fine derivative function
- double derivative(math_func f, double x);
- 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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement