!************************************************************** !* Program to demonstrate integration of a real function F(x) * !* by Romberg's method * !* ---------------------------------------------------------- * !* REFERENCE: "Mathematiques en Turbo-Pascal (Part 1) By * !* Marc Ducamp and Alain Reverchon, Eyrolles, * !* Paris, 1987" [BIBLI 03]. * !* F90 version by J-P Moreau * !* (www.jpmoreau.fr) * !* ---------------------------------------------------------- * !* SAMPLE RUN: * !* (Integrate sin(x) from x=0 to x=1) * !* * !* Integral of a function F(X) by Romberg's method * !* * !* Input begin and end values for x variable: * !* * !* X0 = 0 * !* X1 = 1 * !* * !* Input desired precision: 1e-10 * !* * !* * !* Value of integral : 0.459697694131851 * !* * !* Obtained precision : 9.599082639866197E-011 * !* * !* Number of iterations: 4 * !* * !************************************************************** PROGRAM TEST_ROMBERG real*8 x0, & ! begin x value x1, & ! end x value prec, & ! desired precision integral, & ! result of integral obtprec ! obtained precision integer niter, & ! number of actual iterations itermin, & ! minimal number of iterations itermax ! maximal number of iterations print *,' ' print *,'Integral of a function F(X) by Romberg''s method' print *,' ' print *,'Input begin and end values for x variable:' print *,' ' write(*,"(' X0 = ')",advance='no') read *, x0 write(*,"(' X1 = ')",advance='no') read *, x1 print *,' ' write(*,"(' Input desired precision: ')",advance='no') read *, prec itermin=1; itermax=50 ! call Romberg subroutine integral = RombergIntegral(x0,x1,prec,obtprec,niter,itermin,itermax) print *,' ' print *,' ' print *,'Value of integral : ', integral print *,' ' print *,'Obtained precision : ', obtprec write(*,10) niter stop 10 format(/' Number of iterations: ',i3//) END !of main program ! Given function to integrate real*8 Function FUNC(x) real*8 x FUNC = dsin(x) end !******************************************************* !* Integral of a function FUNC(X) by Romberg's method * !* --------------------------------------------------- * !* INPUTS: * !* a begin value of x variable * !* b end value of x variable * !* prec desired precision * !* itermin minimum number of ietrations * !* itermax maximum number of iterations * !* * !* OUTPUTS: * !* obtprec obtained precision for integral * !* niter number of iterations done * !* integral the integral of FUNC(X) from a to b * !* * !******************************************************* real*8 Function RombergIntegral(a,b,prec,obtprec,niter,itermin,itermax) parameter(MAXITER = 15) real*8 a,b,prec,obtprec,FUNC real*8 t(MAXITER,MAXITER) real*8 pas,r,s,ta if (itermax>MAXITER) itermax=MAXITER r = FUNC(a) ta = (r + FUNC(b)) / 2.d0 niter=0 pas=b-a t(0,0)=ta*pas 100 niter=niter+1 pas=pas/2.d0 s=ta do i=1, 2**niter-1 s = s + FUNC(a+pas*i) end do t(0,niter)=s*pas r=1.d0 do i=1, niter r=r*4.d0 j=niter-i t(i,j)=(r*t(i-1,j+1) - t(i-1,j))/(r-1.d0) end do obtprec = DABS(t(niter,0) - t(niter-1,0)) if (niter > itermax) goto 200 if (niterprec) goto 100 200 RombergIntegral = t(niter,0) return end ! End of file tromberg.f90