Polynomial definite integration

The functionality to carry out definite integration of a Minterpy polynomial may be added to the core.

The integration is based on the given polynomial and it returns a value. The main ingredients are the multi-index set, the generating points, and the bounds.

Below are a couple of sanity checks involving simple polynomials on the Lagrange basis integrated (by integrating the Lagrange basis represented in Newton form).

One-dimensional Gauss-Legendre quadrature

As a sanity check:

• if the roots of the N-th Legendre polynomial are given as the generating points of the Minterpy polynomial,
• then the computed Lagrange-Newton quadrature weights should be equal to the quadrature weights of the Gauss-Legendre quadrature:

\sum_{j = 1}^{N} c_{i}^{(j)} \int_{-1}^{1} N_{j} (x) \, dx = \frac{2}{\left( 1 - \left( x^{(i)} \right)^2 \right) \left[ P'_N ( x^{(i)} ) \right]^2}, \; i = 1, \ldots, N

where P'_N is the derivative of the N-th Legendre polynomial.

And indeed they are the same within some precision:

N	x^{(i)}	\omega^{(i)} Analytical	\omega^{(i)} Lagrange-Newton
1	0.0	2.0	2.0
2	-\frac{1}{\sqrt{3}}	1.0	1.0
	\frac{1}{\sqrt{3}}	1.0	1.0
3	-\sqrt{\frac{3}{5}}	\frac{5}{9}	0.55555556
	0.0	\frac{8}{9}	0.88888889
	\sqrt{\frac{3}{5}}	\frac{5}{9}	0.55555556
4	-\sqrt{\frac{3}{7} - \frac{2}{7}\sqrt{\frac{6}{5}}}	\frac{18 + \sqrt{30}}{36}	0.65214515
	-\sqrt{\frac{3}{7} + \frac{2}{7}\sqrt{\frac{6}{5}}}	\frac{18 - \sqrt{30}}{36}	0.34785485
	\sqrt{\frac{3}{7} - \frac{2}{7}\sqrt{\frac{6}{5}}}	\frac{18 + \sqrt{30}}{36}	0.65214515
	\sqrt{\frac{3}{7} + \frac{2}{7}\sqrt{\frac{6}{5}}}	\frac{18 - \sqrt{30}}{36}	0.34785485

One-dimensional polynomial

An integration of a simple one-dimensional polynomial: \int_{-1}^1 10 x^{100} dx = \frac{20}{101}.

Shown above is the comparison with current minterpy integrator which is based on the canonical polynomial.

We can compare the result from different methods of integrating the Newton monomials, namely: scipy.integrate.quad(), scipy.integrate.quad_vec(), and Gauss-Legendre quadrature (because we know in advance Newton monomials are polynomials). Also shown is the result from directly carrying out Gauss-Legendre quadrature (ref_gauss_quad in the figure) on the polynomial functions of different degrees:

Both functions from SciPy become sensitive to the parameters (i.e., tolerance, number of iterations) for high-degree polynomials. Perhaps they can be made better with the proper choice of the parameters. The Gauss-Legendre quadrature approach to integrate the monomials works fine, but seemingly sub-optimal compared to the directly carrying out Gauss-Legendre quadrature on the function. This is because the generating points (i.e., interpolation nodes) chosen are Chebyshev points. If we were to use Legendre nodes to build the polynomial then integrate it, we will recover the convergence behavior:

Two-dimensional polynomial

An integration of a simple two-dimensional polynomial: \int_{-1}^1 \int_{-1}^1 (-2 x_1^4 + 4 x_2^5) dx_1 dx_2 = - \frac{8}{5}.

Shown above is the comparison with current minterpy integrator which is based on the canonical polynomial.

Three-dimensional polynomial

An integration of a simple two-dimensional polynomial: \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 (-2 x_1 + 4 x_1^3 + 3 x_1^5 x_3^3 + 10 x_3^7) dx_1 dx_2 dx_3= 0.

Shown above is the comparison with current minterpy integrator which is based on the canonical polynomial.

Ten-dimensional polynomial

Finally, the integration of a high-dimensional polynomial function with a simple structure:

\int_{0}^1 \frac{1}{\sqrt{\frac{M}{12}}} \left( \sum_{m = 1}^{10} x_i - \frac{M}{2} \right) \; d\boldsymbol{x} = 0.0

Here we can see that absolute error grows with the number of spatial dimensions, this is to be expected as the magnitude of the integral value scales up with the number of dimension (i.e, 2^M). The scaled error (relative to the volume of the hypercube), however, remains stable.

Proposal: High-level interface

We can imagine for a given Minterpy polynomial instance, a definite integration can be computed as follows:

>>> import minterpy as mp
>>> lag_poly = mp.LagrangePolynomial(mi, lag_coeffs)
>>> lag_poly.integrate()  # By default, integration is on the whole domain [-1, 1]^M
# A value
>>> lag_poly.integrate(bounds=[[-0.5, 0.5], [-0.5, 0.5],...])  # The bound within [-1, 1]^M may also be given
# A value

The method integrate() should be added to the abstract class MultivariatePolynomialSingleABC, while the particular implementation for each of the available concrete classes, namely LagrangePolynomial, NewtonPolynomial, and CanonicalPolynomial should be added to each (such that nwt_poly.integrate() and can_poly.integrate() are valid). Carrying out a definite integration from a given polynomial, therefore, is as routine as getting the differentiated polynomial from the given polynomial (via diff() and partial_diff()).

Relation with the minterpy-integration

There is an on-going project called minterpy-integration whose scope is different from the current proposal. The minterpy integrator is aimed for the integration (adaptive or otherwise) of general functions (polynomial or otherwise). The current proposal only deals with the definite integration of a given minterpy polynomial. In other words, the multivariate polynomial as described by the given basis, multi-index set, and generating points are given and taken for granted. The task is to simply integrate the multivariate polynomial on given bounds and return the value of the integral. It does not do any fancy error estimation and control to obtain the best integral value for the smallest number of function evaluations.