Integrate an M-Dimensional Polynomial#

import numpy as np
import minterpy as mp

import warnings
warnings.filterwarnings('ignore')

A definite integration operation may be carried out for fully-specified Minterpy (multivariate) polynomials. This guide explains how to carry out a definite integration on a three-dimensional polynomial interpolant in different bases.

Motivating example#

Consider the following three-dimensional function:

\[ f(\boldsymbol{x}) = \prod_{m = 1}^3 (m \cos{(m x_m)}),\;\; x \in [-1, 1]. \]

Define the function in Python as follows:

def fun(xx):
    m = np.arange(1, xx.shape[1] + 1)
    return np.prod(m * np.cos(m * xx), axis=1)

Polynomial interpolation#

In this guide, we are going to create a polynomial interpolant in Minterpy from scratch in four steps:

  1. Define the multi-index set

  2. Create the interpolation grid (of unisolvent nodes)

  3. Evaluate the function on the grid

  4. Create a polynomial interpolant in Lagrange basis

A polynomial interpolant of a given degree may be created using the Lagrange basis. First, create the multi-index set:

m = 3
mi = mp.MultiIndexSet.from_degree(spatial_dimension=m, poly_degree=20, lp_degree=1.0)

Because the function is highly irregular, we use a high polynomial degree to interpolate the function. Then, create the interpolation grid given multi-index set:

grd = mp.Grid(multi_index=mi)

The grid contains unisolvent nodes on which the function should be evaluated as the coefficients of a polynomial in the Lagrange basis:

lag_coeffs = fun(grd.unisolvent_nodes)

Finally, a polynomial interpolant in Lagrange basis is created from the multi-index set and the set of coefficients:

lag_poly = mp.LagrangePolynomial(multi_index=mi, coeffs=lag_coeffs)

Note

In Minterpy, polynomials in the Lagrange basis cannot be directly evaluated. To evaluate an interpolating polynomial, the Newton basis is recommended.

Integration over the domain \([-1, 1]^3\)#

The method integrate_over() integrates the polynomial over the default domain of \([-1, 1]^3\):

int_value = lag_poly.integrate_over()
int_value

0.8638208594281641

The analytical value of the integral of the above function is available:

\[ I[f] \equiv \int_{-1}^{1} f(\boldsymbol{x})\; d\boldsymbol{x} = 2^3 \prod_{m = 1}^3 \sin{(m)}. \]
ref_value = 2**m * np.prod(np.sin(np.arange(1, m + 1)))

With polynomial degree of \(20\) and \(l_p\)-degree of \(1.0\), we obtain the following relative error of the integral:

np.abs(ref_value - int_value) / ref_value

np.float64(6.013412886351958e-12)

Integration over specified bounds#

The bounds of the integration may be specified. For instance, to integrate the polynomial over \([-1, 0]^3\), specify the bounds for each dimension as follows:

lag_poly.integrate_over([[-1, 0], [-1, 0], [-1, 0]])

0.10797760744251064

and similarly, to integrate the polynomial over \([0, 1]^3\):

lag_poly.integrate_over([[0, 1], [0, 1], [0, 1]])

0.10797760740695861

Note

If specified, then the bounds for each dimension must be provided and be provided in order (i.e., the first bounds are for the first dimension, etc.).

Integration in different bases#

The integration may also be carried out on the polynomial in different bases. The same method integrate_over() is used. Below is the integration in the Newton basis:

nwt_poly = mp.LagrangeToNewton(lag_poly)()
nwt_poly.integrate_over()

0.8638208594282044

And in the canonical basis:

can_poly = mp.LagrangeToCanonical(lag_poly)()
can_poly.integrate_over()

0.8638208564704958

Warning

The integration of polynomials in the canonical basis having a high polynomial degree is not recommended as it may suffer from a severe instability.