Numerical integration: quadgk vs integral

Hello: I have not been able to find out what is the underlying quadrature formula in Matlab's builtin function integral. I think that quad uses Simpson, quadl a Gauss-Lobato formula and quadgk the Gauss-Kronrod formula G7K15, and all of them use some kind of adaptative scheme. What does integral(f,a,b) do with my poor function f? When should I use integral and when quadgk? Thank you, Mariano

Answer 1

[Christiaan]

Dear Mariano,

The Mathlab function 'integral(fun,xmin,xmax)' numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances.

In the documention of this function (doc integral), this ref can be found:

    %   Portions based on "quadva" by Lawrence F. Shampine.
    %   Ref: L.F. Shampine, "Vectorized Adaptive Quadrature in Matlab",
    %   Journal of Computational and Applied Mathematics 211, 2008, pp.131-140

In the Documentation of quadgk you can find a list of multiple quadrature functions that can be used whith there (dis)advantages.

Kind regards, Christiaan

[Mariano]

Dear Christiaan:

Thank you for the answer, but it does not solve my problem completely.

Both quadgk and integral are based on quadva. From Shampine's paper, we can gather that the difference with the old functions quad and quadl is that now all the work is properly vectorized. All 4 functions use global adaptive quadrature. In the documentation of quadgk, a comparison between quad, quadl and quadgk can be found. Also in the documentation of quadgk one can find an item about the Algorithm (at the end of the page):

quadgk implements adaptive quadrature based on a Gauss-Kronrod pair (15th and 7th order formulas)

Something similar could be found in the old pages for quad and quadl:

(R2011b) quad implements a low order method using an adaptive recursive Simpson's rule.

(R2011b) quadl implements a high order method using an adaptive Gauss/Lobatto quadrature rule.

What is more. In old releases the code of quad and quadl was visible, and you could tell what the specific formulas were.

What I would like to know is the specific quadrature formula used in each of the subintervals to evaluate the integral on that element (and which is used to evaluate the error, although this usually is determined by the previous one).

I would also appreciate a comparison between quadgk and integral like that of quadgk, quad and quadl.

This is mere academic interest, but I find discouraging that in more and more pages in MATLAB documentation, the Algorithm section and proper academic references are being excluded. For instance, the reference to Shampine's paper does not appear in "doc integral" or "help integral". You must "type integral" to see it.

I hope my doubt is clearer now.

Best wishes,

Mariano

[Mark Rupright]

Mariano,

It has been a while since you asked the question, but I had the same one. With some more digging, I discovered that "integral" and "quadgk" use the same algorithm; integral is an updated version renamed to make it easier for new users to find.

Cleve Moler explains it in this blog post.

I suspected this was so because the help files for both functions direct to the same paper by Shamipine, but I wish the help files made it clearer.

[TheStranger]

Strange, I tried both of them on some pretty complicated function and the answers differ by ~10^(-8) Also, for some strange reason, integral does not have an option "MaxIntervalCount", while quadgk does.

[Josh Meyer]

To summarize the differences between the two functions:

• integral and quadgk use the same quadrature method as described in the reference paper
• integral and quadgk use different default error tolerances for RelTol and AbsTol
• quadgk has the MaxIntervalCount option, and uses a relatively small default value compared to integral. This means that with integral you will rarely if ever need to adjust this parameter (hence why it isn't there). This would only be with highly oscillatory integrands that are difficult to evaluate.
• quadgk has an extra output that approxmates the absolute error. This lets you specify MaxIntervalCount, check the error in the calculation, and adjust MaxIntervalCount as needed.

Answer 2

[Mariano]

Thanks a lot for your answers. They have been very helpful.

Best wishes,

Mariano