Paper Menu >>
Journal Menu >>
Interval IntegrationRevisited Sérgio Galdino Polytechnic School ofPernambucoUniversity Rua Benfica, 455- Madalena 50.750-470 - Recife- PE- Brazil Email: galdino@poli.br 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,rdenotes 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=1xi 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 −123 25 cosh(x)−cos(x)dx =0.47942822668880166736 f3=1 −1x4+x2+0.9−1dx =1.5822329637296729331 f4=1 01+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 −11.005 +x2−1dx=1.5643964440690497731 f13 =1 01+(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 |