Iterative Adaptive Simpson and Lobatto Quadrature

Iterative Adaptive Simpson and Lobatto Quadrature

Description

In almost every standard book on numerics quadrature algorithms like the adaptive Simpson or the adaptive Lobatto algorithm are presented in a recursive way. The benefit of the recursive programming is the compact and clear representation. However, recursive quadrature algorithms might be transformed into iterative quadrature algorithms without major modifications in the structure of the algorithm.

We present iterative adaptive quadrature algorithm (adaptiveSimpson and adaptiveLobatto), which preserves the compactness and the clarity of the recursive algorithms (e.g. quad, quadv, and quadl). Our iterative algorithm provides a parallel calculation of the integration function, which leads to tremendous gain in run-time, in general. Our results suggest a general iterative and not a recursive implementation of adaptive quadrature formulas, once the programming language permits parallel access to the integration function.

Acknowledgements

Iterative Adaptive Simpson And Lobatto Quadrature inspired this file.

This file inspired Iterative Adaptive Simpson And Lobatto Quadrature.

• MATLAB release: MATLAB 7.4 (R2007a)

[Martin]

Hi Matthias,
I have tried to use your codes to improve the speed of my functions.
Though if I use the example from your codes:
Q = adaptiveSimpson(@(x) [-cos(50*x); sin(x)], 0, pi, 'tol', 1e-6)
I get the following error:
Undefined function 'adaptiveSimpson' for input arguments of type 'function_handle'.

How come?

[Stephan Koehler]

this works much faster than MATLAB's quadv, because quadv does not vectorize (i.e. it evaluates the function using only single values).

[Matthias Chung]

Reply to the last comment:

It is very easy to pass arguments to the "anonymous" function. As I recall correctly, since Matlab 6 the "@" function is exactly doing what you want. Here is an example

y = c1; f1 = @(x) f(x,y);
int1 = adaptiveLobatto(f1,a,b);

x = c2; f2 = @(y) f(x,y);
int2 = adaptiveLobatto(f2,a,b);

This makes your code more clear than passing all the arguments to the "anonymous" function.

Best,

Matthias

[pedro pablo perez velasco]

Very nice. Congratulations.

It should be very helpful to permit passing arguments to the "anonymous" function (the one that defines the function to be integrated) very much like quad, quadl and the like. For example, to perform integration in several variables.

Regards,
P.P.

[Matthias Conrad]

I updated the files. The requested H1-line have been added. An example in the comment-lines show how to handle algorithm parameters.

[Matthias Conrad]

Reply to the last comment:

It is easy to change the parameters, please use
adaptiveLobatto(fcn, a, b, 'name_of_parameter',value)

example:

Q = adaptiveLobatto(@(x) [-cos(50*x); sin(x)], 0, pi, 'tol', 1e-6)

I will update some comments in the code soon.
Thanks
Matthias

[Matt Fig]

For the adaptiveLobatto function, extra arguments aren't recognized. For example:

f = @(x) 25*exp(-25*t);
int = adaptiveLobatto(f,0,10)
returns the same result as
int = adaptiveLobatto(f,0,10,1e-14)

If, however, I go in and change the tolerance manually, on line 53 (default values), this has the intended outcome. There is something funny about the way varargin is handled. I didn't check the other function.

[John D'Errico]

These seem to be well written codes, often clearly out-performing the tools provided by MATLAB (quad, quadv, quadl.) Well done, and my compliments to the authors for providing these tools.

I'll make only one plea - to add H1 lines to your codes. While the authors will remember the names of these functions next year, or next month, when I want to find them, I won't remember the names. I'll try to use the lookfor utility, that will search through the first line of every function for a logical keyword that I'll provide. But, since the authors have not provided an H1 line here, lookfor will fail. So, I'll never be able to recall the names of these functions when I need them.

Enable the facilities in MATLAB that help you and your users to use it most efficiently. It takes only a minute more to write that H1 line (the very first line of comment in the help block.)

I'll choose not to downgrade these tools for that lack, hoping instead that the authors will repair the flaw.