shenfun.jacobi package¶
Submodules¶
shenfun.jacobi.bases module¶
Module for function spaces of generalized Jacobi type
- class shenfun.jacobi.bases.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
- 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)
- static boundary_condition()¶
- property dim_ortho¶
- eval(x, u, output_array=None)¶
Evaluate
Function
u
at positionx
- 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
atx
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:
- property is_boundary_basis¶
- 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).
- static short_name()¶
- 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 ofx
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()
.
- class shenfun.jacobi.bases.CompactDirichlet(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Dirichlet boundary conditions
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= P^{(\alpha,\beta)}_k + a_k P^{(\alpha,\beta)}_{k+1} + b_k P^{(\alpha,\beta)}_{k+2} \quad k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]where
\[\begin{split}a_k &= \frac{\left(\alpha - \beta\right) \left(n + 1\right) \left(\alpha + \beta + 2 n + 3\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + 2 n + 4\right)}, \\ b_k &= - \frac{\left(n + 1\right) \left(n + 2\right) \left(\alpha + \beta + 2 n + 2\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + 2 n + 4\right)}, \\\end{split}\]and the expansion is
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a \text{ and } u(1) = b.\end{split}\]The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.CompactNeumann(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Neumann boundary conditions
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= P^{(\alpha,\beta)}_k + a_k P^{(\alpha,\beta)}_{k+1} + b_k P^{(\alpha,\beta)}_{k+2} \quad k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]where
\[\begin{split}a_k &= \frac{n \left(\alpha - \beta\right) \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + 2 n + 3\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + n + 2\right) \left(\alpha + \beta + 2 n + 4\right)}, \\ b_k &= - \frac{n \left(n + 1\right) \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + 2 n + 2\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + n + 3\right) \left(\alpha + \beta + 2 n + 4\right)}, \\\end{split}\]and the expansion is
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u'(-1) &= a \text{ and } u'(1) = b.\end{split}\]The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.Generic(N, quad='JG', bc={}, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, alpha=0, beta=0, **kw)[source]¶
Bases:
CompositeBase
Function space for 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 - JG - Jacobi-Gauss
bc (dict, optional) – The dictionary must have keys ‘left’ and ‘right’, to describe boundary conditions on the left and right boundaries, and a list of 2-tuples to specify the condition. Specify Dirichlet on both ends with
{‘left’: {‘D’: a}, ‘right’: {‘D’: b}}
for some values a and b, that will be neglected in the current function. Specify mixed Neumann and Dirichlet as
{‘left’: {‘N’: a}, ‘right’: {‘N’: b}}
For both conditions on the right do
{‘right’: {‘N’: a, ‘D’: b}}
Any combination should be possible, and it should also be possible to use second derivatives N2. See
BoundaryConditions
.domain (2-tuple of numbers, 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
Note
A test function is always using homogeneous boundary conditions.
- class shenfun.jacobi.bases.JacobiBase(N, quad='', alpha=0, beta=0, padding_factor=1, domain=(-1.0, 1.0), dtype=None, dealias_direct=False, coordinates=None)[source]¶
Bases:
SpectralBase
Abstract base class for Jacobi function spaces
- L2_norm_sq(i)[source]¶
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
- property is_jacobi¶
- l2_norm_sq(i=None)[source]¶
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.
- class shenfun.jacobi.bases.LowerDirichlet(N, quad='JG', bc=(0, None), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Dirichlet boundary condition at left edge
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= P^{(\alpha, \beta)}_k + \frac{(k+1)}{\beta+k+1} P^{(\alpha,\beta)}_{k+1} \quad k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= 1\end{split}\]and the expansion is
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) = 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
JG - Jacobi-Gauss
bc (2-tuple of (number, None), optional) – Boundary condition at x=-1.
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.Orthogonal(N, quad='JG', alpha=0, beta=0, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
JacobiBase
Function space for regular (orthogonal) Jacobi functions
The orthogonal basis is
\[\phi_k = P^{(\alpha,\beta)}_k, \quad k = 0, 1, \ldots, N-1,\]where \(P^{(\alpha,\beta)}_k\) is the Jacobi polynomial. The basis is orthogonal with weight \((1-x)^{\alpha}(1+x)^{\beta}\).
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- eval(x, u, output_array=None)[source]¶
Evaluate
Function
u
at positionx
- 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)[source]¶
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)[source]¶
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(x=None, i=0, k=0, output_array=None)[source]¶
Evaluate k’th derivative of basis function
i
atx
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)[source]¶
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_orthogonal(**kwargs)[source]¶
Return orthogonal space (otherwise as self)
- Returns:
The orthogonal space in the same family, and otherwise as self.
- Return type:
- property is_orthogonal¶
- orthogonal_basis_function(i=0, x=x)[source]¶
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)[source]¶
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.
- vandermonde(x)[source]¶
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 ofx
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()
.
- class shenfun.jacobi.bases.Phi1(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Dirichlet boundary conditions
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= \frac{(1-x^2)}{h^{(1,\alpha,\beta)}_{k+1}} \frac{dP^{(\alpha,\beta)}_{k+1}}{dx}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]where
\[\begin{split}h^{(1,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{dP^{(\alpha,\beta)}_{k}}{dx}\right)^2 (1-x)^{\alpha+1}(1+x)^{\beta+1}dx, \\ &= \frac{2^{\alpha+\beta+1}(\alpha+\beta+k+1) \Gamma{(\alpha+k+1)} \Gamma{(\beta+k+1)}}{(\alpha+\beta+2k+1) \Gamma{(k)} \Gamma{(\alpha+\beta+k+1)}},\end{split}\]and the expansion is
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a \text{ and } u(1) = b.\end{split}\]The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.Phi2(N, quad='JG', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for biharmonic equations
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-5\) are
\[\phi_k &= \frac{(1-x^2)^2}{h^{(2,\alpha,\beta)}_{k+2}} \frac{d^2P^{(\alpha,\beta)}_{k+2}}{dx^2},\]where
\[\begin{split}h^{(2,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^2P^{(\alpha,\beta)}_{k}}{dx^2}\right)^2 (1-x)^{\alpha+2}(1+x)^{\beta+2}dx, \\ &= \frac{2^{\alpha+\beta+1}(\alpha+\beta+k+1) \Gamma{(\alpha+k+1)} \Gamma{(\beta+k+1)}}{(\alpha+\beta+2k+1) \Gamma{(k)} \Gamma{(\alpha+\beta+k+1)}},\end{split}\]The 4 boundary basis functions are computed using
jacobi.findbasis.get_bc_basis()
, but they are too messy to print here. We have\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, u'(-1)=b, u(1)=c \text{ and } u'(1) = d.\end{split}\]The last 4 basis functions are for boundary conditions and only used if a, b, c or d are different from 0. In one dimension \(\hat{u}_{N-4}=a\), \(\hat{u}_{N-3}=b\), \(\hat{u}_{N-2}=c\) and \(\hat{u}_{N-1}=d\).
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (4-tuple of numbers, optional) – Boundary conditions.
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.Phi3(N, quad='JG', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 6th order equations
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-7\) are
\[\begin{split}\phi_k &= \frac{(1-x^2)^3}{h^{(3,\alpha,\beta)}_{k+3}} \frac{d^3P^{(\alpha,\beta)}_{k+3}}{dx^3}, \\\end{split}\]where
\[\begin{split}h^{(3,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^3P^{(\alpha,\beta)}_{k}}{dx^3}\right)^2 (1-x)^{\alpha+3}(1+x)^{\beta+3}dx, \\\end{split}\]The 6 boundary basis functions are computed using
jacobi.findbasis.get_bc_basis()
, but they are too messy to print here. We have\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, u'(-1)=b, u''(-1)=c, u(1)=ed u'(1)=e, u''(1)=f.\end{split}\]The last 6 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (6-tuple of numbers, optional) – Boundary conditions.
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.Phi4(N, quad='JG', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 8th order equations
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-9\) are
\[\begin{split}\phi_k &= \frac{(1-x^2)^4}{h^{(4,\alpha,\beta)}_{k+4}} \frac{d^4P^{(\alpha,\beta)}_{k+4}}{dx^4}, \\\end{split}\]where
\[\begin{split}h^{(4,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^4P^{(\alpha,\beta)}_{k}}{dx^4}\right)^2 (1-x)^{\alpha+4}(1+x)^{\beta+4}dx, \\ &=\frac{2^{\alpha + \beta + 1} \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + n + 2\right) \left(\alpha + \beta + n + 3\right) \left(\alpha + \beta + n + 4\right) \Gamma\left(\alpha + n + 1\right) \Gamma\left(\beta + n + 1\right)}{\left(\alpha + \beta + 2 n + 1\right) \Gamma\left(n - 3\right) \Gamma\left(\alpha + \beta + n + 1\right)},\end{split}\]The 8 boundary basis functions are computed using
jacobi.findbasis.get_bc_basis()
, but they are too messy to print here. We have\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, u'(-1)=b, u''(-1)=c, u''''(-1)=d, u(1)=e, u'(1)=f, u''(1)=g, u''''(1)=h.\end{split}\]The last 8 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
JG - Jacobi-Gauss
bc (8-tuple of numbers, optional) – Boundary conditions.
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
- class shenfun.jacobi.bases.UpperDirichlet(N, quad='JG', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Dirichlet boundary conditions at right edge
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= P^{(\alpha, \beta)}_k - \frac{(k+1)}{\alpha+k+1} P^{(\alpha,\beta)}_{k+1} \quad k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= 1\end{split}\]and the expansion is
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(1) = 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
JG - Jacobi-Gauss
bc (2-tuple of (None, number), optional) – Boundary condition at x=1.
alpha (number, optional) – Parameter of the Jacobi polynomial
beta (number, optional) – Parameter of the Jacobi polynomial
domain (2-tuple of numbers, optional) – The computational domain
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
shenfun.jacobi.matrices module¶
- class shenfun.jacobi.matrices.BJJmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Mass matrix \(B=(b_{kj}) \in \mathbb{R}^{M \times N}\), where
\[b_{kj}=(J_j, J_k)_w,\]\(J_k \in\)
jacobi.bases.Orthogonal
and test and trial spaces have dimensions of M and N, respectively.- assemble(method)[source]¶
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 oflegendre.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.
shenfun.jacobi.recursions module¶
Recursions for Jacobi polynomials, or standardized Jacobi polynomials
- shenfun.jacobi.recursions.Lmat(k, q, l, M, N, alf=0, bet=0, gn=1)[source]¶
Return matrix corresponding to
\[(\partial^{k-l}Q_{n}, x^q \phi^{(k)}_m)_{\omega}\quad (1)\]where \(\partial^k\) represents the \(k\)’th derivative and
\[\begin{split}Q_n = g_n^{(\alpha, \beta)} P^{(\alpha, \beta)}_n \\ \phi^{(k)}_m = \frac{(1-x^2)^k \partial^k Q_{m+k}}{h^{(k)}_{m+k}}\end{split}\]- Parameters:
k, q, l (integers) – Numbers in variational form (1)
M, N (int) – Shape of matrix
alf, bet (numbers, optional) – Jacobian parameters
gn (scaling function) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
Example
>>> from shenfun.jacobi.recursions import Lmat, cn, half >>> X = Lmat(2, 2, 0, 10, 12, -half, -half, cn)
- shenfun.jacobi.recursions.ShiftedMatrix(mat, q, r, s, M=0, N=0, k=None, alf=0, bet=0, gn=1)[source]¶
Index-shifted q’th power of matrix mat. Either
\[\begin{split}A^{(k, q)}_{(r, s)} = (a^{(k, q)}_{m+r,n+s})_{m,n=0}^{M,N} \\ B^{(q)}_{(r, s)} = (b^{(q)}_{m+r,n+s})_{m,n=0}^{M,N}\end{split}\]- Parameters:
mat (Python function or SparseMatrix) – The Python function must have signature (k, q, alf, bet, i, j, gn=1) or (alf, bet, i, j, gn=1) if k is None
q (int) – The matrix power
r (int) – Shift in row
s (int) – Shift in column
M, N (int, optional) – The shape of the final matrix
k (int or None, optional) – Parameter of recursion of the k’th derivative This is only used if the recursion matrix is A.
alf, bet (numbers, optional) – Jacobi paramters
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
Note
If mat is a SparseMatrix, then the resulting shifted matrix will be smaller than mat.
- shenfun.jacobi.recursions.a(alf, bet, i, j, gn=1)[source]¶
Recursion matrix \(A\) for standardized Jacobi polynomials
The recursion is
\[x \boldsymbol{Q} = {A}^T \boldsymbol{Q}\]where
\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
i, j (int) – Indices for row and column
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.a_(k, q, alf, bet, i, j, gn=1)[source]¶
Recursion matrix \(\underline{A}\) for standardized Jacobi polynomials
The recursion is
\[x \partial^k \boldsymbol{Q} = \underline{A}^T \partial^k \boldsymbol{Q} \quad (*)\]where \(\partial^k\) represents the \(k\)’th derivative and
\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]- Parameters:
k, q (int) – Parameters in (*)
alf, bet (numbers) – Jacobi parameters
i, j (int) – Indices for row and column
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.a_mat(mat, k, q, alf, bet, M, N, gn=1)[source]¶
Return SparseMatrix of recursion matrix
\[A^{(k,q)} = (a^{(k,q)}_{mn})_{m,n=0}^{M, N}\]- Parameters:
mat (Python function for matrix) – The Python function must have signature (k, q, alf, bet, i, j, gn=1)
k (int) – Parameter
q (int) – Matrix power
alf, bet (numbers) – Jacobi parameters
M, N (int) – Shape of returned matrix
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.b(alf, bet, i, j, gn=1)[source]¶
Recursion matrix \(B\) for standardized Jacobi polynomials
The recursion is
\[\boldsymbol{Q} = {B}^T \partial \boldsymbol{Q}\]where \(\partial\) represents the derivative and
\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
i, j (int) – Indices for row and column
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.bnd_values(alf, bet, k=0, gn=1)[source]¶
Return lambda function for computing boundary values
\[\frac{d^k}{dx^k}Q_n(\pm 1)\]where \(Q^{(\alpha,\beta)}_n = g^{(\alpha,\beta)}_n P_n^{(\alpha,\beta)}\).
- Parameters:
alf, bet (Numbers) – Jacobi parameters
k (int, optional) – Number of derivatives
gn (scaling function, optional)
- shenfun.jacobi.recursions.c(alf, bet, i, j, gn=1)[source]¶
Recursion matrix \(C\) for standardized Jacobi polynomials
The recursion is
\[(1-x^2) \partial \boldsymbol{Q} = {C}^T \boldsymbol{Q}\]where \(\partial\) represents the derivative and
\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
i, j (int) – Indices for row and column
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.cn(alf, bet, n)[source]¶
Scaling function
\[c_n^{(\alpha, \beta)} = \frac{1}{P^{(\alpha,\beta)}_n(1)}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
Note
Used by Legendre and Chebyshev polynomials of first kind.
- shenfun.jacobi.recursions.djt(alf, bet, n, k)[source]¶
k’th derivative of Jacobi polynomial
- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
k (int) – The k’th derivative
- shenfun.jacobi.recursions.dqn(alf, bet, n, k, gn=<function cn>)[source]¶
Derivative of specialized Jacobi polynomial
\[\frac{d^k Q_n^{(\alpha, \beta)}}{dx^k} = g_n^{(\alpha, \beta)} \frac{d^k P^{(\alpha, \beta)}_n}{dx^k}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
k (int) – The k’th derivative
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.gamma(alf, bet, n)[source]¶
Return normalization factor \(h_n\) for inner product of Jacobi polynomials
\[h_n = (P^{(\alpha,\beta)}_n, P^{(\alpha,\beta)}_n)_{\omega^{(\alpha,\beta)}}\]- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
- shenfun.jacobi.recursions.h(alf, bet, n, k, gn=1)[source]¶
Return normalization factor \(h^{(k)}_n\) for inner product of derivatives of Jacobi polynomials
\[\begin{split}Q_n(x) = g_n(x)P^{(\alpha,\beta)}_n(x) \\ h_n^{(k)} = (\partial^k Q_n, \partial^k Q_n)_{\omega^{(\alpha+k,\beta+k)}} \quad (*)\end{split}\]where \(\partial^k\) represents the \(k\)’th derivative.
- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
k (int) – For derivative of k’th order, see (*)
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.jt(alf, bet, n)[source]¶
Jacobi polynomial
- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
- shenfun.jacobi.recursions.matpow(mat, q, alf, bet, i, j, gn=1)[source]¶
Compute and return component of q’th matrix power of mat
- Parameters:
mat (Python function (a, b or c))
q (int) – matrix power
alf, bet (numbers) – Jacobi parameters
i, j (int) – Row and column indices
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.pmat(mat, q, alf, bet, M, N, gn=1)[source]¶
Return SparseMatrix of q’th matrix power of recursion matrix mat
- Parameters:
mat (Python function for matrix) – The Python function must have signature (alf, bet, i, j, gn=1)
q (int) – Matrix power
alf, bet (numbers) – Jacobi parameters
M, N (int) – Shape of returned matrix
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.
- shenfun.jacobi.recursions.psi(alf, bet, n, k)[source]¶
Normalization factor for
\[\partial^k P^{(\alpha, \beta)}_n = \psi^{(k,\alpha,\beta)}_{n} P^{(\alpha+k,\beta+k)}_{n-k}, \quad n \ge k, \quad (*)\]where \(\partial^k\) represents the \(k\)’th derivative
- Parameters:
alf, bet (numbers) – Jacobi parameters
n, k (int) – Parameters in (*)
- shenfun.jacobi.recursions.qn(alf, bet, n, gn=<function cn>)[source]¶
Specialized Jacobi polynomial
\[Q_n^{(\alpha, \beta)} = g_n^{(\alpha, \beta)} P^{(\alpha, \beta)}_n\]- Parameters:
alf, bet (numbers) – Jacobi parameters
n (int) – Index
gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.