[EM] High Order Numerical Cubature (Clean and Simple)

[Forest Simmons]

Forget every article I sent you about numerical quadrature/cubature ...
last night,  in the night, vast simplification came to me ... it is much
easier to generalize the midpoint rule than the trapezoid rule for cubature
in arbitrarily many dimensions.

There is historical precedence for this idea ... ... the
*Gragg–Bulirsch–Stoer
(GBS) * numerical solution method for ODE's uses a combination of the
midpoint rule and Richardson extrapolation for high order accuracy. We are
going to generalize this to multivariate quadrature or "cubature."

Given a simplex (the convex hull of a finite set of points having affine
independence) we form a tree (more of a dendrite, really) whose root node
is the average of the vertices of the simplex (i.e. the "barycenter" of the
simplex), and whose descendants become denser in the interior of the
simplex with each successive generation.

This "space filling" has to be done with sufficient fractal self
symmetrical regularity to ensure that our method can be highly accurate and
efficient. Don't worry ... it turned out to be way simpler than I could
have possibly hoped for!

As we said, the first node of the tree is the baryceenter of the simplex.
The immediate children of that node are the points precisely half way from
the baryceenter to the respective vertices of the simplex. [The simplex is
the convex hull of these vertices.]

The respective arrows (or "vectors") from the barycenter to the vertices
serve as the only directions needed for finding all of the children of
every node.

The only question is "how far?"

Answer: the distance from a parent to a child is half as far as it was in
the previous generation ... i.e. half as far as from grandparent to parent.

Now, let g(k) represent the set of nodes in the kth generation down from
the root node. Then g(1) and g(2), respectively, represent the children and
grandchildren of the barycenter, while g(0) consists of only the barycenter
itself.

Also let G(k) be the union of all of the generations from g(0) to g(k),
inclusive.

Now let f be a real valued function defined at every node of our tree.

And let A(k) be the average of all of the f values of all of the nodes of
the set G(k), i.e. from the barycenter value down to all of the values of
the nodes of the k-th generation.

If f is a reasonable function, then as k increases, A(k) approaches the
integral of f over the simplex, divided by its "capacity" (length, area,
volume, etc ... depending on the dimension of the simplex).

Let's denote this limiting value of A(k) by the Greek letter Alpha.

In particular, when f is a density function, Alpha is the average density
of the simplex, because the integral of density is mass, and mass divided
by capacity is average density ... in units of g/(cm)^d or kg/m^d, where d
is the dimension of the simplex.

Similarly, if f represents temperature or pressure, Alpha will be the
average temperature or pressure experienced by the simplex.

If you need to know the value of the integral of f over the simplex, just
multiply its capacity by Alpha.

To get the capacity of a simplex, let M  be a matrix whose rows are the
respective displacement vectors from one fixed vertex v0 to the others ...
v1-v0, v2-v0, etc.

The capacity of a d-dimensional simplex is the square root of the
determinant of the matrix product of M and its transpose, divided by
d-factorial :

Sqrt(det MM')/d!,

where M' is the transpose of M.

That's everything ... except the Richardson extrapolation to increase the
order of convergence from second order to order 2n, where n depends on the
number of Richardson transformations of the sequence [see penultimate
paragraph below]:

We start with the sequence
A(0), A(1), .... A(j), ...

Then transform the A sequence to a B sequence, where

B(0)=(4A(1)-A(0))/3 ...
B(j)=(4A(j+1)-A(j))/3 ...

Then transform the B's to a sequence of C's:

C(j)=(16B(j+1)-B(j))/15 ...

.......

D(j)=(4^3*C(j+1)-C(j))/(4^3-1)

etc.

The suggested order is ..

A0), A(1),B(0), A(2), B(1), C(0), A(3),B(2), C(1), D(0), etc....

At each step go as far down the alphabet as possible, given the previous
steps.

When you reach a new letter, i.e. the zero_th term in a new sequence,
compare it with the most recent value of the preceding letter to see if the
value is still progressing, or if it has reached a point of diminishing
returns.

Test it on polynomials: by the time the extrapolation reaches the n_th
letter of the alphabet, any polynomial function f of degree 2n or less
should be averaged exactly (up to roundoff error).

If your polynomials have integer coefficients, and the coordinates of the
vertices are integers or precise fractions, then there is no excuse for
roundoff error ... the answers must have infinite precision! (unless you
cannot afford it)

Have Fun!

Forest
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.electorama.com/pipermail/election-methods-electorama.com/attachments/20220203/7f89437c/attachment.html>