shenfun.laguerre package

Submodules

shenfun.laguerre.bases module

class BCGeneric(N, bc=None, domain=None, alpha=0, beta=0, **kw)

Bases: CompositeBase

Function space for setting inhomogeneous boundary conditions

Parameters:

• N (int) – Number of quadrature points in the homogeneous space.

• bc (dict) – The boundary conditions in dictionary form, see BoundaryConditions.

• domain (2-tuple of numbers, optional) – The domain of the homogeneous space.

• alpha (number, optional) – Parameter of the Jacobi polynomial

• beta (number, optional) – Parameter of the Jacobi polynomial

static boundary_condition()

static short_name()

basis_function(i=0, x=x)

Return basis function i

Parameters:

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

• x (sympy Symbol, optional)

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_derivative(x=None, i=0, k=0, output_array=None)

Evaluate k’th derivative of basis function i at x or all quadrature points

Parameters:

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

• i (int, optional) – Basis function number

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

• output_array (array, optional) – return array

Returns:

output_array

Return type:

array

get_orthogonal(**kwargs)

Return orthogonal space (otherwise as self)

Returns:

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

Return type:

SpectralBase

shape(forward_output=True)

Return the allocated shape of arrays used for self

Parameters:

forward_output (bool, optional) – If True then return allocated shape of spectral space (the result of a forward transform). If False then return allocated shape of physical space (the input to a forward transform).

slice()

Return index set of current space

stencil_matrix(N=None)

Return stencil matrix in SparseMatrix format

Parameters:

N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None

to_ortho(input_array, output_array=None)

Project to orthogonal basis

Parameters:

• input_array (array) – Expansion coefficients of input basis

• output_array (array, optional) – Expansion coefficients in orthogonal basis

Returns:

output_array

Return type:

array

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().

property dim_ortho

property is_boundary_basis

class CompactDirichlet(N, bc=(0, ), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Laguerre function space for Dirichlet boundary conditions

The basis \(\{\phi_k\}_{k=0}^{N-1}\) is

\[\begin{split}\phi_k &= (La_k - La_{k+1})\exp(-x/2), \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= L_0\exp(-x/2),\end{split}\]

such that

\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(0) &= a\end{split}\]

The last basis function is for boundary condition and only used if a is different from 0. In one dimension \(\hat{u}_{N-1}=a\).

Parameters:

• N (int) – Number of quadrature points

• quad (str, optional) – Type of quadrature

  • LG - Laguerre-Gauss

• bc (1-tuple of number (a,)) – Boundary value at x=0

• padding_factor (float, optional) – Factor for padding backward transforms.

• dealias_direct (bool, optional) – Set upper 1/3 of coefficients to zero before backward transform

• 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

static boundary_condition()

static short_name()

class CompactNeumann(N, bc=(0, ), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Laguerre function space for Dirichlet boundary conditions

The basis \(\{\phi_k\}_{k=0}^{N-1}\) is

\[\begin{split}\phi_k &= (La_k - \frac{2k+1}{2k+3}La_{k+1})\exp(-x/2), \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= La_0\exp(-x/2),\end{split}\]

such that

\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u'(0) &= a\end{split}\]

Parameters:

• N (int) – Number of quadrature points

• quad (str, optional) – Type of quadrature

  • LG - Laguerre-Gauss

• bc (1-tuple of number (a,)) – Boundary value a = u’(0)

• padding_factor (float, optional) – Factor for padding backward transforms.

• dealias_direct (bool, optional) – Set upper 1/3 of coefficients to zero before backward transform

• 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

static boundary_condition()

static short_name()

class Generic(N, quad='LG', bc={}, domain=(0, oo), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space for Laguerre space with any boundary conditions

Any combination of Dirichlet and Neumann is possible.

Parameters:

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

• quad (str, optional) – Type of quadrature

  • LG - Laguerre-Gauss

• bc (dict, optional) – The dictionary must have key ‘left’ (not ‘right’), to describe boundary conditions on the left boundary. Specify Dirichlet with

  {‘left’: {‘D’: a}}

  for some value a, that will be neglected in the current function. Specify mixed Neumann and Dirichlet as

  {‘left’: {‘D’: a, ‘N’: b}}

  See BoundaryConditions.

• domain (2-tuple of floats, optional) – The computational domain

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

• padding_factor (float, optional) – Factor for padding backward transforms.

• dealias_direct (bool, optional) – Set upper 1/3 of coefficients to zero before backward transform

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

static boundary_condition()

static short_name()

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

Bases: SpectralBase

Function space for a regular Laguerre series

The orthogonal basis is the Laguerre function

\[\phi_k = La_k \exp(-x/2), \quad k = 0, 1, \ldots, N-1,\]

where \(La_k\) is the \(k\)’th Laguerre polynomial.

Parameters:

• N (int) – Number of quadrature points

• padding_factor (float, optional) – Factor for padding backward transforms.

• dealias_direct (bool, optional) – Set upper 1/3 of coefficients to zero before backward transform

• 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 Laguerre functions and not the regular Laguerre polynomials as basis functions.

static bnd_values(k=0, **kw)

static family()

static short_name()

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

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_derivative(x=None, i=0, k=0, output_array=None)

Evaluate k’th derivative of basis function i at x or all quadrature points

Parameters:

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

• i (int, optional) – Basis function number

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

• output_array (array, optional) – return array

Returns:

output_array

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

get_bc_space()

get_orthogonal(**kwargs)

Return orthogonal space (otherwise as self)

Returns:

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

Return type:

SpectralBase

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()

stencil_matrix(N=None)

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().

property is_orthogonal

shenfun.laguerre.matrices module

class ACDCDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

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

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

\(\phi_k \in\) laguerre.bases.CompactDirichlet 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.

class BCDCDmat(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}=(\phi_j, \phi_k)_w,\]

\(\phi_k \in\) laguerre.bases.CompactDirichlet 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.

class BLLmat(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}=(La_j, La_k)_w,\]

\(La_k \in\) laguerre.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