/******************************************************** * Calculate a limited development of a real function * * f(x) at point x0 with step h (at order 5). * * ----------------------------------------------------- * * SAMPLE RUN: * * * * Limited development of f(x) at x0: * * * * Function to develop: 1/(1-x) * * * * At point x0= 0 * * Value of step h= 0.1 * * * * a1 = 1.0000416250 * * a2 = 1.0000416250 * * a3 = 1.0011207928 * * a4 = 1.0011207928 * * a5 = 1.0136468470 * * * * (The correct answer is a1=a2=a3=a4=a5=1) * * * * Function to develop: log(x) * * * * At point x0= 1 * * Value of step h= 0.1 * * * * a1 = 1.0000057552 * * a2 = -0.5000050520 * * a3 = 0.3334508142 * * a4 = -0.2501062334 * * a5 = 0.2011301593 * * * * (the correct answer is a1=1, a2=-0.5, a3=1/3, * * a4=-1/4 and a5=1/5) * * ----------------------------------------------------- * * Ref.: "Mathematiques en Turbo Pascal By Alain * * Reverchon and Marc Ducamp, Editions Eyrolles, * * Paris, 1991" [BIBLI 03]. * * * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * ********************************************************* Note: as you can see this direct method is not suitable to calculate a limited development in an accurate way. The next programs (ltddev1 to ltddev3) will allow us to calculate limited developments at x=0 up to a high order with a good accuracy for functions f*g, f/g or gof, knowing the limited developments of f(x) and g(x) at x=0. ----------------------------------------------------------*/ #include #include double d,h,r,s,x0; int n; //sample function f(x)=1/(1-x) double f(double x) { if (fabs(x-1) >1e-12) return 1.0/(1.0-x); else return 1e12; } /*sample function f(x)=Ln(x) double f(double x) { if (x >1e-12) return log(x); else return -1e12; } */ //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() { printf("\n Limited development of a real function f(x):\n\n"); printf(" Function to develop: 1/(1-x)\n"); printf("\n At point x0= "); scanf("%lf", &x0); printf("\n Value of step h= "); scanf("%lf", &h); printf("\n\n"); s=1; for (n=1; n<6; n++) { s *= n; r = ar_nderiv(x0,n,h) / s; printf(" a%d = %13.10f\n", n, r); } printf("\n"); } // end of file ltddev.cpp