WEEHAAVIE|
begin
comment library A0, A6, A12;

   comment  HAVIE INTEGRATOR. ALGORITHM 257, Robert N. Kubik, CACM 8 (1965) 381;

real  a,b,eps,mask,y,answer; integer m max, count;

real  procedure  havieintegrator(y,a,b,eps,integrand,m);
value  a,b,eps,m;
integer  m;
real  integrand,y,a,b,eps;
   comment  This algorithm performs numerical integration of definite
   integrals using an equidistant sampling of the function and repeated halving of
   the sampling interval. Each halving allows the calculation of a trapezium and a
   tangent formula on a finer grid, but also the calculation of several higher
   order formulas which are defined implicitly. The two families of approximate
   solutions will normally bracket the value of the integral and from these
   convergence is tested on each of the several orders of approximation. The
   algorithm is based on a private communication from F. Haavie of the Institutt
   for Atom-energi Kjeller Research Establishment, Norway. A Fortran version is in
   use on the Philco-2000. ...;
begin
real  h,endpts,sumt,sumu,d;
integer  i,j,k,n;
real  array  t,u,tprev,uprev[1:m];
   count := 0;
   y := a; endpts := integrand;
   y := b; endpts := 0.5 × (integrand + endpts);
   sumt := 0.0; i := n := 1; h := b - a;
estimate:
   count := count + 1;
   t[1] := h × (endpts + sumt); sumu := 0.0;
      comment  t[1] = h×(0.5×f[0]+f[1]+f[2]+...+0.5×f[2 ^ (i-1)]);
   y := a - h/2.0;
   for  j := 1 step  1 until  n do
      begin
      y := y + h; sumu := sumu + integrand;
      end ;
   u[1] := h × sumu; k := 1;
      comment  u[1] = h×(f[1/2]+f[3/2]+...f[(2 ^ i-1)/2],
      k corresponds to approximate solution with truncation
      error term of order 2k;
test:
   if  abs(t[k]-u[k])  <  eps then
      begin
      havieintegrator := 0.5 × (t[k] + u[k]);
      goto  exit
      end ;
   if  k  ±  i then
      begin
      d := 2  ^  (2×k);
      t[k+1] := (d × t[k] - tprev[k]) / (d - 1.0);
      tprev[k] := t[k];
      u[k+1] := (d × u[k] - uprev[k]) / (d - 1.0);
      uprev[k] := u[k];
         comment  This implicit formulation of the higher-order
            integration formulas is given in [ROMBERG, W. ...];
      k := k + 1;
      if  k = m then
         begin
          havieintegrator := 0.5 × (t[k] + u[k]);
         goto  exit
         end ;
      goto  test
      end ;
   h := h / 2.0; sumt := sumt + sumu;
   tprev[k] := t[k]; uprev[k] := u[k];
   i := i + 1; n := 2 × n;
   goto  estimate;
exit:
end  havieintegrator;


eps  := 1.0º-10;
writetext(30, [epsilon*=*]);
write(30, format([d.dddddddddddº+ndc_]), eps);

m max  := 20;
writetext(30, [m max*=*]);
write(30, format([nddc_]), m max);

a := 0; b := 1;
answer := havieintegrator(y,0,b,eps,1.0/(1.0+y×y),m max);
writetext(30, [count*=*]);
write(30, format([nddddddddc_]), count);
write(30, format([+d.dddddddddddº+nd_]), answer);
writetext(30, [*approximates*+7.8539816339745º*-1]);
writetext(30, [,*relative*error*]);
write(30, format([+d.dddddddddddº+ndc_]), (answer - 0.78539816339745)/0.78539816339745);

a := eps; b := 1;
answer := havieintegrator(y,eps,b,eps,(y ^ (-y)),m max);
writetext(30, [count*=*]);
write(30, format([nddddddddc_]), count);
write(30, format([+d.dddddddddddº+nd_]), answer);
writetext(30, [*approximates*+1.2912859970627º*+0]);
writetext(30, [,*relative*error*]);
write(30, format([+d.dddddddddddº+ndc_]), (answer - 1.2912859970627)/1.2912859970627);

a := 0; b := 1;
answer := havieintegrator(y,0,b,eps,ln(1.0+y)/(1.0 + y^2),m max);
writetext(30, [count*=*]);
write(30, format([nddddddddc_]), count);
write(30, format([+d.dddddddddddº+nd_]), answer);
writetext(30, [*approximates*+2.7219826128795º*-1]);
writetext(30, [,*relative*error*]);
write(30, format([+d.dddddddddddº+ndc_]), (answer - 0.27219826128795)/0.27219826128795);


a := 0; b := 1;
answer := havieintegrator(y,0,b,eps,(y+y)/(1.0+y×y),m max);
writetext(30, [count*=*]);
write(30, format([nddddddddc_]), count);
write(30, format([+d.dddddddddddº+nd_]), answer);
writetext(30, [*approximates*+6.9314718055995º*-1]);
writetext(30, [,*relative*error*]);
write(30, format([+d.dddddddddddº+ndc_]), (answer - 0.69314718055995)/0.69314718055995);


a:= 0.0; b:= 4.3;
answer := havieintegrator(y, a, b, eps, exp(-y×y), m max);
writetext(30, [count*=*]);
write(30, format([nddddddddc_]), count);
write(30, format([+d.dddddddddddº+nd_]), answer);
writetext(30, [*approximates*+8.8622692439507º*-1]);
writetext(30, [,*relative*error*]);
write(30, format([+d.dddddddddddº+ndc_]), (answer - 0.88622692439507)/0.88622692439507);

FINISH:

end
|
|