Path: utzoo!attcan!utgpu!watmath!maytag!watstat!dmurdoch
From: dmurdoch@watstat.waterloo.edu (Duncan Murdoch)
Newsgroups: comp.lang.pascal
Subject: Re: Mathematical Functions for Integration
Message-ID: <210@maytag.waterloo.edu>
Date: 24 May 89 13:50:41 GMT
References: <19715@adm.BRL.MIL>
Sender: daemon@maytag.waterloo.edu
Reply-To: dmurdoch@watstat.waterloo.edu (Duncan Murdoch)
Organization: U. of Waterloo, Ontario
Lines: 264

In article <19715@adm.BRL.MIL> PHSTUD8%DKNKURZ1.BITNET@cunyvm.cuny.edu writes:
>I'm searching desperatly for a function which integrates
>(any) function from zero or a point larger than that to infinity.

Numerical Recipes, by Press et al, suggests using the identity
   b             1/a  1    1
  I f(x) dx =   I     -  f(-) dt
   a             1/b  t^2  t
which works provided a and b have the same sign, to turn an upper limit of
infinity into a lower limit of zero.  If a=0, you'll have to break the integral
into two parts, i.e. 0 to c, c to infinity, and sum the results.

To integrate over a finite range, I've had very good luck with the routine
QUANC8, from Forsythe, Malcolm and Moler, Computer Methods for Mathematical
Computations.  Here's a copy, translated to Turbo Pascal v 5.

Duncan Murdoch

---- cut here -----   FMM.PAS  -------------
unit fmm;

interface

type
  float=double;
  integrand = function(var x;info:pointer):float;

procedure quanc8(fun:integrand;info:pointer;a,b,abserr,relerr:float;
                 var result,errest:float;
                 var nofun:integer;var flag:float);

{
    Subroutine from Forsythe, Malcolm and Moler (1977) Computer Methods
    for Mathematical Computations Chapter 5 to estimate the integral
    of fun(x) from a to b to a user provided tolerance. (Translated
    from Fortran to Ratfor and from Ratfor to Turbo Pascal by djm.)
    An automatic adaptive routine based on the 8-panel Newton-Cotes
    rule is used.

  Input:
    fun:      The name of the integrand function subprogram fun(x)
    info:     A pointer which will be passed to fun
    a:        The lower limit of integration.
    b:         "  upper   "   "      "      .  (b may be less than a)
    relerr:   A relative error tolerance (non-negative)
    abserr:   An absolute  "       "           "

  Output:
    result:   An approximation to the integral hopefully satisfying
              the least stringent of the two error tolerances
    errest:   An estimate of the magnitude of the actual error
    nofun:    The number of function evaluations
    flag:     A reliability indicator.  If flag is zero, then result
              probably satisfies the error tolerance.  If flag is
              XXX.YYY, then XXX = the number of intervals which have
              not converged and 0.YYY = the fraction of the interval
              (a,b) left to do when the limit on nofun was approached
}

implementation

function max(x,y:float):float;
begin
  if x>=y then
    max := x
  else
    max := y;
end;

procedure quanc8;
var
  w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp,
               qprev,qnow,qdiff,qleft,esterr,tolerr:float;
  qright: array [1..31] of float;
  f,x:    array[1..16] of float;
  fsave,xsave:  array[1..8,1..30] of float;
  levmin,levmax,levout,nomax,nofin,lev,i,j:integer;
  nim:longint;

label
  break;

begin
{    Stage 1        General Initialization }

{    Set constants   }

levmin := 1;
levmax := 30;
levout := 6;
nomax := 5000;
nofin := nomax - 8*(levmax - levout + 1 shl (levout+1));
{    Trouble when nofun reaches nofin   }

w0 :=   3956.0E0/14175.0E0;
w1 :=  23552.0E0/14175.0E0;
w2 :=  -3712.0E0/14175.0E0;
w3 :=  41984.0E0/14175.0E0;
w4 := -18160.0E0/14175.0E0;

{    Initialize running sums to zero.   }

flag   := 0.0E0;
result := 0.0E0;
cor11  := 0.0E0;
errest := 0.0E0;
area   := 0.0E0;
nofun  := 0;

if (a = b) then
     exit;

{    Stage 2   Initialization for first interval   }

