Double Integration in MATLAB – Understanding Tolerances
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
Adjust the Tolerance(s)
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:
- 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.
- 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.
- Category:
- Integration,
- Numerical Accuracy
