<div dir="auto">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.<div dir="auto"><br></div><div dir="auto">There is historical precedence for this idea ... ... the    <b style="margin:0px;padding:0px;border:0px;line-height:inherit;font-family:-apple-system,blinkmacsystemfont,"segoe ui",roboto,lato,helvetica,arial,sans-serif;font-size:16px;vertical-align:baseline;background:none rgb(255,255,255);color:rgb(32,33,34)">Gragg–Bulirsch–Stoer (GBS) </b> 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."</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">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!</div><div dir="auto"><br></div><div dir="auto">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.]</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">The only question is "how far?"</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">Also let G(k) be the union of all of the generations from g(0) to g(k), inclusive.</div><div dir="auto"><br></div><div dir="auto">Now let f be a real valued function defined at every node of our tree.</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">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).</div><div dir="auto"><br></div><div dir="auto">Let's denote this limiting value of A(k) by the Greek letter Alpha.</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">Similarly, if f represents temperature or pressure, Alpha will be the average temperature or pressure experienced by the simplex.</div><div dir="auto"><br></div><div dir="auto">If you need to know the value of the integral of f over the simplex, just multiply its capacity by Alpha.</div><div dir="auto"><br></div><div dir="auto">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 ...</div><div dir="auto">v1-v0, v2-v0, etc. </div><div dir="auto"><br></div><div dir="auto">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 :</div><div dir="auto"><br></div><div dir="auto">Sqrt(det MM')/d!,</div><div dir="auto"><br></div><div dir="auto">where M' is the transpose of M.</div><div dir="auto"><br></div><div dir="auto">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]:</div><div dir="auto"><br></div><div dir="auto">We start with the sequence </div><div dir="auto">A(0), A(1), .... A(j), ...</div><div dir="auto"><br></div><div dir="auto">Then transform the A sequence to a B sequence, where </div><div dir="auto"><br></div><div dir="auto">B(0)=(4A(1)-A(0))/3 ...</div><div dir="auto">B(j)=(4A(j+1)-A(j))/3 ...</div><div dir="auto"><br></div><div dir="auto">Then transform the B's to a sequence of C's:</div><div dir="auto"><br></div><div dir="auto">C(j)=(16B(j+1)-B(j))/15 ...</div><div dir="auto"><br></div><div dir="auto">.......</div><div dir="auto"><br></div><div dir="auto">D(j)=(4^3*C(j+1)-C(j))/(4^3-1)</div><div dir="auto"><br></div><div dir="auto">etc.</div><div dir="auto"><br></div><div dir="auto">The suggested order is ..</div><div dir="auto"><br></div><div dir="auto">A0), A(1),B(0), A(2), B(1), C(0), A(3),B(2), C(1), D(0), etc....</div><div dir="auto"><br></div><div dir="auto">At each step go as far down the alphabet as possible, given the previous steps.</div><div dir="auto"><br></div><div dir="auto">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.</div><div dir="auto"><br></div><div dir="auto">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). </div><div dir="auto"><br></div><div dir="auto">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)</div><div dir="auto"><br></div><div dir="auto">Have Fun!</div><div dir="auto"><br></div><div dir="auto">Forest</div><div dir="auto"><br></div><div dir="auto"><br></div></div>