'********************************************************* '* Calculate a limited development of a real function * '* f(x) at point x0 with step h (to 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: Ln(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, Armand Colin * '* Editeur, Paris, 1991" [BIBLI 03]. * '* ----------------------------------------------------- * '* Subroutines used: * '* * '* 500: Sample function f(x) * '* 1000: Nth derivative of a discrete function * '* (N=1 to 5). * '* * '* Basic 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. '-----------------------------------------------------------} 'Program Ltddev DEFDBL A-H, O-Z DEFINT I-N DIM itable(5, 7) DATA 3,45,-9,1,0,0,60,3,270,-27,2,0,0,180,4,-488,338,-72,7,0,240 DATA 4,-1952,676,-96,7,0,240,5,1938,-1872,783,-152,13,288 'initialize table of coefficients FOR i = 1 TO 5 FOR j = 1 TO 7 READ itable(i, j) NEXT j NEXT i CLS PRINT PRINT " Limited development of a real function f(x):" PRINT PRINT " Function to develop: 1/(1-x)" PRINT INPUT " At point x0= ", x0 INPUT " Value of step h= ", h PRINT s = 1# FOR n = 1 TO 5 s = s * n GOSUB 1000 'r = arndederiv(x0,n,h) PRINT " a"; n; " = "; PRINT USING "##.##########"; R / s NEXT n PRINT END 'of main program 'sample function 1 f(xx)=1/(1-xx) 500 IF ABS(xx - 1) > 1E-12 THEN f = 1# / (1# - xx) ELSE f = 1E+12 END IF RETURN 'sample function 2 f(xx)=Log(xx) '500 if xx >1e-16 then ' f = LOG(x) ' else ' f = -1e16 ' END IF 'return } '-------------------------------------------------- 'returns the value of nth derivative of function 'f(x) at point x0 with increment step h (n=1 to 5). 'Note that for polynomial functions h must not be 'too small. Here h=1 gives best results than h=0.01 '-------------------------------------------------- 1000 'real*8 Function arnderiv(x0,n,h) IF (n > 5) THEN R = 0# RETURN END IF xx = x0: GOSUB 500: fa = f t = 0# FOR i = 1 TO itable(n, 1) xx = x0 - i * h: GOSUB 500: fb = f xx = x0 + i * h: GOSUB 500: fc = f IF (n MOD 2) <> 0 THEN t = t + itable(n, i + 1) * (fc - fb) ELSE t = t + itable(n, i + 1) * ((fc - fa) + (fb - fa)) END IF NEXT i t = t / itable(n, 7) FOR i = 1 TO n t = t / h NEXT i R = t RETURN 'end of file ltddev.bas