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.
(Adaptive Quadrature - Revisited, ETH Zürich, Departement Informatik,
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, Boole,LobattoGautschi,Adaptive; # [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] Boole := proc(f,a,b,fa,fm,fb) local y1,y2,h,quad,err; h := (b-a)/4; y1 := fa+fb; y2 := f(a+h)+f(b-h); quad := h*(14*y1+64*y2+24*fm)/45; err := h*( -y1+ 4*y2- 6*fm)/1890; [evalf(quad), evalf(abs(err))] end: # 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: Adaptive := proc(f,a,b,fa,fm,fb,abstol) local m,Q,h,fm1,fm2; 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); Q := Adaptive(f,a,m,fa,fm1,fm,0.9*abstol)+Adaptive(f,m,b,fm,fm2,fb,0.9*abstol); 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)); Q := Adaptive(f,a,b,fa,fm,fb,abstol); 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)); Q := AdaptQuad(T[i],R[i],ValidDigits,false,"B"): msg(Q[1],I[i],"B"); Q := AdaptQuad(T[i],R[i],ValidDigits,false,"L"): msg(Q[1],I[i],"L"); od; print(Int(sin(x)-1,x=Pi/2..(Pi/2)/1000),0.56922676416839818929488); AdaptQuad(x -> sin(x)-1, Pi/2..(Pi/2)/1000, 7, true, "L"): end:
TestAdaptQuad(12); # We request 12 valid decimal digits.
![[Math]](temp1.gif)
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
![[Math]](temp2.gif)
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
![[Math]](temp3.gif)
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
![[Math]](temp4.gif)
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
![[Math]](temp5.gif)
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
![[Math]](temp6.gif)
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
![[Math]](temp7.gif)
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
![[Math]](temp9.gif)
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
![[Math]](temp10.gif)
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
![[Math]](temp11.gif)
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
![[Math]](temp12.gif)
[ 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