'************************************************************************ '* Test program: Numerical differentation according to Romberg * '* -------------------------------------------------------------------- * '* SAMPLE RUN: * '* (Find the first derivative of function f(x)=1/x, for x=0.12). * '* * '* Numerical differentation according to Romberg * '* * '* Test function f(x) = 1/x * '* * '* Put in x-value at which you want to evaluate derivative: 0.12 * '* Put in desired accuracy: 1e-10 * '* Maximal number of columns in Romberg scheme: 4 * '* Starting step size: 0.005 * '* * '* Results: * '* * '* er_app= 8.893152880773414D-011 * '* res = -69.44444444444329 * '* nend = 3 * '* hend = .00062 * '* * '* -------------------------------------------------------------------- * '* "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************ 'Program Test_Difrom DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 CLS PRINT PRINT " Numerical differentation according to Romberg" PRINT PRINT " Test function f(x) = 1/x" PRINT INPUT " Put in x-value at which you want to evaluate derivative: ", x0 INPUT " Put in desired accuracy: ", eps INPUT " Maximal number of columns in Romberg scheme: ", n INPUT " Starting step size: ", h GOSUB 1000 'call difrom(x,eps,n,h,res,erapp,nend,hend,ierror) GOSUB 2000 'call aus(ierror,erapp,res,nend,hend) END 'test function to derivate 500 'Function func(xx) func = 1# / xx RETURN 1000 'Subroutine difrom(x0, eps, n, h, res, er_app, nend, hend, error) '************************************************************************* '* Computes an approximation for the first derivative of func at x0 * '* using the ROMBERG method. * '* --------------------------------------------------------------------- * '* Input parameters : * '* * '* real*8 func(real*8) (External) Name of function to be * '* differentiated. * '* real*8 x0 value of abscissa at which derivative is to * '* be found. * '* real*8 eps desired accuracy. * '* integer n max. number of columns in the Romberg scheme * '* (n > 1). * '* real*8 h initial step size. * '* * '* Output parameters: * '* * '* real*8 res approximate derivative * '* real*8 er_app error estimate for res * '* integer nend number of columns actually used in scheme * '* real*8 hend final step size * '* * '* Return value (integer ierror): * '* * '* 0: no error: er_app < eps * '* 1: n < 1 or eps <= 0 or h < MACH_EPS * '* 2: desired accuracy not reached after n steps * '* 3: step size drooped below MACH_EPS * '* 4: lack of sufficient memory (not used here) * '* * '************************************************************************* NMAX = 10: XMACHEPS = 1.2D-16 DIM d(NMAX) IF n <= 1 OR eps <= 0# OR h < XMACHEPS THEN ierror = 1 RETURN END IF h2 = 2# * h xx = x0 + h: GOSUB 500: func1 = func xx = x0 - h: GOSUB 500 d(0) = (func1 - func) / h2 '************************************************************************ '* This loop runs until the maximum of Romberg rows is filled or until * '* the step size use drops below the machine constant or if the desired * '* accuracy is reached. * '************************************************************************ ierror = 2 FOR j = 1 TO n - 1 d(j) = 0# d1 = d(0) h2 = h h = h * .5# nend = j IF h < XMACHEPS THEN 'step size less than machine constant ierror = 3 GOTO 10 END IF xx = x0 + h: GOSUB 500: func1 = func xx = x0 - h: GOSUB 500 d(0) = (func1 - func) / h2 m = 4 FOR i = 1 TO j d2 = d(i) d(i) = (m * d(i - 1) - d1) / (m - 1) d1 = d2 m = m * 4 NEXT i erapp = ABS(d(j) - d(j - 1)) IF erapp < eps THEN ierror = 0 'desired accuracy reached GOTO 10 END IF NEXT j 10 res = d(nend) 'save final values hend = h RETURN 'print results 2000 'Subroutine aus(error, schaetz, res, nend, hend) IF ierror <> 1 THEN PRINT PRINT " Results:" PRINT PRINT " er_app= "; erapp PRINT " res = "; res PRINT " nend = "; nend PRINT " hend = "; hend PRINT END IF RETURN 'end of file tdifrom.bas