Numerical integration methods

[Wes Loewer]

(07-19-2022 05:16 PM)parisse Wrote:  int uses an adaptive method, as described here in French (Hairer: https://www.unige.ch/~hairer/poly/chap1.pdf).

A huge thanks to Parisse for adding this adaptive method to the Prime. The Romberg method made sense back in the day when calculators had very limited memory, but with ample memory now available, adaptive methods are more appropriate. Certain integrands with cusps or asymptotes could slow Romberg to a crawl. The Prime usually makes short work of these.

Since many of us don't read French (myself included), allow me to summarize the article. (Some years ago with the help of Google Translate, I read through this article and I think I've got the gist of it. Corrections welcome.)

The Prime's adaptive method is a variation on Gaussian quadrature. Many are familiar with the more commonly used Gauss-Kronrod Method, so let me use that as a starting point for comparison. Calculators typically use a 15-node Gauss-Kronrod implementation that starts out with a 7-node Gaussian quadrature (exact for polynomials up to degree 13). Then, those 7 nodes are reused along with 8 additional Kronrod nodes using different weights for all 15 nodes to get a more accurate result (exact up to degree 22). The difference between these two calculations is used to estimate the error.

The Prime's method described in the Hairer document starts off with a 15-node Gaussian quadrature (exact up to degree 29). The integral is recalculated a 2nd time using 14 of the original 15 Gaussian nodes with different weights (exact up to degree 13). The integral is then recalculated a 3rd time using 6 of the original Gaussian nodes, again with different weights (exact up to degree 5). These three results are used to estimate the error.

In both methods, the function is evaluated a total of 15 times. If the estimated error is small enough, the more accurate calculation is used. If the error is too large, then the region is cut in half and each half is recursively calculated.

(Note: It's possible that I am mistranslating "ordre". What the paper called "ordre 30" seems to be "degree 29" (ie, 30 terms), but I could be wrong about this. Francophones, please correct me if I am mistaken.)

Since the 15-node Gauss-Kronrod is exact up to degree 22 and the Prime's method is exact up to degree 29, it would seem that the Prime's method is a bit better with only a small added overhead. A slightly more accurate calculation per iteration could occasionally reduce the need for further recursion. It would be interesting to do some comparisons of these two methods with some well-behaved functions as well as more "temperamental" ones.

On a final note, I am curious why Hairer chose only 6 nodes for the 3rd calculation instead of 7. Any insight on that would be welcomed. [Edit: Perhaps it is important to have the 3rd calculation's 6 nodes be a subset of the 2nd calculation's 14 nodes. This would not be the case if 7 nodes were used for the 3rd calculation. It is not clear to me why this subset would be necessary or even desirable.]