Interval IntegrationRevisited
Sérgio Galdino
Polytechnic School ofPernambucoUniversity
Rua Benfica, 455- Madalena
50.750-470 - Recife- PE- Brazil
Abstract—We present an overview of approaches to self-
validating one-dimensional integration quadrature formulas and
a verified numerical integration algorithm with an adaptive
strategy. The new interval integration adaptive algorithm delivers
a desired integral enclosure with an error bounded by a specified
error bound.The adaptive technique isusually much more
efficient than Simpson’s rule and narrow interval results can
be reached.
Index Terms—Interval computing, self-validating methods,
numerical integration & interval arithmetic.
I. INTRODUCTION
Numerical integration(quadrature) formulas are ofthe form
I=b
a
f(x)dx =
n
∑
i=0
w(n)
if(xi)+E,
with a=x0<x1<x2<···<xn=band weights w(n)
i,itis
possible toderive aformulafor thelocalerrorof themethod
based on a higher derivative off(x)at an intermediate point
ξ as
E=C·f(k)( ξ ),where ξ ∈[a,b],
thefactorCdepends onthemethod and usuallyisa power of
the width b−a.
Automatic integration programs mayprint out the ‘theo-
retical’ error which has been achieved and which is used to
provideshutoff.Exist somepitfalls ofpracticalworkingof
automatic integration programs.The tolerance requested may
be impossible to meet. When more and more points are called
the roundoff error enter to worsen progressively the result. On
theoretical side, the estimate of accuracy validity depends on
theoretical information about the function such as theoretically
accurate bounds onderivatives,monotonicity, convexity,etc.
With interval analyses, it is possible to determine veri-
fied bounding of analyticallyderived quadratureformulas for
integration rules. The evaluation code derivativef(k)may
be obtained from automatic differentiation. Different from
classical schemes, they do not require an a-priory derivation
of analytical error bounds, furthermore the rigorous bounds
are calculatedautomatically inparallel tointegralcomputation
allowing better errorcontrol.
Adaptive quadrature is another area in which interval meth-
ods isuseful.This isbecause,due to theform of the error
term, meaningful interval enclosures for the actual integral
can easilybe computed. Replacingheuristic errorestimates by
these rigorousenclosuresresultsina quadraturealgorithmthat
produces guaranteed boundson the actualintegral.A common
interval algorithm istointegrate, with automaticdifferentiation
software, high-orderand variable degree Taylorpolynomial
approximations tof[1],[2], [3].
The next section briefly introduces interval arithmetics.
Section 3describesthe interval integration andpresents the
approach adopted tothe respectiveinterval integrations. Sec-
tion 4 presents numerical experiments and the related analysis
results. Section 5concludes.
II. INTERVAL ARITHMETIC
A form of interval arithmetic perhaps first appeared in 1924
and 1931 in [4], [5], then later in [6]. Modern development of
interval arithmetic began with R. E. Moore’s dissertation [7] as
a methodfordetermining absoluteerrors ofanalgorithm, con-
sidering all data errorsand rounding, afterR.E. Moore intro-
duced interval analysis [8]. Interval arithmetic is an arithmetic
defined onsets of intervals, rather thansets of realnumbers.
The power of interval arithmetic lies in its implementation
on computers. In particular, outwardly rounded computations
allows rigorous enclosures for the ranges of operations and
functions.
A. Notation
Throughout this paper, all scalar variables is denoted by
ordinary lowercase letters (a).Interval variables are enclosedin
square brackets [a]. Underscoresand overscores denote lower
and upper bounds, respectively. Angle brackets ,are used
for defining intervals by midpoints andradius.
A real interval [x]isa nonempty setof realnumbers
[x]=[x,x]={˜x∈R:x≤˜x≤x}(1)
where xand xare called the infimum (inf) and supremum
(sup)
, respectively, and˜xis a point value belonging to an
interval variable [x].
The set ofall intervals Ris denotedby I(R)where
I(R)={[x,x]:x,x∈R:x≤x}(2)
The midpoint of[x]isdefined as,
mid[x]=1
2(x+x),(3)
the width of[x]isdefined as,
wid[x]=( x−x)(4)
Open Journal of Applied Sciences
Supplement：2012 world Congress on Engineering and Technology
108 Cop
y
ri
g
ht © 2012 SciRes.

and the radiusof [x]is definedas,
rad [x]=1
2(x−x)=1 2wid[x](5)
The radius may also be used to define an interval [x]∈I(R).
Then m,rdenotes an interval with midpoint mand radius
r. The [x,x]≡xis called pointinterval or thin interval. A
point or thin interval has zero radius and a thick interval has
a radius greaterthan zero.
III. INTERVAL INTEGRATION
From mean value theorem there exists a ξ ∈[x]=[x,x]
such that
x
x
f(x)dx =wid ([x]) f( ξ ).
Splitting the interval [x] into nsubinterval with break points
{xi},a=x0<x1<···<xn−1,xn=b,then
I=
n
∑
i=1xi
xi−1
f(x)dx =
n
∑
i=1
wid([x]i)f( ξ
i),
where [x]i=[xi−1,xi]and
ξ
i∈[x]i.
Let [f]([x]) bean interval extension off(x)on [x]. Then
Ii≡w([x]i)f(
ξ
i)⊆[R]i≡w([x]i)[ f]([x]i),
and we have an interval inclusion of I,
I⊆[R]≡
n
∑
i=1
[R]i(6)
The width of [R]measures the quality of determined integral
given bounds to both roundingerrors andtruncation errors.
A.IntervalTrapezium Rule
If f(x)istwo timescontinuously differentiable, then
Ti ≡wid ([x]i)
2(f(xi)+ f(xi)) −wid ([x]i)3 12 f(2)(
ξ
i),(7)
where
ξ
i∈[x]i. If we useinterval arithmeticextension[f],
[f](2), thenwe usethese to get
I⊆[T]≡
n
∑
i=1
[T]i(8)
with
[T]i=wid([x]i)
2([ f](xi)+[f](xi)) −wid ([x]i)3 12 [f](2)([x]i).
(9)
Note that thecall of [f] with thininterval arguments return
intervals that includesrounding errors.
B.IntervalSimpson1/3
If f(x) isfour timescontinuously differentiable, then
Si=wid([x]i)
6(f(xi)+4f(mid([x]i)+ f(xi))
−wid([x]i)5
2880 f(4)(
ξ
i)(10)
where ξ
i∈[x]i. If we use interval arithmetic extension [f],
[f](4), thenwe usethese to get
I⊆[S]≡
n
∑
i=1
[S]i(11)
with
[S]i=wid([x]i)
6([ f](xi)+4[f](mid([x]i)+[f](xi))
−wid([x]i)5
2880 [f](4)([x]i)(12)
Note that the call of [f] with thin interval arguments return
intervals that includesrounding errors.
C.Adaptivequadrature
The drawback with any algorithm based on the composite
quadrature rules is that it makes no attempt to respond to
the localformof the integrand. The objectiveof an adaptive
scheme is to use an unequal mesh spacing and todeterminethe
size ofeachsubinterval so astosatisfytheoverall accuracy
requirement withthe minimum numberof subintervals(and
consequently the minimumnumber ofevaluations of the
integrand).
Totalapproximate integral= Sum of approximateintegrals
over subintervals.
Total estimated error = Sum of estimated errors over subin-
tervals:
•Locally Adaptive Quadrature:
Errorforeachsubinterval ≤globalerrortarget ×length of subinterval
b−a
Note thatlocally adaptivequadrature is anatural example of recursion.
•Globally Adaptive Quadrature:
Errorforallsubintervals ≤globalerror target
Adaptive algorithms are just as efficient and effective as
traditional algorithms for "well-behaved" integrands, but can
be also effective for "badly-behaved" integrands for which
traditional algorithms fail.
D.AdaptiveSimpson’s method
Adaptive Simpson’s method, alsocalledadaptive Kuncir’s
Simpson rule, wasproposed by G.F. Kuncir in 1962 [9].
Adaptive Simpson’smethod uses an estimate of the error we
get fromcalculatingadefiniteintegralusingSimpson’s rule.If
the error exceeds auser-specified tolerance, the algorithm calls
for subdividing the interval of integration in twoand applying
adaptive Simpson’s method to each subinterval in a recursive
manner. A criterion fordetermining when to stop subdividing
an interval, suggested by J.N. Lyness[10], is
|S(a,c)+S(c,b)−S(a,b)|/15 < τ (13)
where [a,b]is an interval with midpoint c,S(a,c),S(c,b),
and S(a,b)are the estimatesgiven bySimpson’s rule on the
corresponding intervals and τ is the desired tolerance for the
interval.
Cop
y
ri
g
ht © 2012 SciRes.109

E.AdaptiveIntervalSimpson’s method
If we are using recursion, the implementation of an
interval algorithmis simple. Thefollowing Matlabwith the
Intlab Toolbox [12]functionis basedonadaptsimp.m [13].
1.function int = Iadaptsmp(f,f4,X,tol,lev,fa,fm,fb)
2.global n
3.global D
4.%IADAPTSMP Interval Adaptive Simpson 1/3 rule
5.echo on
6.if nargin == 4
7.lev = 1;
8. a=inf(X);
9. b=sup(X);
10. D=diam(X);
11. fa = feval(f,a);
12. fm = feval(f,(a+b)/2);
13. fb = feval(f,b);
14. n=3;
15. end
16. % recursive calls start here
17. % checking for too many levels of recursion
18. if lev > 1000
19. disp(’1000 levels of recursion reached.’)
20. X=infsup(a,b);
21. wx=diam(X);
22. int = (b-a)*(fa+4*fm+fb)/6 -wx^5/2880*f4(X);
23. else
24. % Divide the interval in half and apply interval
25. % Simpson 1/3 rule on each half. As an verified
26. % error estimate: rad(int) > 2*tol*h/D;
27. a=inf(X);
28. b=sup(X);
29.h=b-a;
30. wx=h/2;
31. flm = feval(f,a+h/4);
32. frm = feval(f,b-h/4);
33. n=n+2;
34. X=infsup(a,a+h/2);
35. simpl = h*(fa + 4*flm + fm)/12 - wx^5/2880*f4(X);
36. X=infsup(a+h/2,b);
37. simpr = h*(fm + 4*frm + fb)/12 - wx^5/2880*f4(X);
38. int = simpl + simpr;
39. X=infsup(a,b);
40. % err =rad(int) > 2*tol*h/D
41. % tolerance not satisfied, recursively refine
42. if rad(int) > 2*tol*h/D;
43.m=(a+b)/2;
44. XL=infsup(a,m);
45. XR=infsup(m,b);
46. int = Iadaptsmp(f,f4,XL,tol,lev+1,fa,flm,fm) ...
47. + Iadaptsmp(f,f4,XR,tol,lev+1,fm,frm,fb);
48. end
49. end
Sample Matlab M-flle with the Intlab Toolboxthat calls the
above functionis followed by running results:
1.format long
2.% intvalinit(’DisplayMidRad’)
3.global n
4.f = @(x)23/25*cosh(x)-cos(x)
5.f4 = @(x)23/25*cosh(x)-cos(x)
6. X=infsup(-1,1)
7. err=1.0E-12;
8. IADAPTSMP(f,f4,X,err)
9. n
10. err1=rad(ans)
>> ADAPTXI
intval ans =
<0.47942822668878, 0.00000000000047>
n=
241
err1 =
4.607425552194400e-013
>>
IV. NUMERICAL EXPERIMENTS
In thefollowing, wewill use the“battery" test offunctions
which are a subset ofthe set used by Gonnet [11].The
battery of testfunctions arelistedin Table I. The resulthas
been calculatedusing 20 digits ofprecision with theMathcad
computer software.
The main mechanism for get a narrow interval results is
to increasethe numberof subdivisions. Nextwe considerthe
progression of the discretizationerror as we increase n.Table
II shows that in the early stages increasing nimproves the
discretization.
Weobserve however that if the number ofsubdivisions were
extremelylarge thismightlead tosignificant round-of (more
terms in the sum, each with a round-off error to contribute).
Verified trapezoidal rulewith105subdivisions is widerthan
104subdivisions andverified Simpson1/3 rulewith 104
subdivisionsis wider than 103. Anarrow interval implies high
precision; we can specify plausible values to within a tiny
range. Verified stepruleinterval results arenarrowing (from
10 to 106subdivisions), but intervals are wider, compared
with verifiedtrapezoidal and Simpson1/3 rules, impliespoor
precision.
In Table III we present the resultsof numerical experiments
with interval adaptive Simpson 1/3 rule. The battery of
test functions are listed in Table I. Forthesetest functions
we will consider the usual metrics of efficiency,i.e. number
of function evaluations required for a given accuracy.We
have tested the code on four radius interval tolerances τ =
10−3,10−6,10−9,10−12.
Interval adaptiveSimpson 1/3 algorithm was moreefficient
than interval Simpson 1/3, the only exception is withf1for
τ =10−3. Should be stressed that interval tolerance τ =10−12
is not reached, by the interval Simpson 1/3 algorithm, for
f9,f10,f11 test functions,with theused precision.
V. CONCLUSION
The design ofintervaliterativeformulas for self-validating
one-dimensional integration quadratureformulas isvery im-
portant andisalsoaninterestingtaskininterval computing.
In thispaper, wehaveproposedanew (supposed) one-
dimensional interval integration adaptive algorithm, toget
rigorous boundsof adesired integral. Numerical testsdemon-
strate that interval adaptive technique is usually more efficient
than interval Simpson’s rule with narrow interval results.
Finding self-validating one-dimensional definite integration by
derivative-free interval methods shouldbe consideredinfuture
studies.
110 Cop
y
ri
g
ht © 2012 SciRes.

TABLE I
THEFUNCTIONSUSEDFOR THE “BATTERY”TEST.
f1=1
0exdx=1.7182818284590452354
f2=1
−123
25 cosh(x)−cos(x)dx =0.47942822668880166736
f3=1
−1x4+x2+0.9−1dx =1.5822329637296729331
f4=1
01+x4−1dx =0.86697298733991103757
f5=1
02(2+sin(10
π x))−1dx=1.154700538379251529
f6=1
0(1+x)−1dx =0.69314718055994530942
f7=1
0(1+ex)−1dx =0.37988549304172247537
f8=1
0.1sin(100
π x)/( π x)dx=0.0090986375391668429156
f9=10
0
√50 e−50
π x2dx=0.5
f10 =10
025 e−25xdx =1.0
f11 =10
050
π
2500x2+1−1dx=0.49936338107645674464
f12 =1
−11.005 +x2−1dx=1.5643964440690497731
f13 =1
01+(230x−30)2−1dx=0.013492485649467772692
TABLE II
ESTIMATES OF 1
−1cosh(x)23
25 −cos(x)dx=sinh(x)23 25 −sin(x)1
−1
.
Analytical Answer Bounds =[0.47942822668880,0.47942822668881]
Step Rule withVerification
SubdivisionsVerifiedResult
10 0.492244876649698,0.191866375632378 
100 0.479556403654468,0.019186637563241 
1000 0.479429508459513,0.001918663756343 
10000 0.479428239506509,1.918663758142536e−004
100000 0.479428226816973,1.918663937550136e−005
1000000 0.479428226690078,1.918681864720995e−006
Trapezoidal Rule withVerification
SubdivisionsVerifiedResult
10 0.479421870930105,6.395545854416262e−004
100 0.479428226049601,16.395545891768606e−007
1000 0.479428226688739,6.395906027023557e−010
10000 0.479428226688802,1.002808946992673e−012
100000 0.479428226687608,3.629152534045943e−012
1000000 0.479428226671977,3.620936883663717e−011
Simpson 1/3 Rule with Verification
SubdivisionsVerifiedResult
10 0.479428217025575,2.131848625963606e−007
100 0.479428226688792,2.139066701545289e−012
1000 0.479428226688796,7.244205235679146e−014
10000 0.479428226688723,7.243095012654521e−013
100000 0.479428226688001,7.241873767327434e−012
1000000 0.479428226680795,7.241879318442557e−011
REFERENCES
[1]G.F.CorlissandL. B.Rall:Adaptive, self-validating numerical quadra-
ture. SIAM J. Sci. Statist. Comput. 8(5):831–47,(1985)
TABLE III
RESULTS OF BATTERY TEST:ADAPTIVEAND SIMPSON 1/3RULESWITH
VERIFICATION
Algorithms:
AdaptiveSimpson 1/3
(Simpson 1/3)
f(x)
τ =10−3 τ =10−6 τ =10−9 τ =10−12
f15(3)9(12)33 (39)129(150)
f25(6)17 (21)65(78)241 (306)
f325 (30)97(108)361 (423)1441 (1701)
f49(12)33 (45)109(174)429 (693)
f5105 (135)393(522)1489(2064)5921 (8718)
f65(6)13 (18)49(63)189 (252)
f75(6)17 (18)65(69)257 (273)
f8445 (516)1797(2070)7077 (8241)28125 (32922)
f949 (684)133(2640)449(10488)1725 (−)
f10 49 (531)125(2106)433 (8376)1681 (−)
f11 89 (2916)305(10833)1193 (42870)4765 (−)
f12 17 (24)57(84)241 (330)945 (1323)
f13 57 (519)153(1872)541 (7323)2161 (29193)
TABLE IV
RESULTS OFBATTERYTEST:ADAPTIVESIMPSON 1/3RULEWITH
VERIFICATION
Integral AdaptiveSimpson1/3 ( τ =10−12 )
fVerifiedResultError
f1<1.71828182845904,0.00000000000029 >0.29 ×10−12
f2<0.47942822668878,0.00000000000047 >0.47 ×10−12
f3<1.58223296372967,0.00000000000048 >0.48 ×10−12
f4<0.86697298733961,0.00000000000072 >0.72 ×10−12
f5<1.15470053839294,0.00000000000056 >0.56 ×10−12
f6<0.69314718055994,0.00000000000056 >0.56 ×10−12
f7<0.37988549304172,0.00000000000019 >0.19 ×10−12
f8<0.00909863753917,0.00000000000044 >0.44 ×10−12
f9<0.50000000000000,0.00000000000004 >0.04 ×10−12
f10 <0.99999999999999,0.00000000000008 >0.08 ×10−12
f11 <0.49936338107645,0.00000000000062 >0.62 ×10−12
f13 <1.56439644405124,0.00000000000057 >0.57 ×10−12
f13 <0.01349248564896,0.00000000000054 >0.54 ×10−12
[2]R.Kelch. Ein adaptivesVerfahren zur numerischen Quadraturmit
automatischer Ergebnisverifikation.PhDthesis, Universität Karlsruhe,
(1989)
[3]U.Storck. Numerical integration in two dimensions with automatic result
verification. InE. Adams and U. Kulisch, editors,Scientific Computing
with Automatic Result Verification,AcademicPress, New York, etc.,
187–224, (1993) .
[4]J.C. Burkill:Functions of intervals. Proceedings ofthe LondonMathe-
matical Society, 22:375-446, (1924)
[5]R.C. Young:The algebra of many-valuedquantities. Math. Ann. 104:260-
290, (1931)
[6]T. Sunaga: TheoryofanInterval AlgebraanditsApplication to
Numerical Analysis. Gaukutsu Bunken Fukeyu-kai, Tokyo, (1958)
[7]R.E.Moore: Interval Arithmetic andAutomatic ErrorAnalysis inDigital
Computing. PhD thesis, Stanford University, October, (1962)
[8]R.E.Moore: Interval Analysis.Prentice Hall, EnglewoodClifs, NJ,USA,
(1966)
[9]G.F. Kuncir: Algorithm 103:Simpson’s rule integrator. Communications
of the ACM 5(6): 347, (1962)
[10]J.N.Lyness: Notes on the adaptive Simpson quadrature routine. Journal
of the ACM 16 (3): 483–495, (1969)
[11]P.Gonnet: Adaptivequadrature re-revisited. ETHZürich ThesisNr.
18347 (2009). http://dx.doi.org/10.3929/ethz-a-005861903
[12]S.M.Rump: INTLAB - INTerval LABoratory. Tibor Csendes, Develop-
ments in Reliable Computing, Kluwer AcademicPublishers, Dordrecht,
77–104, (1999), http://www.ti3.tu-harburg.de/rump/
[13]DouglasN.Arnold: AConcise Introductionto NumericalAnalysis.
Institute forMathematics andits Applications,Universityof Minnesota,
Minneapolis, MN 55455 (2001). http://www.ima.umn.edu/~arnold/
Cop
y
ri
g
ht © 2012 SciRes.111