shenfun.hermite package

Submodules

shenfun.hermite.bases module

class shenfun.hermite.bases.Orthogonal(N, dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: SpectralBase

Function space for Hermite functions

The orthogonal basis is the Hermite function

\[\phi_k = H_k \frac{1}{\pi^{0.25} \sqrt{2^n n!}} \exp(-x^2/2), \quad k = 0, 1, \ldots, N-1,\]

where \(\phi_k\) and \(H_k\) are the Hermite function and Hermite polynomials of order k, respectively.

Parameters:

• N (int, optional) – Number of quadrature points. Should be even for efficiency, but this is not required.

• padding_factor (float, optional) – Factor for padding backward transforms. padding_factor=1.5 corresponds to a 3/2-rule for dealiasing.

• dealias_direct (bool, optional) – True for dealiasing using 2/3-rule. Must be used with padding_factor = 1.

• dtype (data-type, optional) – Type of input data in real physical space. Will be overloaded when basis is part of a TensorProductSpace.

• coordinates (2- or 3-tuple (coordinate, position vector (, sympy assumptions)), optional) – Map for curvilinear coordinatesystem, and parameters to Coordinates

Note

We are using Hermite functions and not the regular Hermite polynomials as basis functions.

L2_norm_sq(i)

Return square of L2-norm

\[\| \phi_i \|^2_{\omega} = (\phi_i, \phi_i)_{\omega} = \int_{I} \phi_i \overline{\phi}_i \omega dx\]

where \(\phi_i\) is the i’th orthogonal basis function for the orthogonal basis in the given family.

Parameters:

i (int) – The number of the orthogonal basis function

static boundary_condition()

domain_factor()

Return scaling factor for domain

Note

The domain factor is the length of the reference domain over the length of the true domain.

eval(x, u, output_array=None)

Evaluate Function u at position x

Parameters:

• x (float or array of floats)

• u (array) – Expansion coefficients or instance of Function

• output_array (array, optional) – Function values at points

Returns:

output_array

Return type:

array

evaluate_basis(x, i=0, output_array=None)

Evaluate basis function i at points x

Parameters:

• x (float or array of floats)

• i (int, optional) – Basis function number

• output_array (array, optional) – Return result in output_array if provided

Returns:

output_array

Return type:

array

evaluate_basis_all(x=None, argument=0)

Evaluate basis at x or all quadrature points

Parameters:

• x (float or array of floats, optional) – If not provided use quadrature points of self

• argument (int) – Zero for test and 1 for trialfunction

Returns:

Vandermonde matrix

Return type:

array

evaluate_basis_derivative_all(x=None, k=0, argument=0)

Return k’th derivative of basis evaluated at x or all quadrature points as a Vandermonde matrix.

Parameters:

• x (float or array of floats, optional) – If not provided use quadrature points of self

• k (int, optional) – k’th derivative

• argument (int) – Zero for test and 1 for trialfunction

Returns:

Vandermonde matrix

Return type:

array

static factor(i)

static family()

get_orthogonal(**kwargs)

Return orthogonal space (otherwise as self)

Returns:

The orthogonal space in the same family, and otherwise as self.

Return type:

SpectralBase

property is_orthogonal

l2_norm_sq(i=None)

Return square of l2-norm

\[\| u \|^2_{N,\omega} = (u, u)_{N,\omega} = \sun_{j=0}^{N-1} u(x_j) \overline{u}(x_j) \omega_j\]

where \(u=\{\phi_i\}_{i=0}^{N-1}\) and \(\phi_i\) is the i’th orthogonal basis function in the given family.

Parameters:

i (None or int) – If None then return the square of the l2-norm for all i=0, 1, …, N-1. Else, return for given i.

orthogonal_basis_function(i=0, x=x)

Return the orthogonal basis function i

Parameters:

• i (int, optional) – The degree of freedom of the basis function

• x (sympy Symbol, optional)

points_and_weights(N=None, map_true_domain=False, weighted=True, **kw)

Return points and weights of quadrature for weighted integral

\[\int_{\Omega} f(x) w(x) dx \approx \sum_{i} f(x_i) w_i\]

Parameters:

• N (int, optional) – Number of quadrature points

• map_true_domain (bool, optional) – Whether or not to map points to true domain

• weighted (bool, optional) – Whether to use quadrature weights for a weighted inner product (default), or a regular, non-weighted inner product.

Note

The weight of the space is included in the returned quadrature weights.

reference_domain()

static short_name()

vandermonde(x)

Return Vandermonde matrix based on the primary (orthogonal) basis of the family.

Evaluates basis function \(\psi_k(x)\) for all wavenumbers, and all x. Returned Vandermonde matrix is an N x M matrix with N the length of x and M the number of bases.

\[\begin{split}\begin{bmatrix} \psi_0(x_0) & \psi_1(x_0) & \ldots & \psi_{M-1}(x_0)\\ \psi_0(x_1) & \psi_1(x_1) & \ldots & \psi_{M-1}(x_1)\\ \vdots & \ldots \\ \psi_{0}(x_{N-1}) & \psi_1(x_{N-1}) & \ldots & \psi_{M-1}(x_{N-1}) \end{bmatrix}\end{split}\]

Parameters:

x (array of floats) – points for evaluation

Note

This function returns a matrix of evaluated orthogonal basis functions for either family. The true Vandermonde matrix of a basis is obtained through SpectralBase.evaluate_basis_all().

shenfun.hermite.matrices module

class shenfun.hermite.matrices.AHHmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where

\[a_{kj}=(H'_j, H'_k)_w,\]

\(H_k \in\) hermite.bases.Orthogonal and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

Subclass SpectralMatrix and overload this method in order to provide a fast and accurate implementation of the matrix representing an inner product. See the matrix modules in either one of

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

\[(T_j, T_i)_{\omega} = \frac{c_i \pi}{2}\delta_{ij},\]

where \(c_0=2\) and \(c_i=1\) for integer \(i>0\). We can implement this as

>>> from shenfun import SpectralMatrix
>>> class Bmat(SpectralMatrix):
...     def assemble(self, method):
...         test, trial = self.testfunction, self.trialfunction
...         ci = np.ones(test[0].N)
...         ci[0] = 2
...         if test[0].quad == 'GL' and method != 'exact':
...             # Gauss-Lobatto quadrature inexact at highest polynomial order
...             ci[-1] = 2
...         return {0: ci*np.pi/2}

Here {0: ci*np.pi/2} is the 0’th diagonal of the matrix. Note that test and trial are two-tuples of (instance of :class:.SpectralBase`, number)`, where the number represents the number of derivatives. For the mass matrix the number will be 0. Also note that the length of the diagonal must be correct.

get_solver()

Return appropriate solver for self

Note

Fall back on generic Solve, which is using Scipy sparse matrices with splu/spsolve. This is still pretty fast.

matvec(v, c, format='cython', axis=0)

Matrix vector product

Returns c = dot(self, v)

Parameters:

• v (array) – Numpy input array of ndim>=1

• c (array) – Numpy output array of same shape as v

• format (str, optional) – Choice for computation

  • csr - Compressed sparse row format

  • dia - Sparse matrix with DIAgonal storage

  • python - Use numpy and vectorization

  • self - To be implemented in subclass

  • cython - Cython implementation that may be implemented in subclass

  • numba - Numba implementation that may be implemented in subclass

  Using config['matrix']['sparse']['matvec'] setting if format is None

• axis (int, optional) – The axis over which to take the matrix vector product

class shenfun.hermite.matrices.BHHmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

Mass matrix \(B=(b_{kj}) \in \mathbb{R}^{M \times N}\), where

\[b_{kj}=(H_j, H_k)_w,\]

\(H_k \in\) hermite.bases.Orthogonal and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

Subclass SpectralMatrix and overload this method in order to provide a fast and accurate implementation of the matrix representing an inner product. See the matrix modules in either one of

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

\[(T_j, T_i)_{\omega} = \frac{c_i \pi}{2}\delta_{ij},\]

where \(c_0=2\) and \(c_i=1\) for integer \(i>0\). We can implement this as

>>> from shenfun import SpectralMatrix
>>> class Bmat(SpectralMatrix):
...     def assemble(self, method):
...         test, trial = self.testfunction, self.trialfunction
...         ci = np.ones(test[0].N)
...         ci[0] = 2
...         if test[0].quad == 'GL' and method != 'exact':
...             # Gauss-Lobatto quadrature inexact at highest polynomial order
...             ci[-1] = 2
...         return {0: ci*np.pi/2}

Here {0: ci*np.pi/2} is the 0’th diagonal of the matrix. Note that test and trial are two-tuples of (instance of :class:.SpectralBase`, number)`, where the number represents the number of derivatives. For the mass matrix the number will be 0. Also note that the length of the diagonal must be correct.

matvec(v, c, format=None, axis=0)

Matrix vector product

Returns c = dot(self, v)

Parameters:

• v (array) – Numpy input array of ndim>=1

• c (array) – Numpy output array of same shape as v

• format (str, optional) – Choice for computation

  • csr - Compressed sparse row format

  • dia - Sparse matrix with DIAgonal storage

  • python - Use numpy and vectorization

  • self - To be implemented in subclass

  • cython - Cython implementation that may be implemented in subclass

  • numba - Numba implementation that may be implemented in subclass

  Using config['matrix']['sparse']['matvec'] setting if format is None

• axis (int, optional) – The axis over which to take the matrix vector product

solve(b, u=None, axis=0, constraints=())

Solve matrix system Au = b

where A is the current matrix (self)

Parameters:

• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multi- dimensional

• constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

Note

Vectors may be one- or multidimensional.

Module contents