begin

library A0, A1, A5, A12, A15;

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

real a, b, eps, y, answer, two pi, pi, half pi; integer m max;

real procedure haaviequadrature(x, a, b, eps, integrand, m);
value a, b, eps, m;
integer m;
real integrand, x, 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];
   x := a; endpts := integrand;
   x := b; endpts := 0.5 × (integrand + endpts);
   sumt := 0.0; i := n := 1; h := b - a;
ESTIMATE:
   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)]);
   x := 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 or k = m then
      begin
      haaviequadrature := 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;
      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 haaviequadrature;

comment Following is a driver program to test haaviequadrature;

writetext(30, {This$is$the$Haavie$quadrature$program.{cc}});

two pi := 6.283185307179588;
pi := two pi / 2;
half pi := pi / 2;
open(30);

m max  := 20;
writetext(30, {{c}mmax$=});
write(30, format({nddc}), m max);

eps  := 0.0000000001;
writetext(30, {{c}epsilon$=});
write(30, format({-d.dddddddddddc}), eps);

comment basic accuracy tests;

answer := haaviequadrature(y, 0, 1, eps, y, m max);
writetext(30, {{c}integral(0, 1, x)$=$$$$$$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1/2) × 2); newline(30, 1);

answer := haaviequadrature(y, 0, 1, eps, y^2, m max);
writetext(30, {{c}integral(0, 1, x$sqr)$=$$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1/3) × 3); newline(30, 1);

answer := haaviequadrature(y, 0, 1, eps, y^9, m max);
writetext(30, {{c}integral(0, 1, x$to$9)$=$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1/10) × 10); newline(30, 1);

answer := haaviequadrature(y, 0, 1, eps, exp(y), m max);
writetext(30, {{c}integral(0, 1, exp(x))$=$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1.718281828459)/1.718281828459); newline(30, 1);

answer := haaviequadrature(y, 1, 2.718281828459045, eps, ln(y), m max);
writetext(30, {{c}integral(1, e, ln(x))$=$$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1)); newline(30, 1);

answer := haaviequadrature(y, 0, 1, eps, arctan(y), m max);
writetext(30, {{c}integral(1, e, arctan(x))$=$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.438824573117)/0.438824573117); newline(30, 1);

comment Following tests are from Kruseman-Aretz;

answer := haaviequadrature(y, 0, half pi, eps, cos(y), m max);
writetext(30, {{c}integral(0, pi/2, cos(x))$=$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1)); newline(30, 1);

answer := haaviequadrature(y, 0, 4.3, eps, exp(-y×y), m max);
writetext(30, {{c}integral(0, 4.3, exp(-x$times$x))$=$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.886226924395)/0.886226924395); newline(30, 1);

answer := haaviequadrature(y, 1, 10, eps, ln(y), m max);
writetext(30, {{c}integral(1, 10, ln(x))$=$$$$$$$$$$$$$$$$$$});
write(30, format({-dd.dddddddddd}), answer);
writetext(30, {:$$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-14.02585092994)/14.02585092994); newline(30, 1);

answer :=  haaviequadrature(y, 0, 20, eps, sqrt(y)/(exp(y-4)+1), m max);
writetext(30, {{c}integral(0, 20, sqrt(x)/(exp(x-4)+1))$=$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-5.77072601204)/5.77072601204); newline(30, 1);

comment Following additional tests;

answer := haaviequadrature(y, 1, 2, eps, exp(y), m max);
writetext(30, {{c}integral(1, 2, exp(x))$=$$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-4.670774270472)/4.670774270472); newline(30, 1);

answer := haaviequadrature(y, 0, half pi / 2, eps, sin(y), m max);
writetext(30, {{c}integral(0, pi/4, sin(x))$=$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.29289321881)/0.29289321881); newline(30, 1);

answer := haaviequadrature(y, half pi / 2, half pi, eps, sin(y), 16);
writetext(30, {{c}integral(pi/4, pi/2, sin(x))$=$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.707106781186)/0.707106781186); newline(30, 1);

answer := haaviequadrature(y, 0, half pi, eps, sin(y), m max);
writetext(30, {{c}integral(0, pi/2, sin(x))$=$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1)); newline(30, 1);

answer := haaviequadrature(y, half pi, pi, eps, sin(y), m max);
writetext(30, {{c}integral(pi/2, pi, sin(x))$=$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1)); newline(30, 1);

answer := haaviequadrature(y, 0, pi, eps, sin(y), m max);
writetext(30, {{c}integral(0, pi, sin(x))$=$$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-2)/2); newline(30, 1);

answer := haaviequadrature(y, 0, two pi, eps, sin(y), m max);
writetext(30, {{c}integral(0, 2pi, sin(x))$=$$$$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0)); newline(30, 1);

answer := haaviequadrature(y,0,1,eps,1/(1+y×y),m max);
writetext(30, {{c}integral(0, 1, 1/(1+x$times$x))$=$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.78539816339)/0.78539816339); newline(30, 1);

answer := haaviequadrature(y,eps,1,eps,(y ^ (-y)),m max);
writetext(30, {{c}integral(eps, 1, x$to$(-x))$=$$$$$$$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-1.29128599706)/1.29128599706); newline(30, 1);

answer := haaviequadrature(y,0,1,eps,ln(1+y)/(1 + y^2),m max);
writetext(30, {{c}integral(0, 1, ln(1+x)/(1 + x$sqr))$=$$$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.272198261288)/0.272198261288); newline(30, 1);

answer := haaviequadrature(y,0,1,eps,(y+y)/(1+y×y),m max);
writetext(30, {{c}integral(0, 1, (x+x)/(1+x$times$x))$=$$$$$$});
write(30, format({-d.ddddddddddd}), answer);
writetext(30, {:$RELATIVE$ERROR:$});
write(30, format({-d.dddddddddddº+dd}), (answer-0.693147180559)/0.693147180559); newline(30, 1);

FINISH:

newline(30, 1);
end
|