!************************************************************************ !* 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.891731795301894E-011 * !* res = -69.4444444444433 * !* nend = 3 * !* hend = 6.250000000000000E-004 * !* * !* -------------------------------------------------------------------- * !* "Numerical Algorithms with C, By Gisela Engeln-Muellges * !* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * !* * !* F90 Release 1.0 By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !************************************************************************ Program Test_Difrom integer error, n, nend real*8 x,prec,h, res,schaetz, hend print *,' ' print *,' Numerical differentation according to Romberg' print *,' ' print *,' Test function f(x) = 1/x' print *,' ' write(*,10,advance='no'); read *, x write(*,20,advance='no'); read *, prec write(*,30,advance='no'); read *, n write(*,40,advance='no'); read *, h call difrom(x,prec,n,h,res,schaetz,nend,hend,error) call aus(error,schaetz,res,nend,hend) stop 10 format(' Put in x-value at which you want to evaluate derivative: ') 20 format(' Put in desired accuracy: ') 30 format(' Maximal number of columns in Romberg scheme: ') 40 format(' Starting step size: ') END !print results Subroutine aus(error, schaetz, res, nend, hend) integer error, nend real*8 schaetz,res,hend if (error.ne.1) then print *,' ' print *,' Results:' print *,' ' print *,' er_app= ', schaetz print *,' res = ', res print *,' nend = ', nend print *,' hend = ', hend print *,' ' end if return End !test function to derivate real*8 Function func(x) real*8 x func = 1.d0 / x return end Subroutine difrom(x0, eps, n, h, res, er_app, nend, hend, error) real*8 x0,eps,h,res,er_app,hend integer 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 error): * !* * !* 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) * !* * !************************************************************************* parameter(NMAX=10, MACH_EPS = 1.2d-16) real*8 h2, d1, d2, func real*8 d(0:NMAX) if (n <= 1.or.eps <= 0.0.or.h < MACH_EPS) then !check input error=1 return end if h2 = 2.d0 * h d(0) = (func(x0 + h) - func(x0 - h)) / 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. * !************************************************************************ error=2 do j = 1, n - 1 d(j) = 0.d0 d1 = d(0) h2 = h h = h * 0.5d0 nend = j if (h < MACH_EPS) then error = 3 !step size less than machine constant goto 10 end if d(0) = (func(x0 + h) - func(x0 - h)) / h2 m=4 do i = 1, j d2 = d(i) d(i) = (m * d(i-1) - d1) / (m-1) d1 = d2 m = m * 4 end do er_app = DABS(d(j) - d(j-1)) if (er_app < eps) then error = 0 !desired accuracy reached goto 10 end if end do 10 res = d(nend) !save final values hend = h Return END !end of file tdifrom.f90