Using OpenQuad#
This guide describes the basic steps for using the key functionality of OpenQuad, illustrated with a concrete example. If you are not sure how to use this package, this guide is a good starting point.
Initialization#
As a first step, import the class for the desired geometry. For example,
openquad.S2
for the integral over the surface of the 2d unit sphere,
\(\mathrm{S}^2\):
>>> import numpy as np
>>> from openquad import S2
We will also use NumPy later on.
Initiate the class with the method specification (see S2
) in
the form of a list of tuples, containing the a string with the name of the method
and a dictionary with required parameters.
Either use a method tailored to this geometry
>>> quad = S2([('S2-Design-Womersley', degree=7)])
or combine lower-dimensional methods with compatible geometries
>>> quad = S2([
... ('GaussLegendre', dict(degree=7))]),
... ('Trapezoid', dict(size=8)),
... ], polar_sampling='cos')
Some classes take additional parameters, like polar_sampling
for
openquad.S2
.
Common functionality#
All geometry classes share a set of common attributes and functions, supplemented by functionality specific to certain geometries.
Sample points and weights#
Sample points and weights are stored as NumPy arrays in points
and
weights
, respecitvely.
>>> quad.weights.shape
(32,)
>>> quad.points.shape
(2, 32)
Weights are normalized such that their sum equals the volume of the integration domain. For \(\mathrm{S}^2\), this is the area \(4\pi\):
>>> np.sum(quad.weights) / (4 * np.pi)
np.float(1.0)
Sample points in points
are given in the default
coordinates of the selected integration domain. For S2
,
these are spherical polar angles. Other coordinates might be available, e.g.
angles
or xyz
.
>>> np.array_equal(quad.angles, quad.points)
True
>>> quad.xyz.shape
(3, 32)
Exporting quadratures#
You can save quadrature points and weights as a textfile with
savetxt()
.
>>> quad.savetxt('points_and_weights.txt')
Integration#
Each class is equipped with the integrate()
function, which can handle
arrays and Python callables.
Suppose the integrand \(f(x)\) is a Python function, e.g.
>>> def f(theta, phi):
... """Spherical harmonic |Y_{2, 1}|^2."""
... return ((np.sin(2*theta) * np.cos(phi))*np.sqrt(15/(16*np.pi)))**2
To perform the integral of this function directly
>>> quad.integrate(f)
np.float64(1.0000000000000002)
In some situations it may be desirable or necessary to access the function values available on the quadrature grid.
>>> f_values = f(*quad.angles)
>>> f_values.shape
(32,)
You can perform the integration on the array data at a later point with
>>> quad.integrate(f_values)
np.float64(1.0000000000000002)
Other parameters#
Other attributes that are available for all top-level classes include:
dim
: the dimension of the domain \(\mathcal{D}\).size
: the number of sample points.shape
: the shape ofpoints
.
See the API reference for details.