/******************************************************** * Estimate the nth derivative of a real function f(x) * * at point x with step h (n from 1 to 5). * * ----------------------------------------------------- * * SAMPLE RUN: * * * * Nth derivative of a real function f(x): * * * * Function to derivate: x^3 * * * * At point x0= 10 * * Order of derivative (1 to 5): 1 * * Value of step h= 0.1 * * * * f'(x0)=300.000000 * * * * At point x0= 10 * * Order of derivative (1 to 5): 2 * * Value of step h= 0.1 * * * * f''(x0)=60.000000 * * * * ----------------------------------------------------- * * Ref.: "Analyse numérique en C By Alain Reverchon and * * Marc Ducamp, Armand Colin Editeur, Paris, * * 1990" [BIBLI 14]. * * * * C++ (simplified) version by J-P Moreau. * * (www.jpmoreau.fr) * ********************************************************/ #include #include //sample function f(x)=x^3 double f(double x) { return x*x*x; } //returns the value of nth derivative of function //f(x) at point x0 with increment step h //Note that for polynomial functions h must not be //too small. Here h=1 gives best results than h=0.01 double ar_nderiv(double x,int n,double h) { static int table[5][7] = { {3, 45, -9, 1, 0, 0, 60}, {3, 270, -27, 2, 0, 0, 180}, {4, -488, 338, -72, 7, 0, 240}, {4, -1952, 676, -96, 7, 0, 240}, {5, 1938, -1872, 783, -152, 13, 288} }; int i, *coeff; double fa,fb,fc,t; if(n>5) return 0; coeff = table[n-1]; fa = f(x); for (i=1, t=0; i<=coeff[0]; i++) { fb = f(x-i*h); fc = f(x+i*h); t = t + coeff[i]*((n%2) ? (fc-fb) : ((fc-fa) + (fb-fa))); } t /= coeff[6]; for (i=1; i<=n; i++) t /= h; return t; } void main() { double d,h,x; int n; printf("\n Nth derivative of a real function f(x):\n\n"); printf(" Function to derivate: x^3\n\n"); printf(" At point x0= "); scanf("%lf",&x); printf(" Order of derivative (1 to 5): "); scanf("%d",&n); printf(" Value of step h= "); scanf("%lf",&h); printf("\n"); d = ar_nderiv(x,n,h); switch(n) { case 1: printf(" f'(x0)=%f\n\n",d); break; case 2: printf(" f''(x0)=%f\n\n",d); break; case 3: printf(" f'''(x0)=%f\n\n",d); break; case 4: printf(" f''''(x0)=%f\n\n",d); break; case 5: printf(" f'''''(x0)=%f\n\n",d); break; default: printf(" n must be between 1 and 5 !\n\n"); } } //end of file nderiv.cpp