Double Integration in MATLAB – Understanding Tolerances4

Posted by Loren Shure,

In today's post, I am joined by Mike Hosea, a developer who occasionally works on integration routines for MATLAB. In recent releases, we added new integration routines, including integral2 for double integrals. Today we'll talk about some nuances in using this routine and adjusting the absolute and relative tolerances.

Contents

Requirement for Calling Integration Functions

There is at least one requirement on the integrand when you call any integration function, and that is that it can be called with multiple inputs and it will return corresponding outputs, 1-for-1, with the inputs. So, to integrate $x * y$, you need to supply the vectorized form, i.e., x .* y.

Instead of creating a file with the integrand definition, we instead create an anonymous function. If you are unfamiliar with anonymous functions and want to learn more about them, please check the categories to the right of the blog for discussions and examples.

F0 = @(x,y) x .* y;


We can now integrate F between 0 and 1 for both directions, x and y.

Q = integral2(F0,0,1,0,1)

Q =
0.250000000000000


Exploring Tolerances

Here is another example. First let's integrate the function F1

F = @(x,y)exp(10*x).*tan(y)

F =
@(x,y)exp(10*x).*tan(y)


between 0 and 1 for x and 0 and $\pi/4$ for y.

Q1 = integral2(F,0,1,0,pi/4)

Q1 =
7.633444730985692e+02


The analytic answer for this integral is

Q0 = log(sqrt(2))/10*(exp(10) - 1)

Q0 =
7.633444758094897e+02


and you can see that Q1 differs from Q0, perhaps more than we might want. We first adjust the numeric display to see enough decimal digits, before displaying the abolute error of the integration Q1.

format long
intabserr1 = abs(Q1-Q0)

intabserr1 =
2.710920512072335e-06


One of the advantages of the integral2 routine (and its companions for other dimensions) over most of the older routines is that they provide mixed relative and absolute error control, which is to say, you can supply two tolerances instead of just one. Not only can you, you should. You actually need both, or rather sometimes the one and sometimes the other, but it is good practice to think about what you want for both of them and to be aware of the default values in case you choose not to set them. What you would rarely want to do is change only the absolute error tolerance, AbsTol.

Now suppose we want more accuracy and try to crank down the absolute tolerance.

Q2 = integral2(F,0,1,0,pi/4,'AbsTol',1e-14)
intabserr2 = abs(Q2 - Q0)

Q2 =
7.633444730985692e+02
intabserr2 =
2.710920512072335e-06


We got the same answer as before! The reason is that integral2 uses a mixture of relative and absolute error control, and it is always the least restrictive of the two that matters. The default AbsTol is 1e-10, and the default RelTol is 1e-6. The absolute error is bigger than the absolute tolerance, true, but the relative error is

intrelerr2 = abs(Q2 - Q0)/abs(Q0)

intrelerr2 =
3.551372411777180e-09


and that's smaller than RelTol, so integral2 returned that answer. Since the relative error is between 1e-8 and 1e-9, there should be about 8 correct significant digits, and indeed, we have that.

disp([Q0;Q2])

   1.0e+02 *
7.633444758094897
7.633444730985692


The default RelTol is 1e-6, so let's try making that smaller, instead of AbsTol, to get more accuracy.

Q3 = integral2(F,0,1,0,pi/4,'RelTol',1e-12)
intrelerr3 = abs(Q3 - Q0)/abs(Q0)

Q3 =
7.633444758094937e+02
intrelerr3 =
5.212639177138190e-15


Now since the relative error is between 1e-14 and 1e-15, we should have about 14 significant digits correct, and we do

disp([Q0;Q3])

   1.0e+02 *
7.633444758094897
7.633444758094937


Usually RelTol is the more important tolerance to manipulate, and the more useful one, because it indirectly controls the number of correct significant digits. But there is always a catch. Suppose the correct answer is zero. Then relative error is undefined! That's why integral2 has to use the least restrictive of AbsTol and RelTol. The actual error condition it tries to meet is

     ErrEst <= max(AbsTol,RelTol*abs(Q))

where ErrEst is an approximate upper bound on the absolute error.

It turns out that the ratio R = AbsTol/RelTol is a cutoff. When abs(Q) < R, AbsTol is the tolerance that matters. When abs(Q) > R, RelTol is the one that matters. So it is good practice to set both tolerances to something reasonable by asking yourself two questions:

1. If the true solution is small in magnitude, how many decimal places do I care about? The answer to that question tells you what AbsTol should be. If the answer is 10 places, then make AbsTol 1e-10 or perhaps 5e-11, so that the last digit should be correct when rounding to that point.
2. Assuming the answer is not too small, how many correct significant digits do I need? The answer to that question tells you what RelTol should be. Here again, if you would like 10 significant digits, make RelTol 1e-10 or perhaps 5e-11.

So when you choose your own tolerances, the calls should look something like

Q4 = integral2(F,0,1,0,pi/4,'AbsTol',1e-12,'RelTol',1e-10)

Q4 =
7.633444758095101e+02


Have You Been Successful with the Newer Integration Routines?

In this post, we've only covered the topic of setting tolerances to insure a reasonable estimate for your integrals. We plan to talk more about which method to choose, order of integration (e.g., x then y), and dealing with singularities and discontinuities.

In the meantime, do you have any experiences to share about performing your own numerical integrations with MATLAB? Let us know here.

Get the MATLAB code

Published with MATLAB® R2013b

Note

John D'Errico replied on : 1 of 4

Hi Loren,

I just learned about these tools from this blog, so thanks for the updates. Sometimes we miss the new features in each release, so these are a great way to learn. A few thoughts…

I’d like to see more output information, for example, how many function evals were taken? A nice idea is to return an exit flag, much as the optimization tools do. That flag might indicate whether the result meets the stopping tolerances, or if some problem was identified. So if an exitflag output is returned, then no warning messages need be generated, allowing the user to gather that information via an output argument.

John

Loren Shure replied on : 2 of 4

Thanks for the suggestions, John. Will enter them for the team.

–Loren

Mike Hosea replied on : 3 of 4

Hi, John. Thanks for the suggestions. There’s a little bit of “keep it simple” in the interface design of the new integral functions, and I agree that some additional features for programming uses would be helpful. In particular, a return flag or some data structure with information about how the integration went is sometimes what one needs. Conversely, when programming, warnings are more likely to be a nuisance that you have to deal with than something you can use.

As an aside, the number of function evaluations in particular is a traditional metric, but here it is difficult to interpret or use. This is because integral2 houses more than one method, and both methods are highly vectorized.

Ben Petschel replied on : 4 of 4

John, you might like to try the evaluation monitor utility PEEK (FEX #44120). For example:

F = @(x,y)exp(10*x).*tan(y)
Fpeek = @(x,y)peek(F(x,y),x,y)
Q1 = integral2(Fpeek,0,1,0,pi/4)
[Fe,Xe,Ye] = peek()

shows that integral2 used 6 function calls on mesh grids of size 14×14, 2×3, 14×14 etc. in R2013b, and you can plot the mesh grids coloured according to their iteration number with:

figure;view(3);colorbar;hold all;
cellfun(@surf,Xe,Ye,Fe,arrayfun(@(k)repmat(k,size(Fe{k})),1:numel(Fe),’Unif’,0)’);

Ben