shenfun.chebyshev package¶
Submodules¶
shenfun.chebyshev.bases module¶
Module for defining function spaces in the Chebyshev family.
A function is approximated in the Chebyshev basis as
where \(T_i(x)\) is the i’th Chebyshev polynomial of the first kind. The Chebyshev polynomials are orthogonal with weight \(\omega=1/\sqrt{1-x^2}\)
where \(c_0=2\) and \(c_i=1\) for \(i>0\).
All other bases defined in this module are combinations of \(T_i\)’s. For example, a Dirichlet basis is
The basis is implemented using a stencil matrix \(K \in \mathbb{R}^{N-2 \times N}\), such that
where \(\boldsymbol{\phi}=(\phi_0, \phi_1, \ldots, \phi_{N-3})\) and \(\boldsymbol{T}=(T_0, T_1, \ldots, T_{N-1})\). For the Dirichlet basis \(K = (\delta_{i, j} - \delta_{i+2, j})_{i,j=0}^{N-2, N}\).
All composite bases make use of the fast transforms that exists for \(\boldsymbol{T}\) through fast cosine transforms. The stencil matrix is used to transfer any composite basis back and forth to the orthogonal basis.
- class shenfun.chebyshev.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.chebyshev.bases.CombinedShenNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, 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 = \begin{cases} T_0, \quad &k=0, \\ T_1 - T_3/9, \quad &k=1, \\ T_2/4 - T_4/16, \quad &k=2, \\ -\frac{T_{k-2}}{(k-2)^2} +2\frac{T_k}{k^2} - \frac{T_{k+2}}{(k+2)^2}, &k=3, 4, \ldots, N-3, \\ \frac{1}{8}(4T_1-T_2), \quad &k=N-2 \\ \frac{1}{8}(4T_1+T_2), \quad &k=N-1 \end{cases}\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary condition values at, respectively, x=(-1, 1).
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
- stencil_matrix(N=None)[source]¶
Return stencil matrix in
SparseMatrix
format- Parameters:
N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None
- sympy_stencil(i=i, j=j)[source]¶
Return stencil matrix as a Sympy matrix
- Parameters:
i, j (Sympy symbols) – indices for row and column
implicit (bool or str, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the stencil. This makes the matrix prettier, and it can still be evaluated. If implicit is not False, then it must be a string of length one. This string represents the diagonals of the matrix. If implicit=’a’ and the stencil matrix has two diagonals in the main diagonal and the first upper diagonal, then these are called ‘a0’ and ‘a1’.
Example
>>> from shenfun import FunctionSpace >>> import sympy as sp >>> i, j = sp.symbols('i,j', integer=True) >>> D = FunctionSpace(8, 'L', bc=(0, 0), scaled=True) >>> D._stencil {0: 1/sqrt(4*n + 6), 2: -1/sqrt(4*n + 6)} >>> D.sympy_stencil() KroneckerDelta(i, j)/sqrt(4*i + 6) - KroneckerDelta(j, i + 2)/sqrt(4*i + 6) >>> D.sympy_stencil(implicit='a') KroneckerDelta(i, j)*a0(i) + KroneckerDelta(j, i + 2)*a2(i)
Get the main diagonal
>>> D.sympy_stencil(implicit=False).subs(j, i) 1/sqrt(4*i + 6)
- class shenfun.chebyshev.bases.Compact3(N, quad='GC', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 6’th order equation
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-7\) are
\[\begin{split}\phi_k &= \frac{h_k}{b^{(3)}_{k+3,k}}\frac{(1-x^2)^3}{h^{(3)}_{k+3}} T^{(3)}_{k+3} \\ h^{(3)}_{k+3} &= \frac{\pi (k+3) \Gamma (k+6)}{2k!} = \int_{-1}^1 T^{(3)}_k T^{(3)}_k (1-x^2)^{2.5} dx.\end{split}\]where \(T^{(3)}_k\) is the 3rd derivative of \(T_k\). This is
Phi3
scaled such that the main diagonal of the stencil matrix is unity.The boundary basis for inhomogeneous boundary conditions is too messy to print, but can be obtained using
get_bc_basis()
. 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.\end{split}\]The last 6 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.
- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (6-tuple of numbers)
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
- class shenfun.chebyshev.bases.Compact4(N, quad='GC', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 8’th order equation
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-9\) are
\[\begin{split}\phi_k &= \frac{h_k}{b^{(4)}_{k+4,k}}\frac{(1-x^2)^4}{h^{(4)}_{k+4}} T^{(4)}_{k+4} \\ h^{(4)}_{k+4} &= \int_{-1}^1 T^{(4)}_k T^{(4)}_k (1-x^2)^{3.5} dx.\end{split}\]where \(T^{(4)}_k\) is the 4rd derivative of \(T_k\). This is
Phi4
scaled such that the main diagonal of the stencil matrix is unity.The boundary basis for inhomogeneous boundary conditions is too messy to print, but can be obtained using
get_bc_basis()
. 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (8-tuple of numbers)
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
- class shenfun.chebyshev.bases.DirichletNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for mixed Dirichlet/Neumann boundary conditions
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} + \frac{4(k+1)}{2k^2+6k+5}T_{k+1} - \frac{2k^2+2k+1}{2k^2+6k+5}T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= T_0, \\ \phi_{N-1} &= T_0+T_1,\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of numbers) – Boundary condition values at x=-1 and x=1
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
- class shenfun.chebyshev.bases.Generic(N, quad='GC', bc={}, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **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
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-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.chebyshev.bases.Heinrichs(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, scaled=False, padding_factor=1, dealias_direct=False, 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 &= \alpha_k(1-x^2) T_{k}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{1}{2}(T_0-T_1), \\ \phi_{N-1} &= \frac{1}{2}(T_0+T_1),\end{split}\]such that
\[\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\). If the parameter scaled=True, then \(\alpha_k=1/(k+1)/(k+2)\), otherwise, \(\alpha_k=1\).
- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary conditions at, respectively, x=(-1, 1).
domain (2-tuple of numbers, optional) – The computational domain
dtype (data-type, optiona) – Type of input data in real physical space. Will be overloaded when basis is part of a
TensorProductSpace
.scaled (bool, optional) – Whether or not to use scaled basis
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
- stencil_matrix(N=None)[source]¶
Return stencil matrix in
SparseMatrix
format- Parameters:
N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None
- sympy_stencil(i=i, j=j)[source]¶
Return stencil matrix as a Sympy matrix
- Parameters:
i, j (Sympy symbols) – indices for row and column
implicit (bool or str, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the stencil. This makes the matrix prettier, and it can still be evaluated. If implicit is not False, then it must be a string of length one. This string represents the diagonals of the matrix. If implicit=’a’ and the stencil matrix has two diagonals in the main diagonal and the first upper diagonal, then these are called ‘a0’ and ‘a1’.
Example
>>> from shenfun import FunctionSpace >>> import sympy as sp >>> i, j = sp.symbols('i,j', integer=True) >>> D = FunctionSpace(8, 'L', bc=(0, 0), scaled=True) >>> D._stencil {0: 1/sqrt(4*n + 6), 2: -1/sqrt(4*n + 6)} >>> D.sympy_stencil() KroneckerDelta(i, j)/sqrt(4*i + 6) - KroneckerDelta(j, i + 2)/sqrt(4*i + 6) >>> D.sympy_stencil(implicit='a') KroneckerDelta(i, j)*a0(i) + KroneckerDelta(j, i + 2)*a2(i)
Get the main diagonal
>>> D.sympy_stencil(implicit=False).subs(j, i) 1/sqrt(4*i + 6)
- class shenfun.chebyshev.bases.LowerDirichlet(N, quad='GC', bc=(0, None), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space with single Dirichlet boundary condition
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} + T_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= T_0,\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a.\end{split}\]The last basis funciton 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
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (tuple of (number, None)) – Boundary conditions at edges of domain.
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.chebyshev.bases.LowerDirichletNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for both Dirichlet and Neumann boundary conditions on the left hand side
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} + \frac{4(k+1)}{2k+3}T_{k+1} + \frac{2k+1}{2k+3}T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= T_0, \\ \phi_{N-1} &= T_0+T_1,\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of numbers) – Boundary condition values at the left edge of domain
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
- class shenfun.chebyshev.bases.MikNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for Neumann boundary conditions
The basis functions \(\phi_k\) for \(k=0,1, \ldots, N-3\) are
\[\phi_k &= \frac{2}{k+1}\int (T_{k-1}-T_{k+1}),\]which (with also boundary functions) leads to the basis
\[\begin{split}\phi_k = \frac{1}{k+1} \begin{cases} T_0, &k=0, \\ 3T_1-T_3/3, &k=1, \\ T_2-T_4/4, &k=2, \\ -\frac{T_{k-2}}{k-2} + 2\frac{T_k}{k} - \frac{T_{k+2}}{k+2} , &k=3, 4, \ldots, N-3, \\ \frac{1}{8}(4T_1-T_2), \quad &k=N-2 \\ \frac{1}{8}(4T_1+T_2), \quad &k=N-1 \end{cases}\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary condition values at, respectively, x=(-1, 1).
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
- stencil_matrix(N=None)[source]¶
Return stencil matrix in
SparseMatrix
format- Parameters:
N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None
- sympy_stencil(i=i, j=j)[source]¶
Return stencil matrix as a Sympy matrix
- Parameters:
i, j (Sympy symbols) – indices for row and column
implicit (bool or str, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the stencil. This makes the matrix prettier, and it can still be evaluated. If implicit is not False, then it must be a string of length one. This string represents the diagonals of the matrix. If implicit=’a’ and the stencil matrix has two diagonals in the main diagonal and the first upper diagonal, then these are called ‘a0’ and ‘a1’.
Example
>>> from shenfun import FunctionSpace >>> import sympy as sp >>> i, j = sp.symbols('i,j', integer=True) >>> D = FunctionSpace(8, 'L', bc=(0, 0), scaled=True) >>> D._stencil {0: 1/sqrt(4*n + 6), 2: -1/sqrt(4*n + 6)} >>> D.sympy_stencil() KroneckerDelta(i, j)/sqrt(4*i + 6) - KroneckerDelta(j, i + 2)/sqrt(4*i + 6) >>> D.sympy_stencil(implicit='a') KroneckerDelta(i, j)*a0(i) + KroneckerDelta(j, i + 2)*a2(i)
Get the main diagonal
>>> D.sympy_stencil(implicit=False).subs(j, i) 1/sqrt(4*i + 6)
- class shenfun.chebyshev.bases.NeumannDirichlet(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for mixed Neumann/Dirichlet boundary conditions
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} - \frac{4(k+1)}{2k^2+6k+5}T_{k+1} - \frac{2k^2+2k+1}{2k^2+6k+5}T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= -T_0+T_1, \\ \phi_{N-1} &= T_0,\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u'(-1) &= a, 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of numbers) – Boundary condition values at x=-1 and x=1
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
- class shenfun.chebyshev.bases.Orthogonal(N, quad='GC', domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
JacobiBase
Function space for regular Chebyshev series
The orthogonal basis is
\[T_k, \quad k = 0, 1, \ldots, N-1,\]where \(T_k\) is the \(k\)’th Chebyshev polynomial of the first kind.
- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
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
- 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
- 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_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¶
- 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.
- 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)
- plan(shape, axis, dtype, options)[source]¶
Plan transform
Allocate work arrays for transforms and set up methods forward, backward and scalar_product with or without padding
- Parameters:
shape (array) – Local shape of global array
axis (int) – This base’s axis in global
TensorProductSpace
dtype (numpy.dtype) – Type of array
options (dict) – Options for planning transforms
- 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.
- to_ortho(input_array, output_array=None)[source]¶
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)[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.chebyshev.bases.Phi1(N, quad='GC', bc=(0, 0), domain=(-1.0, 1.0), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, 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{T_k - T_{k+2}}{\pi (k+1)} = \frac{2(1-x^2)}{\pi k(k+1)} T'_{k+1}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{1}{2}(T_0-T_1), \\ \phi_{N-1} &= \frac{1}{2}(T_0+T_1),\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary conditions at, respectively, x=(-1, 1).
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
- class shenfun.chebyshev.bases.Phi2(N, quad='GC', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for biharmonic equation.
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-5\) are
\[\phi_k = \frac{2(1-x^2)^2 T''_{k+2}}{\pi (k+1)(k+2)^2(k+3)} ,\]which (along with boundary functions) gives the basis
\[\begin{split}\phi_k &= \frac{1}{2 \pi (k+1)(k+2)}(T_k - \frac{2(k+2)}{k+3}T_{k+2} + \frac{k+1}{k+3}T_{k+4}), \, k=0, 1, \ldots, N-5, \\ \phi_{N-4} &= \frac{1}{16}(8T_0-9T_1+T_3), \\ \phi_{N-3} &= \frac{1}{16}(2T_0-T_1-2T_2+T_3), \\ \phi_{N-2} &= \frac{1}{16}(8T_0+9T_1-T_3), \\ \phi_{N-1} &= \frac{1}{16}(2T_0-T_1+2T_2+T_3),\end{split}\]such that
\[\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.\end{split}\]The last four bases 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (4-tuple of numbers) – The values of the 4 boundary conditions at x=(-1, 1). The two conditions at x=-1 first and then x=1. With (a, b, c, d) corresponding to bc = {‘left’: {‘D’: a, ‘N’: b}, ‘right’: {‘D’: c, ‘N’: d}}
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
- class shenfun.chebyshev.bases.Phi3(N, quad='GC', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 6’th order equation
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-7\) are
\[\begin{split}\phi_k &= \frac{(1-x^2)^3}{h^{(3)}_{k+3}} T^{(3)}_{k+3} \\ h^{(3)}_{k+3} &= \frac{\pi (k+3) \Gamma (k+6)}{2k!} = \int_{-1}^1 T^{(3)}_k T^{(3)}_k (1-x^2)^{2.5} dx.\end{split}\]where \(T^{(3)}_k\) is the 3rd derivative of \(T_k\). The boundary basis for inhomogeneous boundary conditions is too messy to print, but can be obtained using
get_bc_basis()
. 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.\end{split}\]The last 6 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.
- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (6-tuple of numbers)
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
- class shenfun.chebyshev.bases.Phi4(N, quad='GC', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space with 4 Dirichlet and 4 Neumann boundary conditions
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-9\) are
\[\begin{split}\phi_k &= \frac{(1-x^2)^4}{h^{(4)}_{k+4}} T^{(4)}_{k+4} \\ h^{(4)}_k &= \frac{\pi k \Gamma (k+4)}{2(k-4)!} = \int_{-1}^1 T^{(4)}_k T^{(4)}_k (1-x^2)^{3.5} dx,\end{split}\]where \(T^{(4)}_k\) is the 4th derivative of \(T_k\). The boundary basis for inhomogeneous boundary conditions is too messy to print, but can be obtained using
get_bc_basis()
.- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (8-tuple of numbers)
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
- class shenfun.chebyshev.bases.Phi6(N, quad='GC', bc=(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for 12th order equation
The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-13\) are
\[\begin{split}\phi_k &= \frac{(1-x^2)^6}{h^{(6)}_{k+6}} T^{(6)}_{k+6} \\ h^{(6)}_k &= \frac{\pi k (k+5)!}{2(k-6)!} = \int_{-1}^1 T^{(6)}_k T^{(6)}_k (1-x^2)^{5.5} dx,\end{split}\]where \(T^{(6)}_k\) is the 6th derivative of \(T_k\). The boundary basis for inhomogeneous boundary conditions is too messy to print, but can be obtained using
get_bc_basis()
.- Parameters:
N (int, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (12-tuple of numbers)
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
- class shenfun.chebyshev.bases.PolarDirichlet(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for polar coordinates.
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_k - T_{k+2}, \, k=0, 1 \phi_k &= T_{k-2} - T_{k+2}, \, k=2, 3, \ldots, N-3, \\ \phi_{N-2} &= \frac{1}{2}(T_0-T_1), \\ \phi_{N-1} &= \frac{1}{2}(T_0+T_1),\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1)&=a, u(1)=b\end{split}\]The last four bases 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (4-tuple of numbers) – The values of the 4 boundary conditions at x=(-1, 1). The two conditions on x=-1 first, and then x=1. With (a, b, c, d) corresponding to bc = {‘left’: [(‘D’, a), (‘N’, b)], ‘right’: [(‘D’, c), (‘N’, d)]}
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
- stencil_matrix(N=None)[source]¶
Return stencil matrix in
SparseMatrix
format- Parameters:
N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None
- sympy_stencil(i=i, j=j)[source]¶
Return stencil matrix as a Sympy matrix
- Parameters:
i, j (Sympy symbols) – indices for row and column
implicit (bool or str, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the stencil. This makes the matrix prettier, and it can still be evaluated. If implicit is not False, then it must be a string of length one. This string represents the diagonals of the matrix. If implicit=’a’ and the stencil matrix has two diagonals in the main diagonal and the first upper diagonal, then these are called ‘a0’ and ‘a1’.
Example
>>> from shenfun import FunctionSpace >>> import sympy as sp >>> i, j = sp.symbols('i,j', integer=True) >>> D = FunctionSpace(8, 'L', bc=(0, 0), scaled=True) >>> D._stencil {0: 1/sqrt(4*n + 6), 2: -1/sqrt(4*n + 6)} >>> D.sympy_stencil() KroneckerDelta(i, j)/sqrt(4*i + 6) - KroneckerDelta(j, i + 2)/sqrt(4*i + 6) >>> D.sympy_stencil(implicit='a') KroneckerDelta(i, j)*a0(i) + KroneckerDelta(j, i + 2)*a2(i)
Get the main diagonal
>>> D.sympy_stencil(implicit=False).subs(j, i) 1/sqrt(4*i + 6)
- class shenfun.chebyshev.bases.ShenBiPolar(N, quad='GC', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for the Biharmonic equation in polar coordinates
u(-1)=a, u(1)=c, u’(-1)=b and u’(1)=d
- Parameters:
N (int) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (4-tuple of numbers)
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
- stencil_matrix(N=None)[source]¶
Return stencil matrix in
SparseMatrix
format- Parameters:
N (int or None, optional) – Shape (N, N) of the stencil matrix. Using self.N if None
- sympy_stencil(i=i, j=j)[source]¶
Return stencil matrix as a Sympy matrix
- Parameters:
i, j (Sympy symbols) – indices for row and column
implicit (bool or str, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the stencil. This makes the matrix prettier, and it can still be evaluated. If implicit is not False, then it must be a string of length one. This string represents the diagonals of the matrix. If implicit=’a’ and the stencil matrix has two diagonals in the main diagonal and the first upper diagonal, then these are called ‘a0’ and ‘a1’.
Example
>>> from shenfun import FunctionSpace >>> import sympy as sp >>> i, j = sp.symbols('i,j', integer=True) >>> D = FunctionSpace(8, 'L', bc=(0, 0), scaled=True) >>> D._stencil {0: 1/sqrt(4*n + 6), 2: -1/sqrt(4*n + 6)} >>> D.sympy_stencil() KroneckerDelta(i, j)/sqrt(4*i + 6) - KroneckerDelta(j, i + 2)/sqrt(4*i + 6) >>> D.sympy_stencil(implicit='a') KroneckerDelta(i, j)*a0(i) + KroneckerDelta(j, i + 2)*a2(i)
Get the main diagonal
>>> D.sympy_stencil(implicit=False).subs(j, i) 1/sqrt(4*i + 6)
- class shenfun.chebyshev.bases.ShenBiharmonic(N, quad='GC', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for biharmonic equation.
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_n - \frac{2(k+2)}{k+3}T_{k+2}+\frac{k+1}{k+3}T_{k+4}, \, k=0, 1, \ldots, N-5, \\ \phi_{N-4} &= \frac{1}{16}(8T_0-9T_1+T_3), \\ \phi_{N-3} &= \frac{1}{16}(2T_0-T_1-2T_2+T_3), \\ \phi_{N-2} &= \frac{1}{16}(8T_0+9T_1-T_3), \\ \phi_{N-1} &= \frac{1}{16}(2T_0-T_1+2T_2+T_3),\end{split}\]such that
\[\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.\end{split}\]The last four bases 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (4-tuple of numbers) – The values of the 4 boundary conditions at x=(-1, 1). The two conditions on x=-1 first, and then x=1. With (a, b, c, d) corresponding to bc = {‘left’: [(‘D’, a), (‘N’, b)], ‘right’: [(‘D’, c), (‘N’, d)]}
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
- class shenfun.chebyshev.bases.ShenDirichlet(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, 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 &= T_k - T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{1}{2}(T_0-T_1), \\ \phi_{N-1} &= \frac{1}{2}(T_0+T_1),\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary conditions at, respectively, x=(-1, 1).
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
- class shenfun.chebyshev.bases.ShenNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, 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 &= T_{k} - \left(\frac{k}{k+2}\right)^2 T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{1}{8}(4T_1-T_2), \\ \phi_{N-1} &= \frac{1}{8}(4T_1+T_2),\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of floats, optional) – Boundary condition values at, respectively, x=(-1, 1).
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
- class shenfun.chebyshev.bases.UpperDirichlet(N, quad='GC', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space with single Dirichlet on upper edge
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} - T_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= T_0,\end{split}\]such that
\[\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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of (None, number), optional) – The number is the boundary condition value
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
- class shenfun.chebyshev.bases.UpperDirichletNeumann(N, quad='GC', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]¶
Bases:
CompositeBase
Function space for both Dirichlet and Neumann boundary conditions on the right hand side.
The basis \(\{\phi_k\}_{k=0}^{N-1}\) is
\[\begin{split}\phi_k &= T_{k} - \frac{4(k+1)}{2k+3}T_{k+1} + \frac{2k+1}{2k+3}T_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= T_0, \\ \phi_{N-1} &= -T_0+T_1,\end{split}\]such that
\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(1) &= a, 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, optional) – Number of quadrature points
quad (str, optional) – Type of quadrature
GL - Chebyshev-Gauss-Lobatto
GC - Chebyshev-Gauss
bc (2-tuple of numbers) – Boundary condition values at the right edge of domain
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
shenfun.chebyshev.la module¶
- class shenfun.chebyshev.la.ADDSolver(mats)[source]¶
Bases:
SparseMatrixSolver
- Poisson_Solve_ADD = <shenfun.optimization.runtimeoptimizer object>¶
- solve(b, u, axis, lu)[source]¶
Solve Au=b
Solve along axis if b and u are multidimensional arrays.
- Parameters:
b, u (arrays of rhs and output) – Both can be multidimensional
axis (int) – The axis we are solving over
lu (LU-decomposition) – Can be either the output from splu, or a dia-matrix containing the L and U matrices. The latter is used in subclasses.
- class shenfun.chebyshev.la.ANNSolver(mats)[source]¶
Bases:
SparseMatrixSolver
- __call__(b, u=None, axis=0, constraints=((0, 0),))[source]¶
Solve matrix problem Au = b along axis
This routine also applies boundary conditions and constraints, and performes LU-decomposition on the fully assembled matrix.
- 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 multidimensional
constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val
Note
If u is not provided, then b is overwritten with the solution and returned
- class shenfun.chebyshev.la.Biharmonic(*args)[source]¶
Bases:
object
Multidimensional Biharmonic solver for
\[a_0 u'''' + \alpha u'' + \beta u = b\]where \(u\) is the solution, \(b\) is the right hand side and \(a_0, \alpha\) and \(\beta\) are scalars, or arrays of scalars for a multidimensional problem.
The user must provide mass, stiffness and biharmonic matrices with associated scale arrays \((a_0/\alpha/\beta)\). The matrices and scales can be provided in any order
- Parameters:
S (
TPMatrix
orSpectralMatrix
)A (
TPMatrix
orSpectralMatrix
)B (
TPMatrix
orSpectralMatrix
)scale_S (array, optional)
scale_A (array, optional)
scale_B (array, optional)
If only three arguments are passed, then we decide which matrix is which through inspection. The three scale arrays must then be available as S.scale, A.scale, B.scale. If six arguments are provided they must be in order S, A, B, scale S, scale A, scale B.
Variables are extracted from the matrices
The solver can be used along any axis of a multidimensional problem. For example, if the Chebyshev basis (Biharmonic) is the last in a 3-dimensional TensorProductSpace, where the first two dimensions use Fourier, then the 1D equation listed above arises when one is solving the 3D biharmonic equation
\[\nabla^4 u = b\]With the spectral Galerkin method we multiply this equation with a test function (\(v\)) and integrate (weighted inner product \((\cdot, \cdot)_w\)) over the domain
\[(v, \nabla^4 u)_w = (v, b)_w\]See the Poisson problem for details, since it is actually quite involved. But basically, one obtains a linear algebra system to be solved along the \(z\)-axis for all combinations of the two Fourier indices \(k\) and \(l\)
\[(S_{mj} - 2(k^2 + l^2) A_{mj}) + (k^2 + l^2)^2 B_{mj}) \hat{u}[k, l, j] = (v, b)_w[k, l, m]\]Note that \(k\) only varies along \(x\)-direction, whereas \(l\) varies along \(y\). To allow for Numpy broadcasting these two variables are stored as arrays of shape
\[ \begin{align}\begin{aligned}k : (N, 1, 1)\\l : (1, M, 1)\end{aligned}\end{align} \]Here it is assumed that the solution array \(\hat{u}\) has shape (N, M, P). Now, multiplying \(k\) array with \(\hat{u}\) is achieved as
\[k \cdot \hat{u}\]Numpy will then take care of broadcasting \(k\) to an array of shape (N, M, P) before performing the elementwise multiplication. Likewise, the constant scale \(1\) in front of the \(A_{mj}\) matrix is stored with shape (1, 1, 1), and multiplying with \(\hat{u}\) is performed as if it was a scalar (as it here happens to be).
This is where the scale arrays in the signature to the Helmholt solver comes from. \(a_0\) is here \(1\), whereas \(\alpha\) and \(\beta\) are \(-2(k^2+l^2)\) and \((k^2+l^2)^2\), respectively. Note that \(k+l\) is an array of shape (N, M, 1).
- Biharmonic_Solve = <shenfun.optimization.runtimeoptimizer object>¶
- Biharmonic_factor_pr = <shenfun.optimization.runtimeoptimizer object>¶
- Biharmonic_matvec = <shenfun.optimization.runtimeoptimizer object>¶
- LU_Biharmonic = <shenfun.optimization.runtimeoptimizer object>¶
- class shenfun.chebyshev.la.Helmholtz(*args)[source]¶
Bases:
object
Helmholtz solver
\[\alpha u'' + \beta u = b\]where \(u\) is the solution, \(b\) is the right hand side and \(\alpha\) and \(\beta\) are scalars, or arrays of scalars for a multidimensional problem.
The user must provide mass and stiffness matrices with scale arrays \((\alpha/\beta)\) to each matrix. The matrices and scales can be provided as instances of
TPMatrix
, orSpectralMatrix
.- Parameters:
A (
SpectralMatrix
orTPMatrix
) – mass or stiffness matrixB (
SpectralMatrix
orTPMatrix
) – mass or stiffness matrixscale_A (array, optional) – Scale array to stiffness matrix
scale_B (array, optional) – Scale array to mass matrix
The two matrices must be one stiffness and one mass matrix. Which is which will be found by inspection if only two arguments are provided. The scales \(\alpha\) and \(\beta\) must then be available as A.scale and B.scale. If four arguments are provided they must be in the order
stiffness matrix, mass matrix, scale stiffness, scale mass
The solver can be used along any axis of a multidimensional problem. For example, if the Chebyshev basis (Dirichlet or Neumann) is the last in a 3-dimensional TensorProductSpace, where the first two dimensions use Fourier, then the 1D Helmholtz equation arises when one is solving the 3D Poisson equation
\[\nabla^2 u = b\]With the spectral Galerkin method we multiply this equation with a test function (\(v\)) and integrate (weighted inner product \((\cdot, \cdot)_w\)) over the domain
\[(v, \nabla^2 u)_w = (v, b)_w\]See Poisson’s equation for details, since it is actually quite involved. But basically, one obtains a linear algebra system to be solved along the \(z\)-axis for all combinations of the two Fourier indices \(k\) and \(l\)
\[(A_{mj} - (k^2 + l^2) B_{mj}) \hat{u}[k, l, j] = (v, b)_w[k, l, m]\]Note that \(k\) only varies along \(x\)-direction, whereas \(l\) varies along \(y\). To allow for Numpy broadcasting these two variables are stored as arrays of shape
\[ \begin{align}\begin{aligned}k : (N, 1, 1)\\l : (1, M, 1)\end{aligned}\end{align} \]Here it is assumed that the solution array \(\hat{u}\) has shape (N, M, P). Now, multiplying k array with \(\hat{u}\) is achieved as an elementwise multiplication
\[k \cdot \hat{u}\]Numpy will then take care of broadcasting \(k\) to an array of shape (N, M, P) before performing the elementwise multiplication. Likewise, the constant scale \(1\) in front of the \(A_{mj}\) matrix is stored with shape (1, 1, 1), and multiplying with \(\hat{u}\) is performed as if it was a scalar (as it here happens to be).
This is where the scale arrays come from. \(\alpha\) is here \(1\), whereas \(\beta\) is \((k^2+l^2)\). Note that \(k+l\) is an array of shape (N, M, 1).
- Helmholtz_Neumann_matvec = <shenfun.optimization.runtimeoptimizer object>¶
- Helmholtz_matvec = <shenfun.optimization.runtimeoptimizer object>¶
- LU_Helmholtz = <shenfun.optimization.runtimeoptimizer object>¶
- Solve_Helmholtz = <shenfun.optimization.runtimeoptimizer object>¶
- __call__(b, u=None, constraints=())[source]¶
Solve matrix problem
- Parameters:
b (array) – Array of right hand side on entry and solution on exit unless u is provided.
u (array) – Output array
If b and u are multidimensional, then the axis over which to solve for is determined on creation of the class.
shenfun.chebyshev.matrices module¶
This module contains specific inner product matrices for the different bases in the Chebyshev family.
A naming convention is used for the first capital letter for all matrices. The first letter refers to type of matrix.
Mass matrices start with B
One derivative start with C
Two derivatives (Laplace) start with A
Four derivatives (Biharmonic) start with S
A matrix may consist of different types of test and trialfunctions. The next letters in the matrix name uses the ‘short_name’ method for all these different bases, see chebyshev.bases.py.
So a mass matrix using ShenDirichlet test and ShenNeumann trial is named BSDSNmat.
All matrices in this module may be looked up using the ‘mat’ dictionary, which takes test and trialfunctions along with the number of derivatives to be applied to each. As such the mass matrix BSDSDmat may be looked up as
>>> from shenfun.chebyshev.matrices import mat
>>> from shenfun.chebyshev.bases import ShenDirichlet as SD
>>> B = mat[((SD, 0), (SD, 0))]
and an instance of the matrix can be created as
>>> B0 = SD(10)
>>> BM = B((B0, 0), (B0, 0))
Check that this matrix corresponds to the matrix ‘d’ hardcoded below:
>>> import numpy as np
>>> d = {-2: -np.pi/2,
... 0: np.array([ 1.5*np.pi, np.pi, np.pi, np.pi, np.pi, np.pi, np.pi, np.pi]),
... 2: -np.pi/2}
>>> [np.all(BM[k] == v) for k, v in d.items()]
[True, True, True]
However, this way of creating matrices is not reccommended use. It is far more elegant to use the TrialFunction/TestFunction interface, and to generate the matrix as an inner product:
>>> from shenfun import TrialFunction, TestFunction, inner
>>> u = TrialFunction(B0)
>>> v = TestFunction(B0)
>>> BM = inner(u, v)
>>> [np.all(BM[k] == v) for k, v in d.items()]
[True, True, True]
To see that this is in fact the BSDSDmat:
>>> print(BM.__class__)
<class 'shenfun.chebyshev.matrices.BSDSDmat'>
- class shenfun.chebyshev.matrices.ACNCNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\phi''_j, \phi_k)_w\]where \(\phi_k \in\)
chebyshev.bases.CombinedShenNeumann
, 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.
- class shenfun.chebyshev.matrices.AHHHHmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\phi''_j, \phi_k)_w\]where \(\phi_k \in\)
chebyshev.bases.Heinrichs
, 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.
- class shenfun.chebyshev.matrices.ASBSBmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\phi''_j, \phi_k)_w\]where \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.ASBSDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, \phi_k)_w\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, the trial function \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- class shenfun.chebyshev.matrices.ASBTmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (T''_j, \phi_k)_w\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, the trial function \(T_j \in\)chebyshev.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.
- class shenfun.chebyshev.matrices.ASDHHmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, \phi_k)_w\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial function \(\psi_j \in\)chebyshev.bases.Heinrichs
, 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.
- class shenfun.chebyshev.matrices.ASDMNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, \phi_k)_w\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial function \(\psi_j \in\)chebyshev.bases.MikNeumann
, 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.
- class shenfun.chebyshev.matrices.ASDSDmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\phi''_j, \phi_k)_w\]where \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, 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.
- get_solver()[source]¶
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=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.ASNCNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, \phi_k)_w\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenNeumann
, the trial function \(\psi_j \in\)chebyshev.bases.CombinedShenNeumann
, 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.
- class shenfun.chebyshev.matrices.ASNSNmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\phi''_j, \phi_k)_w\]where \(\phi_k \in\)
chebyshev.bases.ShenNeumann
, 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.
- get_solver()[source]¶
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=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.ATSDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, T_k)_w\]where the test function \(T_k \in\)
chebyshev.bases.Orthogonal
, the trial function \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- class shenfun.chebyshev.matrices.ATSNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (\psi''_j, T_k)_w\]where the test function \(T_k \in\)
chebyshev.bases.Orthogonal
, the trial function \(\psi_j \in\)chebyshev.bases.ShenNeumann
, 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.
- class shenfun.chebyshev.matrices.ATTmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Stiffness matrix \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where
\[a_{kj} = (T''_j, T_k)_w\]where \(T_k \in\)
chebyshev.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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BLDLDmat(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}=(\phi_j, \phi_k)_w,\]where \(\phi_k \in\)
chebyshev.bases.LowerDirichlet
, 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.
- class shenfun.chebyshev.matrices.BSBSBmat(test, trial, scale=1, 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}=(\phi_j, \phi_k)_w,\]where \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, 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.
- get_solver()[source]¶
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=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BSBSDmat(test, trial, scale=1, 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}=(\psi_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BSBTmat(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}=(T_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, the trial \(T_j \in\)chebyshev.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.
- class shenfun.chebyshev.matrices.BSDSDmat(test, trial, scale=1, 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}=(\phi_j, \phi_k)_w,\]where \(\phi_j \in\)
chebyshev.bases.ShenDirichlet
, 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.
- get_solver()[source]¶
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=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BSDSNmat(test, trial, scale=1, 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}=(\psi_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial \(\psi_j \in\)chebyshev.bases.ShenNeumann
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BSDTmat(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}=(T_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial \(\psi_j \in\)chebyshev.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.
- class shenfun.chebyshev.matrices.BSNSDmat(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}=(\psi_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenNeumann
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- class shenfun.chebyshev.matrices.BSNSNmat(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}=(\phi_j, \phi_k)_w,\]where \(\phi_k \in\)
chebyshev.bases.ShenNeumann
, 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.
- get_solver()[source]¶
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=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BTLmat(test, trial, scale=1, 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}=(L_j, T_k)_w,\]where \(T_k \in\)
chebyshev.bases.Orthogonal
, \(L_j \in\)shenfun.legendre.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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.BTSDmat(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}=(\psi_j, T_k)_w,\]where the test function \(T_k \in\)
chebyshev.bases.Orthogonal
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- class shenfun.chebyshev.matrices.BTSNmat(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}=(\psi_j, T_k)_w,\]where the test function \(T_k \in\)
chebyshev.bases.Orthogonal
, the trial \(\psi_j \in\)chebyshev.bases.ShenNeumann
, 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.
- class shenfun.chebyshev.matrices.BTTmat(test, trial, scale=1, 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}=(T_j, T_k)_w,\]where \(T_k \in\)
chebyshev.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.
- matvec(v, c, format='csr', axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- solve(b, u=None, axis=0, constraints=())[source]¶
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 is u is None
u (array, optional) – Solution array if provided
axis (int, optional) – The axis over which to solve for if u is multi- dimensional
constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val
- class shenfun.chebyshev.matrices.CSBSDmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\psi'_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.CSDSBmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\psi'_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial \(\psi_j \in\)chebyshev.bases.ShenBiharmonic
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.CSDSDmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\phi'_j, \phi_k)_w,\]where \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.CSDSNmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\psi'_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial \(\psi_j \in\)chebyshev.bases.ShenNeumann
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.CSDTmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(T'_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenDirichlet
, the trial \(T_j \in\)chebyshev.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.
- class shenfun.chebyshev.matrices.CSNSDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\psi'_j, \phi_k)_w,\]where the test function \(\phi_k \in\)
chebyshev.bases.ShenNeumann
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- class shenfun.chebyshev.matrices.CTSDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(\psi'_j, T_k)_w,\]where the test function \(T_k \in\)
chebyshev.bases.Orthogonal
, the trial \(\psi_j \in\)chebyshev.bases.ShenDirichlet
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.CTTmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Derivative matrix \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where
\[c_{kj}=(T'_j, T_k)_w,\]where \(T_j \in\)
chebyshev.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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
- class shenfun.chebyshev.matrices.SSBSBmat(test, trial, scale=1, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]¶
Bases:
SpectralMatrix
Biharmonic matrix \(S=(s_{kj}) \in \mathbb{R}^{M \times N}\), where
\[s_{kj}=(\phi''''_j, \phi_k)_w,\]where \(\phi_k \in\)
chebyshev.bases.ShenBiharmonic
, 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.
- matvec(v, c, format=None, axis=0)[source]¶
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 Noneaxis (int, optional) – The axis over which to take the matrix vector product
Module contents¶
Functionality for working with Chebyshev bases