Error propagation in adaptive Simpson algorithm

[Dieter]

(07-30-2018 06:42 PM)Claudio L. Wrote:  The first thing I notice is that this algorithm halves the requested tolerance every time it splits. Basically it's saying Error=Error_Left+Error_Right, therefore they want Error_Left and Error_Right to be less than half the requested tolerance to make sure it meets the requirement.
This seems to me quite conservative, once you are doing hundreds or thousands of steps, shouldn't the total error be proportional to the sqrt(N)*AvgError with N being the number of intervals considered?

Maybe a few general thoughts on this kind of algorithm may help.

First of all, adaptive Simpson methods have already been discussed here. Take a look at this thread and at this program.

The idea basically is to estimate the error of Simpson approximations with n and 2n intervals. This yields an even more accurate result with the value (16·S_2n – S_n)/15. That's where the 15 comes from.

For more details on this take a look at the German Wikipedia article where this is explained. The error term E^(N)(f) with the 16/15 factor is discussed earlier in that article.

As far as the threshold for quitting the iteration is concerned: this is a general problem with any method that produces successive approximations of the true result. Since the latter is unknown all you can do is compare the last two approximations, i.e. the last and the previous one. If these two agree within the desired accuracy, the iteration stops and the result is returned. But – if the final and the previous value agree in the desired, say, 8 places, this means that already the previous value was sufficiently accurate! And since the last approximation is more accurate than this, it may be correct to e.g. 11 places. So its true accuracy is greater than requred. That's what you observe here.

(07-30-2018 06:42 PM)Claudio L. Wrote:  For example (from another thread): Integral of exp(-x) between 0 and 500.
If I request 1e-8 tolerance, the result is:
1.00000000000704832....
This case has 2 more digits than requested.

Yes, that's the ex-post accuracy when comparing the result with the true value. But the algorithm doesn't know the true result. All it can do is compare the final and the previous approximation.

Here and there the exit condition can be improved. Assume you know for sure (!) that the iteration converges quadratically and the difference between the last two approximations was 1E–2, 1E–4 and 1E–8. If the next approximation can be expected to change the result only in the 16th digit this means that the current (!) approximation already is this precise. So if you need 12 places you may stop now although the last two approximation only agreed in 8 places. Because one more iteration will only change the result in the 16th place.

BTW, integrating exp(–x) up to 500 doesn't make sense. At x=20 the function value has dropped to about 2E–9 so anything beyond this will not affect the first 8 decimals.

Dieter