numerical integration algorithm details

[Wes Loewer]

A few years ago Bernard Parisse shared this link to an article that explains the numerical integration algorithm used by the Prime. The algorithm described works quite well. I'm curious about the mathematical reasons behind a couple of decisions made in the algorithm.

The method uses a standard 15-point Gaussian Quadrature and then a 14-point and a 6-point calculation to estimate the error. The 14 and 6 point calculations reuse nodes from the 15-point method but with different weights. I've tried to show this visually below using X to represent the nodes used. (df=degrees of freedom)

Code:

XXXXXXXXXXXXXXX - 15 nodes, df=30, exact for degree 29 polynomial
XXXXXXX XXXXXXX - 14 nodes, df=14, exact for degree 13 polynomial
 X X X   X X X  -  6 nodes, df=6,  exact for degree 5 polynomial

My question relates to why 6 nodes were chosen. As the least accurate part of the calculation it seems like you'd want to use more nodes. Using 7 nodes seems more natural. Why skip the middle node when you've already skipped the two adjacent nodes? The original 15 nodes are already distributed non-uniformly towards the endpoints.

Code:

 X X X X X X X  -  7 nodes, df=7,  exact for degree 6 polynomial

Or why not use even more, say 13 nodes

Code:

XXXXXX X XXXXXX - 13 nodes, df=13, exact for degree 12 polynomial

So my question boils down to, "Why 6?"


My second question pertains to the error calculation. Using ERR1 as the difference between the 15 and 14 point calculation, and ERR2 as the difference between the 15 and 6 point calculation, the error is estimated as

error = ERR1*(ERR1/ERR2)^2.

This seems to work well, but I don't see how it was derived. Why squared? Or perhaps it was simply found by experimenting and they found that this worked reasonable well.

Thanks.