# 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.

## 评论

要发表评论，请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。