Skip to content

Ahrens-Beylkin 552-point grid is not a valid order-39 rule #108

Description

@susilehtola

The tabulated 552-point Ahrens-Beylkin grid
(include/integratorxx/quadratures/s2/ahrens_beylkin/ahrens_beylkin_552.hpp)
is not a valid quadrature rule. It is advertised as algebraic order 39, but:

  • its weights sum to 12.104473167902746, i.e. 0.9632 * 4pi, so it does
    not even integrate a constant correctly;
  • after renormalising to 4pi it integrates monomials up to degree 4
    exactly, but is still wrong at degree 6 (7.5e-03 relative) and degree 8
    (1.4e-02 relative).

So it is not merely mis-normalised -- renormalising would produce a grid
that looks correct at low order and is silently wrong at high order, which
is worse than the current state.

Every other size in the family is exact to machine precision:

  AB 432   sum(w)=12.566370614359144   sum/4pi=1.000000000
  AB 492   sum(w)=12.566370614359053   sum/4pi=1.000000000
  AB 552   sum(w)=12.104473167902746   sum/4pi=0.963243369   <---
  AB 612   sum(w)=12.566370614359178   sum/4pi=1.000000000
  AB 672   sum(w)=12.566370614359140   sum/4pi=1.000000000

I swept all 230 grid sizes advertised by the npts_by_algebraic_order
tables across all four angular families, checking exact integration of even
monomials up to each grid's own claimed algebraic order (capped at degree 8
for run time). This is the only failure:

Lebedev        : 32 sizes checked, 0 bad
Delley         : 18 sizes checked, 0 bad
AhrensBeylkin  : 55 sizes checked, 1 bad   <-- 552
Womersley      : 125 sizes checked, 0 bad

Reproduction:

#include <integratorxx/quadratures/s2/ahrens_beylkin.hpp>
#include <numeric>
#include <cstdio>
int main() {
  IntegratorXX::AhrensBeylkin<double> q(552);
  auto& w = q.weights();
  printf("sum(w) = %.15f  (4pi = %.15f)\n",
         std::accumulate(w.begin(), w.end(), 0.0), 4*M_PI);
}

All 552 points are distinct and of unit norm, and the weights are
individually of plausible magnitude, so the table looks like it lost or
duplicated a subset of the rule rather than being corrupted wholesale. I
have not been able to determine the exact provenance of the error.

Interim action taken. PR #101 withdraws the size from the
dispatch and from the algebraic order tables, so requesting 552 points now
throws rather than returning a 3.7%-wrong quadrature, and
next_algebraic_order rounds order 39 up to 41. The data file is left in
place.

What is needed. Regenerating the table from the original Ahrens-Beylkin
data (Proc. R. Soc. A 465, 3103 (2009)) and re-enabling the size. The
polynomial-exactness test added in the same PR will confirm a correct
replacement.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions