Adaptive Simpson and Romberg integration methods

[robve]

(03-25-2021 03:09 PM)Albert Chan Wrote:  Trapezoids seems to be twice as good as mid-points (in the opposite direction).
When we extrapolate with Romberg, trapezoids are *much* better

Yes and no, as I think you are pointing out. When you use them in isolation then generally yes and when you use them in succession to sum up a quadrature. No when you used them as building blocks for the "extended" open and closed formulas used in Romberg. The best choices for Romberg are the extended trapezoidal rule for closed intervals

$$
\int_{x1}^{x_N} f(x)\,dx = h\left[ \frac12 f_1 + f_2 + f_3 + \cdots + f_{N-1} + \frac12 f_N \right] + O\left(\frac{(b-a)^3f''}{N^2}\right)
$$

and the extended midpoint rule for open intervals

$$
\int_{x1}^{x_N} f(x)\,dx = h\left[ f_{3/2} + f_{5/2} + f_{7/2} + \cdots + f_{N-3/2} + f_{N-1/2} \right] + O\left(\frac{1}{N^2}\right)
$$

So there is definitely not a "twice as good" comparison when we consider Romberg integration. Romberg combines the previous interpolated result with the next using extrapolation to produce even better results by cancelling out leading error terms by the combination, in a more general way than who Simpson's rule cancelled out leading error terms, which is by:

$$ S = \frac43 S_{2N} - \frac13 S_N $$

to get \( O(1/N^4) \) instead of \( O(1/N^2) \). But Simpson's is a special case, while Romberg is a general. That is the whole point of using Romberg. In fact, for \( k = 2 \) refinements, Romberg is the same as Simpson's rule and much better for larger \( k \).

Romberg's \( k \) successive refinements by subdividing the interval has an error series up to but not including \( O(1/N^{2k}) \) for trapezoidal. With the extended midpoint rule, there is a problem with doubling the steps which is not possible to preserve the benefits of the method. Steps are tripled. This does \( \sqrt{3} \) additional work that may be unnecessary instead of \( \sqrt{2} \) for extended trapezoidal rule, but we are now able to reuse the previous results to obtain an error series up to but not including \( O(1/N^{3k}) \) for extended midpoint. Remember that each step triples the number of points (extended midpoint) rather than doubles (extended trapezoidal), so a lower max \( k \) is applicable to midpoint to get the same number of points as trapezoidal by a ratio of 2:3, which to some extent can be experimentally confirmed in the table and in other runs I've tried. The ratio 2:3 is of course idealized. In practice these subtle differences in the evaluated points can throw this off quite easily*)

There is a bit more theory behind this, which would make this post a lot longer. An excellent summary can be found in Numerical Recipes Ch. Integration of Functions.

Sorry if it was not clear what "trapezoidal" and "midpoint" meant in the original posts.

Edit: *) and there is of course the difference in the little bit of unnecessary work that is thrown away. Sometimes the difference can be big, however such as:

$$
\begin{array}{cccccc}
& f(x) & \int_0^1 f(x)\,dx & \texttt{qromb} & \texttt{qromo} & \texttt{qasi} \\
\hline
19 & x^7-3x^6+5x^5-7x^4+7x^3-5x^2+3x-1 & -0.28690476190 & 17^= & 81^= & 205 \\
\end{array}
$$

Edit 2: I briefly read your post in the other thread on trapezoidal and midpoint. I could be mistaken, but I don't think this is correct to do to get (extended) midpoints:
for i = 1, n, 2 do m = m + f(a + i*h) end
m = m*h -- midpoint area / 2
t = t*0.5 + m -- trapezoid area

- Rob