Question:
Do you know some numerical integration rule (e.g. like Simpson rule, Gauss-Legendre etc ) which provides some rough estimate of integral (e.g. error $\epsilon < 10^{-2}-10^{-3}$) with minimal number of points for any reasonably smooth function (i.e. with derivatives bound by some number $\delta$ ).
Do you know some general algorithm how to derive optimal integration rule (i.e. position of sampling points and weights) for given domain (e.g. rectangle, triangle, sphere ... ) and expected distribution of the functions (e.g. assuming most of function variation will be in the center of integration domains)
Motivation
Most integration rules like Newton-Cotes, Gauss-Legendre or Clenshaw–Curtis are developed to be exact for polynomial of given degree. Or to provide high asymptotic order of convergence for small integration step (large number of sampling points).
But in practice in many numerical applications you very often need something else. You are fine with rather low precision ($10^{-2} - 10^{-3}$) but you want to achieve it with as little sampling points as possible (especially in 3D and more-dimensional integrals the number of integration point become quickly unfeasible).
Also, you have no reason to expect the functions are polynominals. The very reason why to use numerical integration rule is that you have have some arbitrary function for which you cannot calculate the integral analytically.
More rigorous definition of the problem
Lets assume space of all possible analytic functions ${f_i(\vec x)}$ on $N$-dimensional interval $(0 ... 1)^N$ which has finite 1st and 2nd derivative everywhere. (In practice $N\leq 3$, i.e. 1D,2D,3D are most interesting). The value of 1st and 2nd derivative is also bound by some value $\delta_1$ and $\delta_2$.
We want to find sampling points (nodes) $\vec x_j$ and associated weights $w_j$ such that for any function ${f_i(\vec x)}$ the
$S_n = \sum_j^n w_i f_i(\vec x_j)$
approximates the $N$-diminsional integral $I_i = \int_{{\vec x} \in (0..1)^N} f_i(\vec x) $ with accuracy better than choosen $\epsilon$ where $\epsilon$ is in th order $10^{-2}$-$10^{-3}$.
In short
$| \sum_j^n w_i f_i(\vec x_j) - \int_{{\vec x} \in (0..1)^N} f_i(\vec x) |<\epsilon$
Sketch of the solution
The space of such smooth functions ${f_i(\vec x)}$ can be well approximated e.g. by finite Fourier series of rather small number of terms. The bound of 1st and 2nd derivatives $\delta_1$ and $\delta_2$ guarantie rapid convergence of the series with error $<\epsilon/2$. This means we can transform the problem into relatively small finite space of lets say 10^6 orthogonal functions. Due to symmetry it should be possible to reduce it just to 1 dimensions and half of the period, i.e. to ~50 Fourier components. Than it should be possible to find optimal sapling points and weights for fourier series of 50 terms by some algorithm or even brute force.
Is there anything better than Midpoint rule ?
From numerical experiments I did up to now it seems that simple mid-point rule (i.e. equidistant points offset by half of spacing step all with same weights) is the solution, at least for 1D case.
This would make lot of sense:
- Since we don't know anything about the function other than bounds on derivatives, the peaks of the function (where most of density is concentrated) can be anywhere. Therefore is only reasonable to distribute the sampling points as equally as possible, all with equal weights. This gives you best chance that by minimum density of sampling point at least one sampling point hit all the important peaks of the function. (see e.g. discussion of Simpson's_rules_in_the_case_of_narrow_peaks )
In more-dimensional space it is less trivial
However even if conclude that all weights the same is optimal is optimal and nodes should be on equally spaced, geometry of the lattice can still play a role. I suppose that Face-Centered-Cubic (FCC) o Body-Centered-Cubic will be better than Simple Cubic (i.e. regular rectangular grid).
Boundary conditions
I'm entirely not sure how to approach the problem when it is defined on different domain than N-dimensional interval. E.g. on sphere, tetrahedron (simplex) etc. In case of sphere I would also expect that most of function variation is in the center of the sphere as is the case in most practical applications (e.g. center of an atom, or of galaxy) which should be reflected in optimal distribution of the points.
Clenshaw–Curtis and Gauss-Legendre are quite inefficient in that regard since they put most sampling points close to the end of integration domain, while in practical numerical applications it is only natural to choose integration domain such that region of interest is in the center and almost nothing happens close to the boundary.