Gauss-Lobatto-Kronrod formulae and adaptive numerical integration Radu T. Trˆ ımbit ¸as ¸ Department of Applied Mathematics “Babes ¸-Bolyai” University Cluj-Napoca, Romania Email: tradu@math.ubbcluj.ro Maria Gabriela Trˆ ımbit ¸as ¸ Department of Computer Science “Babes ¸-Bolyai” University Cluj-Napoca, Romania Email: gabitr@cs.ubbcluj.ro Abstract—The aim of this paper is to develop a MATLAB function for one dimensional numerical integration based on adaptive algorithms and Gauss-Lobatto-Kronrod formulas. Using Maple, we find a triple of formulas, and then we use it to code an analogous of MATLAB quadl function with a higher degree of exactness. Finally, some examples and tests which compare our function and quadl are given. Our function is a good alternative to quadl when the accuracy and reliability requirements are hard. I. I NTRODUCTION We want to approximate I =  b a f (x)dx, where f :[a, b] → R is an integrable function, and the interval [a, b] is bounded. Adaptive algorithms [1] decide dynamically how many function evaluations are needed. The efficiency and reliability of such algorithms depends upon the subdivision strategy. We shall integrate f using two methods which provide us the approximations I 1 and I 2 . If the difference of this two approximations is less than a given tolerance, we accept the better of them, say I 2 , as approximate value of integral. Otherwise, we divide [a, b] into two (or more) subintervals; then proceed recursively on each part. Carl deBoor [2] and William Kahan [8] showed that, it is impossible to construct a correct program that integrates each given function. Hence, the task of a software developer is to code programs which function correctly for a class of function as large as possible. The papers [3], [4] describe the basic ideas of MATLAB quad and quadl. They also emphasize the importance of a good termination criterion. Maple is used as a tool to find convenient quadrature formulas for integral and error estimation. In this paper we use these ideas to find a combination of formulae with a higher degree of exactness and to turn it into an adaptive algorithm and a MATLAB 1 implementation. We hope our function behaves better when the accuracy requirements are hard. 1 MATLAB is a trademark of Mathworks Inc, Natick, MA II. THE TRIPLE OF FORMULAE We look for a combination of three quadrature rules • A Gauss-Lobatto formula with five nodes for the evalua- tion of the integral • A Kronrod extension of the previous formula with nine nodes for the error estimation; • A Kronrod extension of the previous extension to esti- mate how much more accurate is the previous extension compared to the basic rule. For details on Gauss-Kronrod extensions see [6]. The (variable) nodes are the roots of an orthogonal polynomial with a convenient weight. Once we have the nodes, we can compute the coefficients using the exactness of formula for x k , k =0,...,d, where d is the degree of exactness. In order to generate the triple of formulae, we coded a Maple procedure to compute the nodes and a Maple procedure to compute the coefficients. Here is the code. > GN:=proc(we,deg,x) > local c,pk,k,ec,d1,sc,nodes; > description "find nodes"; > d1:=deg-1; > pk:=sum(c[j] * xˆj,j=0..d1)+xˆdeg; > for k from 0 to d1 do > ec[k]:=int(pk * xˆk * we,x=-1..1)=0: > end do; > sc:=simplify(solve({seq(ec[j], > j=0..d1)},{seq(c[j],j=0..d1)})); > assign(sc); > nodes:=solve(eval(pk),x); > return nodes; > end proc; > GC:=proc(nodes) > local k,d1,eco,co,F,Qk; > description "find coefficients"; > d1:=nops(nodes)-1; > Qk:=0; > for k from 0 to d1 do > Qk:=Qk+F[k] * f(nodes[k+1]); > end do; > Qk:=unapply(Qk,f); > for k from 0 to d1 do > eco[k]:=Qk(x->xˆk)=int(xˆk,x=-1..1); > end do; > co:=solve([seq(eco[j],j=0..d1)], > [seq(F[j],j=0..d1)]); > assign(co); > return [seq(F[j],j=0..d1)]; > end proc; The interior nodes of the basic Gauss-Lobatto formula are the roots of the polynomial π ∈ P 3 , orthogonal on [−1, 1] with respect to the weight w(1) = (1−x) α (1+x) β ,