```KEYWORDS:  Numerical Integration, Adaptive Quadrature, Boole's Rule,
Extrapolated Simpson Rule, Lobatto Rule, Quadrature Formula, Maple. ```

The following routine for numerical integration is designed for
low accuracy, say in the range of 6 to 16 exact decimal digits.

There are two procedures, one is based on Boole's quadrature rule
and the other is a Lobatto formula with a Kronrod extension which
was designed by Walter Gautschi and Walter Gander.

Institut für Wissenschaftliches Rechnen, tech-report 306, 1998)

```AdaptQuad := proc(
###########################################################
f::procedure,    # Integrand
interval::range, # Interval of integration
valDig::posint,  # Number of valid decimal digits requested
trace::boolean,  # Trace the computation?
formula::string  # Name of quadrature formula: "B" or "L"
)##########################################################

local Q,abstol,a,b,fa,fb,fm,h,q,subdiv,savedDigits,

# [x1,x2,x3,x4,x5]    = [0,1/4,1/2,3/4,1]
# [w1,w2,w3,w4,w5;ws] = [7,32,12,32,7,90]
# [e1,e2,e3,e4,e5;es] = [-1,4,-6,4,-1,7560]

h := (b-a)/4; y1 := fa+fb; y2 := f(a+h)+f(b-h);
err  := h*(  -y1+ 4*y2- 6*fm)/1890;

# sq6 = sqrt(6)/3, sq5 = sqrt(5)/5.
# [x1,x2,x3,x4,x5,x6,x7]    = [0,(1-sq6)/2,(1-sq5)/2,1/2,(1+sq5)/2,(1+sq6)/2,1]
# [w1,w2,w3,w4,w5,w6,w7;ws] = [ 77,432, 625,672, 625,432, 77; 2940]
# [e1,e2,e3,e4,e5,e6,e7;es] = [-28, 72,-100,112,-100, 72,-28;  245]

LobattoGautschi := proc(f,a,b,fa,fm,fb) local sq6,sq5,h,m,y1,y2,y3,q1,q2;
h := b-a; m := (a+b)/2; sq5 := h*sqrt(5)/10; sq6 := h*sqrt(6)/6;
y1 := fa + fb; y2 := f(m-sq5)+f(m+sq5); y3 := f(m-sq6)+f(m+sq6);
q1 := (y1+5*y2)*h/12; q2 := (77*y1+432*y3+625*y2+672*fm)*h/2940;
[evalf(q2), evalf(abs(q2-q1))] end:

Q := `if`(formula="B", Boole(f,a,b,fa,fm,fb), LobattoGautschi(f,a,b,fa,fm,fb));
if abstol < Q[2] then subdiv := subdiv + 1;
m := (a+b)/2; h := (b-a)/4; fm1 := f(a+h); fm2 := f(b-h);
fi;
if trace then printf("[%14.8f, %14.8f] -> %18.16f\n", evalf(a),evalf(b),Q[1]) fi;
Q end:

# --- Main begins here! ---
savedDigits := Digits; Digits := 42; subdiv := 0;
a := op(1, interval); b := op(2, interval);
fa := f(a); fb := f(b); fm := f((a+b)/2);
h := b-a; q := h*(fa+fb+fm+f(a+0.2311*h)+f(a+0.4860*h)+
f(a+0.6068*h)+f(a+0.8913*h)+f(a+0.9501*h))/8;
abstol := evalf(0.01*10^(-valDig)*abs(q));

print(`Integral: `,evalf(Q[1],valDig),`Error-estimate: `,evalf(Q[2],4),
`Subdivisions: `, subdiv); Digits := savedDigits;
Q end:
```

```TestAdaptQuad := proc(ValidDigits)

T1  := proc(x) `if`(x=0,1,sin(x)/x) end; R1 := 0..3; I1 := Si(3);
T2  := proc(x) sqrt(x*(4-x)) end; R2 := 0..2; I2 := Pi;
T3  := proc(x) exp(x)*cos(x) end; R3 := 0..Pi; I3 := -1/2*exp(Pi)-1/2;
T4  := proc(x) sqrt(x) end; R4 := 0..1; I4 := 2/3;
T5  := proc(x) exp(-(x*x)) end; R5 := 0..30; I5 := 1/2*erf(30)*sqrt(Pi);
T6  := proc(x) `if`(x=0,0, log(x)) end; R6 := 0..1; I6 := -1;
T7  := proc(x) sqrt(x*x+power(10,-10)) end; R7 := -3..5; I7 := 17;
T8  := proc(x) `if`(x=0,0,cos(log(x))) end; R8 := 0..1; I8 := 1/2;
T9  := proc(x) sin(sqrt(x)) end; R9 := 0..1; I9 := 2*sin(1)-2*cos(1);
T10 := proc(x) `if`(x=0 or x=1,0,sqrt(x)/(x-1)-1/(log(x))) end;
R10 := 0..1; I10 := 2-gamma-2*ln(2);
T11 := proc(x) exp(x^2)*(1+erf(x)) end;
R11 := -8..2; I11 := 33.659167338320206189432;
T12 := proc(x) (x*(x-88)*(x+88)*(x-47)*(x+47)*(x-117)*(x+117))^2 end;
R12 := -256..256; I12 := 0.1310268955222668656627*10^29;

T := [T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12];
R := [R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12];
I := [I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12];

msg := proc(q,i,n) name := `if`(n="B",`Boole   `,`LobattoG`);
edd := -log[10](evalf(abs(1-q/i),60));
valid := evalb(edd >= ValidDigits);
if valid then print(name,`: Success! `,evalf(edd,4),` valid decimal digits.`);
else print(name,`: Failed! Only `,evalf(edd,4),` valid decimal digits but `,
ValidDigits,` requested.`);
fi; end;

for i from 1 to nops(T) do
print(Int(T[i](x),x=R[i]),` = `, evalf(I[i],ValidDigits));
od;

print(Int(sin(x)-1,x=Pi/2..(Pi/2)/1000),0.56922676416839818929488);
end:
```

