'********************************************************* '* 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 numerique en C By Alain Reverchon and * '* Marc Ducamp, Armand Colin Editeur, Paris, * '* 1990" [BIBLI 14]. * '* * '* Basic Version By J-P Moreau. * '* (www.jpmoreau.fr) * '********************************************************* 'Program Nderiv DEFDBL A-H, O-Z DEFINT I-N 'Main variables: 'double d: nth derivative at x ' h: increment step ' x: abscissa of function f(x) 'integer n: order of derivative (1 to 5) CLS PRINT PRINT " Nth derivative of a real function f(x):" PRINT PRINT " Function to derivate: x^3" PRINT INPUT " At point x0= ", x INPUT " Order of derivative (1 to 5): ", n INPUT " Value of step h= ", h PRINT GOSUB 1000 'd = nderiv(x,n,h) IF n = 1 THEN PRINT USING " f'(x0) = ####.######"; d ELSEIF n = 2 THEN PRINT USING " f''(x0) = ####.######"; d ELSEIF n = 3 THEN PRINT USING " f'''(x0) = ####.######"; d ELSEIF n = 4 THEN PRINT USING " f''''(x0) = ####.######"; d ELSEIF n = 5 THEN PRINT USING " f'''''(x0) = ####.######"; d END IF IF n < 1 OR n > 5 THEN PRINT " n must be between 1 and 5 !" END IF PRINT END 'sample function f(xx)=xx^3 500 f = xx * xx * xx RETURN '--------------------------------------------------- ' 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 same results than h=0.01 '--------------------------------------------------- 1000 'Function nderiv(x, n, h) DIM itable(5, 7) itable(1, 1) = 3: itable(1, 2) = 45: itable(1, 3) = -9: itable(1, 4) = 1 itable(1, 5) = 0: itable(1, 6) = 0: itable(1, 7) = 60 itable(2, 1) = 3: itable(2, 2) = 270: itable(2, 3) = -27: itable(2, 4) = 2 itable(2, 5) = 0: itable(2, 6) = 0: itable(2, 7) = 180 itable(3, 1) = 4: itable(3, 2) = -488: itable(3, 3) = 338: itable(3, 4) = -72 itable(3, 5) = 7: itable(3, 6) = 0: itable(3, 7) = 240 itable(4, 1) = 4: itable(4, 2) = -1952: itable(4, 3) = 676: itable(4, 4) = -96 itable(4, 5) = 7: itable(4, 6) = 0: itable(4, 7) = 240 itable(5, 1) = 5: itable(5, 2) = 1938: itable(5, 3) = -1872: itable(5, 4) = 783 itable(5, 5) = -152: itable(5, 6) = 13: itable(5, 7) = 288 IF n > 5 THEN d = 0# RETURN END IF xx = x: GOSUB 500: fa = f 'fa=f(x) d = 0# FOR i = 1 TO itable(n, 1) xx = x - i * h: GOSUB 500: fb = f 'fb=f(x-i*h) xx = x + i * h: GOSUB 500: fc = f 'fc=f(x+i*h) IF (n MOD 2) <> 0 THEN d = d + itable(n, i + 1) * (fc - fb) ELSE d = d + itable(n, i + 1) * ((fc - fa) + (fb - fa)) END IF NEXT i d = d / itable(n, 7) FOR i = 1 TO n d = d / h NEXT i RETURN 'end of file nderiv.bas