In this tutorial we describe how to use Equadratures capabilities to numerically compute integrals using efficient quadratures.
import numpy as np
import matplotlib.pyplot as plt
import equadratures as eq
import plotly.graph_objects as go
import plotly.express as px
import pandas as pd
order = 100
parameter = eq.Parameter(lower=-1, upper=1, order=order, distribution='uniform')
We will first consider univariate polynomials, evaluating using Gauss-Legendre Quadrature.
Above, we being by defining a uniform paramter from -1 to 1 and of order 100, meaning that we are only considering the first 100 orthogonal polynomials in our quadrature computations.
$$ \int_{-1}^{1}f(x)dx $$
Utilizing Gauss-Legendre Quadrature (since we are using a uniform weight distribution), we can evaluate the above integral by evaluating the function at certain points, that correspond to the roots of the Legendre polynomial, and dot those evaluations with weights associated with those points.
basis = eq.Basis('univariate')
poly = eq.Poly(parameter, basis, method='numerical-integration')
points, weights = poly.get_points_and_weights()
Equadratures allows simple retrieval of these weights and points, allowing us to evaluate computationally efficient integrals.
$$ \int_{-1}^{1}f(x)dx \approx \sum_{n=1}^{i}f(points_i)*weights_i $$
def f(x):
return np.cos(x**2);
Consider the function: $$ f(x) = cos(x^2) $$
Utilizing the weights and points we previously collected, evaluation is easily done through the following.
eval = float( 2 * np.dot(weights , eq.evaluate_model(points, f) ) )
print(eval)
1.8090484758005432
$$ \int_{-1}^{1}cos(x^2)dx \approx \sum_{n=1}^{i}f(points_i)*weights_i = 1.809048475800536 $$
This approach can be altered to help compute a wide variety of integrals. Including switching away from Gauss-Legendre, by specifying a non uniform distribution.
We will now quickly showcase how to use equadratures to compute integrals with non uniform weight distributions such as a Gaussian distribution. Consider the following function:
$$ L(x) = \int_{-\infty}^{\infty}sin(x)\rho(x) \approx \sum_{n=1}^{i}f(points_i)*weights_i $$ where... $$ \rho(x) = \mathcal{N}(\mu, \sigma^2) $$
def l(x):
return np.sin(x)
#Gaussian Distribution with mean -1 and variance 0.8
gaussianParameter = eq.Parameter(lower=-np.inf, upper=np.inf, order=order, distribution='Gaussian', shape_parameter_A = -1, shape_parameter_B=0.8)
uniBasis = eq.Basis('univariate')
gaussianPoly = eq.Poly(gaussianParameter, uniBasis, method='numerical-integration')
#Returning points and weights as previously done.
gaussianPoints, gaussianWeights = gaussianPoly.get_points_and_weights()
evalGaussian = float(np.dot(gaussianWeights , eq.evaluate_model(gaussianPoints, l) ) )
print(evalGaussian)
-0.5640548692740877
And thus... $$ L(x) = \int_{-\infty}^{\infty}sin(x)\rho(x) \approx -0.5640548692740783 $$
Given that in the real world, we are typically interested in functions of multiple variables.
Equadratures allows this through the specification of the multi-index, the possible multi-indices are called index sets.
The following are the 6 index sets in Equadratures: tensor-grid, total-order, sparse-grid, hyperbolic-basis, and euclidean degree.
For the pruposes of this tutorial, we will begin with a tensor-grid.
tensor = eq.Basis('tensor-grid', [20,20])
biParam = eq.Parameter(lower=-1, upper=1, order=20, distribution='Uniform')
biPoly = eq.Poly([biParam, biParam], tensor, method='numerical-integration')
biPoints, biWeights = biPoly.get_points_and_weights()
Now consider the multivariable function: $$ G(x) = \int_{-1}^{1}\int_{-1}^{1}16x^3 + 8y^2 - 8y+6 \;dxdy $$
def g(x):
return 16*x[0]**3+8*x[1]**2-8*x[1]+6;
xValuesBi = np.linspace(-1., 1., 1000)
yValuesBi = np.linspace(-1., 1., 1000)
x, y, z = np.transpose(biPoints)[0], np.transpose(biPoints)[1], biWeights
fig = go.Figure(data=[go.Scatter3d(
x=x,
y=y,
z=z,
mode='markers',
marker=dict(
size=8,
color=z,
colorscale='Viridis',
opacity=0.8
)
)])
fig.show()
The numerical integral is now easily calculated using the points and weights we had earlier computed.
evalBi = float( 2 * np.dot(biWeights , eq.evaluate_model(biPoints, g) ) )
print(evalBi);
17.33333333333335
$$ G(x) = \int_{-1}^{1}\int_{-1}^{1}16x^3 + 8y^2 - 8y+6 \;dxdy = 17.\overline{33} $$
Now we will similarly approach the problem using a sparse grid.
sparseGrid = eq.Basis('sparse-grid', level=5, growth_rule='linear')
sparsePoly = eq.Poly([biParam, biParam], sparseGrid, method='numerical-integration')
sparsePoints, sparseWeights = sparsePoly.get_points_and_weights();
x1, y1, z1 = np.transpose(sparsePoints)[0], np.transpose(sparsePoints)[1], sparseWeights
fig = go.Figure(data=[go.Scatter3d(
x=x1,
y=y1,
z=z1,
mode='markers',
marker=dict(
size=5,
color=z,
colorscale='Viridis',
opacity=0.8
)
)])
fig.update_layout(title='Weights vs Points')
fig.show()
Once more, the numeric value is then easily computed
evalSparse = float( 2 * np.dot(sparseWeights , eq.evaluate_model(sparsePoints, g) ) )
print(evalBi);
17.33333333333335
$$ G(x) = \int_{-1}^{1}\int_{-1}^{1}16x^3 + 8y^2 - 8y+6 \;dxdydz = 17.\overline{33} $$
When considering the same bivariate function, we can utilize the true strengths of equadratures to compute a more computationally efficient integral. Using our effectively subsampled quadratures grid we are able to select a sampling grid that is able to both accurately and efficiently compute the very same integral, as we will demonstrate below.
$$ G(x) = \int_{-1}^{1}\int_{-1}^{1}16x^3 + 8y^2 - 8y+6 \;dxdy $$
X1 = eq.Parameter(lower=-1, upper=1, order = 3)
basis = eq.Basis('total-order', orders=[5,10])
#basis.plot_index_set()
effSubPoly = eq.Poly([X1,X1], basis, method='least-squares',
sampling_args = {'mesh': 'tensor-grid', 'subsampling-algorithm': 'qr'})
effSubPoints, effSubWeights = effSubPoly.get_points_and_weights()
evalEffSub = float( 2 * np.dot(effSubWeights , eq.evaluate_model(effSubPoints, g) ) )
print(evalEffSub)
17.33333333333333
Resulting in:
$$ G(x) = \int_{-1}^{1}\int_{-1}^{1}16x^3 + 8y^2 - 8y+6 \;dxdy = 17.\overline{33} $$
As seen above, allowing equadratures to select a subsampling grid will enable it to more efficiently compute the desired integral.