```  Boole   :  Success!  16.55 valid decimal digits.
Integral:  1.84865252800 Error-estimate:  .1257 10^(-12) Subdivisions: 108

LobattoG:  Success!  23.61  valid decimal digits.
Integral:  1.84865252800 Error-estimate:  .1163 10^(-12) Subdivisions: 30
```

```  Boole   :  Success!  12.59 valid decimal digits.
Integral:  3.14159265359 Error-estimate:  .7200 10^(-12) Subdivisions: 226

LobattoG:  Success!  16.06 valid decimal digits.
Integral:  3.14159265359 Error-estimate:  .3583 10^(-12) Subdivisions: 162
```

```  Boole   :  Success!  14.56 valid decimal digits.
Integral:  -12.0703463164 Error-estimate:  .3338 10^(-11) Subdivisions: 162

LobattoG:  Success!  24.21 valid decimal digits.
Integral:  -12.0703463164 Error-estimate:  .9952 10^(-12) Subdivisions: 61
```

```  Boole   :  Success!  12.63 valid decimal digits.
Integral:  .666666666667 Error-estimate:  .1456 10^(-12) Subdivisions: 212

LobattoG:  Success!  16.14 valid decimal digits.
Integral:  .666666666667 Error-estimate:  .6880 10^(-13) Subdivisions: 159
```

```  Boole   :  Success!  12.86 valid decimal digits.
Integral:  .886226925453 Error-estimate:  .5565 10^(-12) Subdivisions: 192

LobattoG:  Success!  16.66 valid decimal digits.
Integral:  .886226925453 Error-estimate:  .1687 10^(-12) Subdivisions: 109
```

```  Boole   :  Success!  14.54 valid decimal digits.
Integral:  -1.00000000000 Error-estimate:  .1708 10^(-12) Subdivisions: 558

LobattoG:  Success!  17.39 valid decimal digits.
Integral:  -1.00000000000 Error-estimate:  .7552 10^(-13) Subdivisions: 373
```

```  Boole   :  Success!  12.53 valid decimal digits.
Integral:  17.0000000014 Error-estimate:  .7559 10^(-13) Subdivisions: 31

LobattoG:  Success!  17.21 valid decimal digits.
Integral:  17.0000000014 Error-estimate:  .1598 10^(-12) Subdivisions: 36
```

```  Boole   :  Success!  12.62 valid decimal digits.
Integral:  .602337357879 Error-estimate:  .1349 10^(-12) Subdivisions: 219

LobattoG:  Success!  16.09 valid decimal digits.
Integral:  .602337357880 Error-estimate:  .6748 10^(-13) Subdivisions: 160
```

```  Boole   :  Success!  13.40 valid decimal digits.
Integral:  .0364899739786 Error-estimate:  .4821 10^(-14) Subdivisions: 399

LobattoG:  Success!  17.10 valid decimal digits.
Integral:  .0364899739786 Error-estimate:  .2464 10^(-14) Subdivisions: 280
```

```  Boole   :  Success!  13.77 valid decimal digits.
Integral:  33.6591673383 Error-estimate:  .4157 10^(-10) Subdivisions: 262

LobattoG:  Success!  18.08 valid decimal digits.
Integral:  33.6591673383 Error-estimate:  .1844 10^(-10) Subdivisions: 135
```

```[    1.57079633     1.37464314] ->  .0012554497912985
[    1.37464314     1.17848994] ->  .0087303366026269
[    1.57079633     1.17848994] ->  .0099857863939254
[    1.17848994      .98233675] ->  .0233934268989034
[     .98233675      .78618356] ->  .0446823491543437
[    1.17848994      .78618356] ->  .0680757760532471
[    1.57079633      .78618356] ->  .0780615624471725
[     .78618356      .59003037] ->  .0717806122030891
[     .59003037      .39387718] ->  .1036489200148414
[     .78618356      .39387718] ->  .1754295322179305
[     .39387718      .00157080] ->  .3157356695032952
[     .78618356      .00157080] ->  .4911652017212257
[    1.57079633      .00157080] ->  .5692267641683982

Integral:  .5692268 Error-estimate:  .2204 10^(-9) Subdivisions: 13

```