Error propagation in adaptive Simpson algorithm

[Dieter]

This is a reply to Claudio L's post #4 in this thread. Since the "Quote" button currently doesn't seem to work I have to do it this way:

Quote:Thanks, now I know where the 15 comes from, yet we are talking very different algorithms from the one you provided links to. I'm not even sure why they still call the adaptive one "Simpson", since you are not using the summation that is the hallmark of the Simpson method, but the error analysis still applies.

If you refer to the HP-41 program I linked to, or the VBA version (post #8) in the other thread: they both use the classic Simpson method. The sums with weight 4 and 2 are summed separately and a Simpson approximation is calculated from these. Then the number of intervals is doubled and the two sums are added together to get the new weight-2-sum. This way only the new points have to be calculated which sum up to the new weight-4-sum.

The additional "non-Simpson"-step then is the extrapolation (16·S_new – S_old)/15 which gives an improved, more accurate value for the integral. When two successive values of this improved approximation match with the desired accuracy the iteration exits.

Take a look at the VBA code: simp_old and simp_new are the regular Simpson approximations for n/2 and n intervals, whereas improved_old and improved_new are the extrapolated improved approximations.

Both the HP-41 program and the VBA version in the 2016 thread only use the last two Simpson approximations to calculate a new, improved integral. In his recent answer Albert Chan suggested using all previous Simpson approximations, similar to the Romberg method. This is another worthwile option, but I didn't want to work with arrays, so only the last two approximations were used.

Quote:You do have a good point, perhaps the overshoot in error is because halving the step increases the accuracy a lot more than what we needed.

The essential point here is that we cannot and do not know the error of the calculated approximation. We can't say for sure if it's 7, 8 or 9 digits that match the true result. All we got is an upper bound, and even this is only reliable if we assume that the iteration converges to the true result. In this case the upper error bound is the difference between the last two approximations. If they match in 6 decimals this can mean that the result is accurate to 6, 7, 8 or even more digits. But we can't say anything more than this.

To make this more clear, take a look at the Simpson approximations for the integral of f(x)=1/x from 1 to 2:

Code:

   n      Standard Simpson    Error
 -------------------------------------
   2      0,694444444444      1,3 E-03
   4      0,693253968254      1,1 E-04
   8      0,693154530655      7,4 E-06
  16      0,693147652819      4,7 E-07
  32      0,693147210290      3,0 E-08
  64      0,693147182421      1,9 E-09
 128      0,693147180676      1,2 E-10

Assume you want an accuracy of six decimals. The Error column shows that n=16 already met that target. But we can't know this during the calculation unless we also know the true value of the integral! The Simpson approximations for n=16 matches the previous one (n=8) in (nearly) 5 decimals, so after calculating the n=16 iapproximation there are only two things we can say:

1. The n=8 approximation agrees with the (most probably) more accurate n=16 approximation when rounded to 5 decimals (0,69315). So already the n=8 approximation probably had 5 valid decimals. Which is what we know after having calculated the n=16 approximation.

2. The n=16 approximaion most probably is more accurate than this, but we can't say by how much. It may have 6 valid decimals or more.

Dieter