'**************************************************************
'* 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].                       *
'*                            Basic 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   :  .4596976941318507                   *
'*                                                            *
'* Obtained precision  :  9.599099293211566E-011              *
'*                                                            *
'* Number of iterations:  4                                   *
'*                                                            *
'**************************************************************
defint i-n
defdbl a-h,o-z

'       x0      : double;     {begin x value}
'       x1      : double;     {end   x value}
'       prec    : double;     {desired precision}
'      xintegral: double;     {result of integral}
'       obtprec : double;     {obtained precision}

'       niter   : integer;    {number of actual iterations}

  cls
  print
  print " Integral of a function F(X) by Romberg''s method"
  print
  print " Input begin and end values for x variable:"
  print
  input "         X0 = ", a
  input "         X1 = ", b
  print
  input " Input desired precision: ", prec

  itermin=1 : itermax=50

  'call Romberg subroutine
  gosub 2000

  print
  print
  print " Value of integral   : "; xintegral
  print
  print " Obtained precision  : "; obtprec
  print
  print " Number of iterations: "; niter

END 'of main program


1000 'Given function to integrate
    FUNC = SIN(x)
return


  '*******************************************************
  '* 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           *
  '*  xintegral      the integral of FUNC(X) from a to b *
  '*                                                     *
  '*******************************************************
2000 'Function RombergIntegral
  MAXITER = 15
  dim t(MAXITER,MAXITER)
    if itermax>MAXITER then itermax=MAXITER
    x=a : gosub 1000 : r=FUNC
    x=b : gosub 1000
    ta = (r + FUNC) / 2#
    niter=0
    pas=b-a
    t(0,0)=ta*pas
2100 niter=niter+1
    pas=pas/2#
    s=ta
    for i=1 to 2^niter-1
      x=a+pas*i : gosub 1000
      s = s + FUNC
    next i
    t(0,niter)=s*pas
    r=1#
    for i=1 to niter
      r=r*4#
      j=niter-i
      t(i,j)=(r*t(i-1,j+1) - t(i-1,j))/(r-1#)
    next i
    obtprec = ABS(t(niter,0) - t(niter-1,0))
    if niter > itermax then goto 2200
    if niter<itermin then goto 2100
    if obtprec>prec then goto 2100
2200 xintegral = t(niter,0)
return

'End of file tromberg.bas