lev := 0;
nim := 1;
x0 := a;
x[16] := b;
qprev := 0.0;
f0 := fun(x0,info);
stone := (b-a)/16.0E0;
x[ 8] := (  x0  + x[16])/2.0E0;
x[ 4] := (  x0  + x[ 8])/2.0E0;
x[12] := (x[ 8] + x[16])/2.0E0;
x[ 2] := (  x0  + x[ 4])/2.0E0;
x[ 6] := (x[ 4] + x[ 8])/2.0E0;
x[10] := (x[ 8] + x[12])/2.0E0;
x[14] := (x[12] + x[16])/2.0E0;

for j:=1 to 8 do
     f[2*j] := fun(x[2*j],info);
nofun := 9;


repeat
{    Stage 3   Central Calculation   }
     {         Requires qprev,x0,x2,...,x16,f0,f2,...,f16.   }

     {         Calculates x1,x3,...,x15, f1,f3,...,f15,qleft,qright,qnow   }
     {              qdiff,area.   }
     x[1] := (x0 + x[2])/2.0E0;
     f[1] := fun(x[1],info);
     for j:=2 to 8 do    {3 to 15 by 2}
     begin
          x[2*j-1] := (x[2*j-2] + x[2*j])/2.0E0;
          f[2*j-1] := fun(x[2*j-1],info);
     end;
     nofun := nofun + 8;
     step := (x[16] - x0)/16.0E0;

     qleft := (w0*(f0+f[8]) + w1*(f[1]+f[7]) + w2*(f[2]+f[6])
            + w3*(f[3]+f[5]) + w4*f[4])*step;
     qright[lev+1] := (w0*(f[8]+f[16]) + w1*(f[9]+f[15])
            + w2*(f[10]+f[14]) + w3*(f[11]+f[13]) + w4*f[12])*step;
     qnow := qleft + qright[lev+1];

     qdiff := qnow - qprev;
     area := area + qdiff;

     {    Stage 4   Interval Convergence test   }

     esterr := abs(qdiff)/1023.0E0;
     tolerr := max(abserr,relerr*abs(area))*(step/stone);

     if (lev < levmin) or
         ((lev < levmax) and (nofun <= nofin) and (esterr > tolerr)) then
     begin
          {    Stage 5   No convergence   }
          {    Locate next interval.   }

          nim := 2*nim;
          lev := lev + 1;

          {    Store right hand elements for future use   }

          for i:=1 to 8 do
          begin
               fsave[i,lev] := f[i+8];
               xsave[i,lev] := x[i+8];
          end;
          {    Assemble left hand elements for immediate use   }

          qprev := qleft;
          for  i:=1 to 8 do
          begin
               j := -i;
               f[2*j+18] := f[j+9];
               x[2*j+18] := x[j+9];
          end;
          {    Go on to next iteration   }
     end
     else
     begin
          if (nofun > nofin) then
          begin
               {    Stage 6   Trouble Section   }
               {    Number of function evaluations is about to exceed   }
               {         limit   }

               nofin := 2*nofin;
               levmax := levout;
               flag := flag + (b-x0)/(b-a);
          end
          else if (lev >= levmax) then
          begin
               {    Current level is levmax   }

               flag := flag + 1.0;
          end;

          {    Stage 7   Interval converged   }
          {         Add contributions into running sums   }

          result := result + qnow;
          errest := errest + esterr;
          cor11 := cor11 + qdiff/1023.0E0;

          {    Locate next interval   }

          while (nim <> 2*(nim div 2)) do
          begin
               nim := nim div 2;
               lev := lev - 1;
          end;
          nim := nim + 1;

          if (lev <= 0) then
               goto break;          { Completely done   }

          {    Assemble elements required for the next interval   }

          qprev := qright[lev];
          x0 := x[16];
          f0 := f[16];
          for i:=1 to 8 do
          begin
               f[2*i] := fsave[i,lev];
               x[2*i] := xsave[i,lev];
          end;
     end;
until false;

{    Stage 8   Finalize and return   }

break:

result := result + cor11;

{    Make sure errest not less than roundoff level   }

if (errest = 0.0E0) then
     exit;

temp := abs(result) + errest;
while (temp = abs(result)) do
begin
     errest := 2.0E0*errest;
     temp := abs(result) + errest;
end;

end;

end.
--------- end of FMM.PAS -----------