shenfun.jacobi package

Submodules

shenfun.jacobi.bases module

Module for function spaces of generalized Jacobi type

class shenfun.jacobi.bases.BCGeneric(N, bc=None, domain=None, alpha=0, beta=0, **kw)

Bases: CompositeBase

Function space for setting inhomogeneous boundary conditions

Parameters:

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

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

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

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

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

basis_function(i=0, x=x)

Return basis function i

Parameters:

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

  • x (sympy Symbol, optional)

static boundary_condition()
property dim_ortho
eval(x, u, output_array=None)

Evaluate Function u at position x

Parameters:

  • x (float or array of floats)

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

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

Returns:

output_array

Return type:

array

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

Evaluate basis function i at points x

Parameters:

  • x (float or array of floats)

  • i (int, optional) – Basis function number

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

Returns:

output_array

Return type:

array

evaluate_basis_derivative(x=None, i=0, k=0, output_array=None)

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

Parameters:

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

  • i (int, optional) – Basis function number

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

  • output_array (array, optional) – return array

Returns:

output_array

Return type:

array

get_orthogonal(**kwargs)

Return orthogonal space (otherwise as self)

Returns:

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

Return type:

SpectralBase

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 of x and M the number of bases.

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

x (array of floats) – points for evaluation

Note

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

class shenfun.jacobi.bases.CompactDirichlet(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Dirichlet boundary conditions

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

\[\begin{split}\phi_k &= P^{(\alpha,\beta)}_k + a_k P^{(\alpha,\beta)}_{k+1} + b_k P^{(\alpha,\beta)}_{k+2} \quad k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]

where

\[\begin{split}a_k &= \frac{\left(\alpha - \beta\right) \left(n + 1\right) \left(\alpha + \beta + 2 n + 3\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + 2 n + 4\right)}, \\ b_k &= - \frac{\left(n + 1\right) \left(n + 2\right) \left(\alpha + \beta + 2 n + 2\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + 2 n + 4\right)}, \\\end{split}\]

and the expansion is

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

The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.CompactNeumann(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Neumann boundary conditions

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

\[\begin{split}\phi_k &= P^{(\alpha,\beta)}_k + a_k P^{(\alpha,\beta)}_{k+1} + b_k P^{(\alpha,\beta)}_{k+2} \quad k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]

where

\[\begin{split}a_k &= \frac{n \left(\alpha - \beta\right) \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + 2 n + 3\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + n + 2\right) \left(\alpha + \beta + 2 n + 4\right)}, \\ b_k &= - \frac{n \left(n + 1\right) \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + 2 n + 2\right)}{\left(\alpha + n + 1\right) \left(\beta + n + 1\right) \left(\alpha + \beta + n + 3\right) \left(\alpha + \beta + 2 n + 4\right)}, \\\end{split}\]

and the expansion is

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

The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.Generic(N, quad='JG', bc={}, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, alpha=0, beta=0, **kw)[source]

Bases: CompositeBase

Function space for space with any boundary conditions

Any combination of Dirichlet and Neumann is possible.

Parameters:

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

  • quad (str, optional) – Type of quadrature - JG - Jacobi-Gauss

  • bc (dict, optional) – The dictionary must have keys ‘left’ and ‘right’, to describe boundary conditions on the left and right boundaries, and a list of 2-tuples to specify the condition. Specify Dirichlet on both ends with

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

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

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

    For both conditions on the right do

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

    Any combination should be possible, and it should also be possible to use second derivatives N2. See BoundaryConditions.

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

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

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

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

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

Note

A test function is always using homogeneous boundary conditions.

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.JacobiBase(N, quad='', alpha=0, beta=0, padding_factor=1, domain=(-1.0, 1.0), dtype=None, dealias_direct=False, coordinates=None)[source]

Bases: SpectralBase

Abstract base class for Jacobi function spaces

L2_norm_sq(i)[source]

Return square of L2-norm

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

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

Parameters:

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

static bnd_values(k=0, alpha=None, beta=None, gn=None)[source]
get_bc_space()[source]
property is_jacobi
l2_norm_sq(i=None)[source]

Return square of l2-norm

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

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

Parameters:

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

lagrange_poly(i)[source]

Return i’th Lagrange polynomial of self

Parameters:

i (int) – For the i’th Lagrange polynomial

reference_domain()[source]
unweighted_quadrature_weights()[source]

Return quadrature weights for unweighted integrals

\[\int_{-1}^{1} f dx \approx \sum_{j=0}^{N} f(x_j) w_j,\]

where w_j are the quadrature weights.

class shenfun.jacobi.bases.LowerDirichlet(N, quad='JG', bc=(0, None), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Dirichlet boundary condition at left edge

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

\[\begin{split}\phi_k &= P^{(\alpha, \beta)}_k + \frac{(k+1)}{\beta+k+1} P^{(\alpha,\beta)}_{k+1} \quad k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= 1\end{split}\]

and the expansion is

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

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

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (2-tuple of (number, None), optional) – Boundary condition at x=-1.

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.Orthogonal(N, quad='JG', alpha=0, beta=0, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]

Bases: JacobiBase

Function space for regular (orthogonal) Jacobi functions

The orthogonal basis is

\[\phi_k = P^{(\alpha,\beta)}_k, \quad k = 0, 1, \ldots, N-1,\]

where \(P^{(\alpha,\beta)}_k\) is the Jacobi polynomial. The basis is orthogonal with weight \((1-x)^{\alpha}(1+x)^{\beta}\).

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

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

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

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

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

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

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

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

derivative_jacobi(x, alpha, beta, k=1)[source]
eval(x, u, output_array=None)[source]

Evaluate Function u at position x

Parameters:

  • x (float or array of floats)

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

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

Returns:

output_array

Return type:

array

evaluate_basis(x, i=0, output_array=None)[source]

Evaluate basis function i at points x

Parameters:

  • x (float or array of floats)

  • i (int, optional) – Basis function number

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

Returns:

output_array

Return type:

array

evaluate_basis_all(x=None, argument=0)[source]

Evaluate basis at x or all quadrature points

Parameters:

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

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

Returns:

Vandermonde matrix

Return type:

array

evaluate_basis_derivative(x=None, i=0, k=0, output_array=None)[source]

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

Parameters:

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

  • i (int, optional) – Basis function number

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

  • output_array (array, optional) – return array

Returns:

output_array

Return type:

array

evaluate_basis_derivative_all(x=None, k=0, argument=0)[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

static family()[source]
get_orthogonal(**kwargs)[source]

Return orthogonal space (otherwise as self)

Returns:

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

Return type:

SpectralBase

property is_orthogonal
jacobi(x, alpha, beta, N)[source]
orthogonal_basis_function(i=0, x=x)[source]

Return the orthogonal basis function i

Parameters:

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

  • x (sympy Symbol, optional)

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

Return points and weights of quadrature for weighted integral

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

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

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

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

Note

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

static short_name()[source]
stencil_matrix(N=None)[source]
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 of x and M the number of bases.

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

x (array of floats) – points for evaluation

Note

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

class shenfun.jacobi.bases.Phi1(N, quad='JG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Dirichlet boundary conditions

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

\[\begin{split}\phi_k &= \frac{(1-x^2)}{h^{(1,\alpha,\beta)}_{k+1}} \frac{dP^{(\alpha,\beta)}_{k+1}}{dx}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \frac{\alpha + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0 - \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1, \\ \phi_{N-1} &= \frac{\beta + 1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_0+ \frac{1}{\alpha + \beta + 2}P^{(\alpha,\beta)}_1,\end{split}\]

where

\[\begin{split}h^{(1,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{dP^{(\alpha,\beta)}_{k}}{dx}\right)^2 (1-x)^{\alpha+1}(1+x)^{\beta+1}dx, \\ &= \frac{2^{\alpha+\beta+1}(\alpha+\beta+k+1) \Gamma{(\alpha+k+1)} \Gamma{(\beta+k+1)}}{(\alpha+\beta+2k+1) \Gamma{(k)} \Gamma{(\alpha+\beta+k+1)}},\end{split}\]

and the expansion is

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

The last two bases are for boundary conditions and only used if a or b are different from 0. In one dimension \(\hat{u}_{N-2}=a\) and \(\hat{u}_{N-1}=b\).

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (2-tuple of numbers, optional) – Boundary conditions at, respectively, x=(-1, 1).

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.Phi2(N, quad='JG', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for biharmonic equations

The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-5\) are

\[\phi_k &= \frac{(1-x^2)^2}{h^{(2,\alpha,\beta)}_{k+2}} \frac{d^2P^{(\alpha,\beta)}_{k+2}}{dx^2},\]

where

\[\begin{split}h^{(2,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^2P^{(\alpha,\beta)}_{k}}{dx^2}\right)^2 (1-x)^{\alpha+2}(1+x)^{\beta+2}dx, \\ &= \frac{2^{\alpha+\beta+1}(\alpha+\beta+k+1) \Gamma{(\alpha+k+1)} \Gamma{(\beta+k+1)}}{(\alpha+\beta+2k+1) \Gamma{(k)} \Gamma{(\alpha+\beta+k+1)}},\end{split}\]

The 4 boundary basis functions are computed using jacobi.findbasis.get_bc_basis(), but they are too messy to print here. We have

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

The last 4 basis functions are for boundary conditions and only used if a, b, c or d are different from 0. In one dimension \(\hat{u}_{N-4}=a\), \(\hat{u}_{N-3}=b\), \(\hat{u}_{N-2}=c\) and \(\hat{u}_{N-1}=d\).

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (4-tuple of numbers, optional) – Boundary conditions.

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.Phi3(N, quad='JG', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for 6th order equations

The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-7\) are

\[\begin{split}\phi_k &= \frac{(1-x^2)^3}{h^{(3,\alpha,\beta)}_{k+3}} \frac{d^3P^{(\alpha,\beta)}_{k+3}}{dx^3}, \\\end{split}\]

where

\[\begin{split}h^{(3,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^3P^{(\alpha,\beta)}_{k}}{dx^3}\right)^2 (1-x)^{\alpha+3}(1+x)^{\beta+3}dx, \\\end{split}\]

The 6 boundary basis functions are computed using jacobi.findbasis.get_bc_basis(), but they are too messy to print here. We have

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

The last 6 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (6-tuple of numbers, optional) – Boundary conditions.

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.Phi4(N, quad='JG', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for 8th order equations

The basis functions \(\phi_k\) for \(k=0, 1, \ldots, N-9\) are

\[\begin{split}\phi_k &= \frac{(1-x^2)^4}{h^{(4,\alpha,\beta)}_{k+4}} \frac{d^4P^{(\alpha,\beta)}_{k+4}}{dx^4}, \\\end{split}\]

where

\[\begin{split}h^{(4,\alpha,\beta)}_k&=\int_{-1}^1 \left(\frac{d^4P^{(\alpha,\beta)}_{k}}{dx^4}\right)^2 (1-x)^{\alpha+4}(1+x)^{\beta+4}dx, \\ &=\frac{2^{\alpha + \beta + 1} \left(\alpha + \beta + n + 1\right) \left(\alpha + \beta + n + 2\right) \left(\alpha + \beta + n + 3\right) \left(\alpha + \beta + n + 4\right) \Gamma\left(\alpha + n + 1\right) \Gamma\left(\beta + n + 1\right)}{\left(\alpha + \beta + 2 n + 1\right) \Gamma\left(n - 3\right) \Gamma\left(\alpha + \beta + n + 1\right)},\end{split}\]

The 8 boundary basis functions are computed using jacobi.findbasis.get_bc_basis(), but they are too messy to print here. We have

\[\begin{split}u(x) &= \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \\ u(-1) &= a, u'(-1)=b, u''(-1)=c, u''''(-1)=d, u(1)=e, u'(1)=f, u''(1)=g, u''''(1)=h.\end{split}\]

The last 8 basis functions are for boundary conditions and only used if there are nonzero boundary conditions.

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (8-tuple of numbers, optional) – Boundary conditions.

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]
class shenfun.jacobi.bases.UpperDirichlet(N, quad='JG', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, beta=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Dirichlet boundary conditions at right edge

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

\[\begin{split}\phi_k &= P^{(\alpha, \beta)}_k - \frac{(k+1)}{\alpha+k+1} P^{(\alpha,\beta)}_{k+1} \quad k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= 1\end{split}\]

and the expansion is

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

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

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • JG - Jacobi-Gauss

  • bc (2-tuple of (None, number), optional) – Boundary condition at x=1.

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

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

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

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

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

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

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

static boundary_condition()[source]
static short_name()[source]

shenfun.jacobi.matrices module

class shenfun.jacobi.matrices.BJJmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]

Bases: SpectralMatrix

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

\[b_{kj}=(J_j, J_k)_w,\]

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

assemble(method)[source]

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

  • ‘exact’

  • ‘quadrature’

Note

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

  • legendre.matrix

  • chebyshev.matrix

  • chebyshevu.matrix

  • ultraspherical.matrix

  • fourier.matrix

  • laguerre.matrix

  • hermite.matrix

  • jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

shenfun.jacobi.recursions module

Recursions for Jacobi polynomials, or standardized Jacobi polynomials

shenfun.jacobi.recursions.Lmat(k, q, l, M, N, alf=0, bet=0, gn=1)[source]

Return matrix corresponding to

\[(\partial^{k-l}Q_{n}, x^q \phi^{(k)}_m)_{\omega}\quad (1)\]

where \(\partial^k\) represents the \(k\)’th derivative and

\[\begin{split}Q_n = g_n^{(\alpha, \beta)} P^{(\alpha, \beta)}_n \\ \phi^{(k)}_m = \frac{(1-x^2)^k \partial^k Q_{m+k}}{h^{(k)}_{m+k}}\end{split}\]
Parameters:

  • k, q, l (integers) – Numbers in variational form (1)

  • M, N (int) – Shape of matrix

  • alf, bet (numbers, optional) – Jacobian parameters

  • gn (scaling function) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

Example

>>> from shenfun.jacobi.recursions import Lmat, cn, half
>>> X = Lmat(2, 2, 0, 10, 12, -half, -half, cn)
shenfun.jacobi.recursions.ShiftedMatrix(mat, q, r, s, M=0, N=0, k=None, alf=0, bet=0, gn=1)[source]

Index-shifted q’th power of matrix mat. Either

\[\begin{split}A^{(k, q)}_{(r, s)} = (a^{(k, q)}_{m+r,n+s})_{m,n=0}^{M,N} \\ B^{(q)}_{(r, s)} = (b^{(q)}_{m+r,n+s})_{m,n=0}^{M,N}\end{split}\]
Parameters:

  • mat (Python function or SparseMatrix) – The Python function must have signature (k, q, alf, bet, i, j, gn=1) or (alf, bet, i, j, gn=1) if k is None

  • q (int) – The matrix power

  • r (int) – Shift in row

  • s (int) – Shift in column

  • M, N (int, optional) – The shape of the final matrix

  • k (int or None, optional) – Parameter of recursion of the k’th derivative This is only used if the recursion matrix is A.

  • alf, bet (numbers, optional) – Jacobi paramters

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

Note

If mat is a SparseMatrix, then the resulting shifted matrix will be smaller than mat.

shenfun.jacobi.recursions.a(alf, bet, i, j, gn=1)[source]

Recursion matrix \(A\) for standardized Jacobi polynomials

The recursion is

\[x \boldsymbol{Q} = {A}^T \boldsymbol{Q}\]

where

\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • i, j (int) – Indices for row and column

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.a_(k, q, alf, bet, i, j, gn=1)[source]

Recursion matrix \(\underline{A}\) for standardized Jacobi polynomials

The recursion is

\[x \partial^k \boldsymbol{Q} = \underline{A}^T \partial^k \boldsymbol{Q} \quad (*)\]

where \(\partial^k\) represents the \(k\)’th derivative and

\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]
Parameters:

  • k, q (int) – Parameters in (*)

  • alf, bet (numbers) – Jacobi parameters

  • i, j (int) – Indices for row and column

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.a_mat(mat, k, q, alf, bet, M, N, gn=1)[source]

Return SparseMatrix of recursion matrix

\[A^{(k,q)} = (a^{(k,q)}_{mn})_{m,n=0}^{M, N}\]
Parameters:

  • mat (Python function for matrix) – The Python function must have signature (k, q, alf, bet, i, j, gn=1)

  • k (int) – Parameter

  • q (int) – Matrix power

  • alf, bet (numbers) – Jacobi parameters

  • M, N (int) – Shape of returned matrix

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.b(alf, bet, i, j, gn=1)[source]

Recursion matrix \(B\) for standardized Jacobi polynomials

The recursion is

\[\boldsymbol{Q} = {B}^T \partial \boldsymbol{Q}\]

where \(\partial\) represents the derivative and

\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • i, j (int) – Indices for row and column

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.bnd_values(alf, bet, k=0, gn=1)[source]

Return lambda function for computing boundary values

\[\frac{d^k}{dx^k}Q_n(\pm 1)\]

where \(Q^{(\alpha,\beta)}_n = g^{(\alpha,\beta)}_n P_n^{(\alpha,\beta)}\).

Parameters:

  • alf, bet (Numbers) – Jacobi parameters

  • k (int, optional) – Number of derivatives

  • gn (scaling function, optional)

shenfun.jacobi.recursions.c(alf, bet, i, j, gn=1)[source]

Recursion matrix \(C\) for standardized Jacobi polynomials

The recursion is

\[(1-x^2) \partial \boldsymbol{Q} = {C}^T \boldsymbol{Q}\]

where \(\partial\) represents the derivative and

\[\begin{split}Q_n(x) = g_n^{(\alpha,\beta)} P^{(\alpha,\beta)}_n(x) \\ \boldsymbol{Q} = (Q_0, Q_1, \ldots, Q_{N})^T\end{split}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • i, j (int) – Indices for row and column

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.cn(alf, bet, n)[source]

Scaling function

\[c_n^{(\alpha, \beta)} = \frac{1}{P^{(\alpha,\beta)}_n(1)}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

Note

Used by Legendre and Chebyshev polynomials of first kind.

shenfun.jacobi.recursions.djt(alf, bet, n, k)[source]

k’th derivative of Jacobi polynomial

Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

  • k (int) – The k’th derivative

shenfun.jacobi.recursions.dqn(alf, bet, n, k, gn=<function cn>)[source]

Derivative of specialized Jacobi polynomial

\[\frac{d^k Q_n^{(\alpha, \beta)}}{dx^k} = g_n^{(\alpha, \beta)} \frac{d^k P^{(\alpha, \beta)}_n}{dx^k}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

  • k (int) – The k’th derivative

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.gamma(alf, bet, n)[source]

Return normalization factor \(h_n\) for inner product of Jacobi polynomials

\[h_n = (P^{(\alpha,\beta)}_n, P^{(\alpha,\beta)}_n)_{\omega^{(\alpha,\beta)}}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

shenfun.jacobi.recursions.h(alf, bet, n, k, gn=1)[source]

Return normalization factor \(h^{(k)}_n\) for inner product of derivatives of Jacobi polynomials

\[\begin{split}Q_n(x) = g_n(x)P^{(\alpha,\beta)}_n(x) \\ h_n^{(k)} = (\partial^k Q_n, \partial^k Q_n)_{\omega^{(\alpha+k,\beta+k)}} \quad (*)\end{split}\]

where \(\partial^k\) represents the \(k\)’th derivative.

Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

  • k (int) – For derivative of k’th order, see (*)

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.jt(alf, bet, n)[source]

Jacobi polynomial

Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

shenfun.jacobi.recursions.matpow(mat, q, alf, bet, i, j, gn=1)[source]

Compute and return component of q’th matrix power of mat

Parameters:

  • mat (Python function (a, b or c))

  • q (int) – matrix power

  • alf, bet (numbers) – Jacobi parameters

  • i, j (int) – Row and column indices

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.pmat(mat, q, alf, bet, M, N, gn=1)[source]

Return SparseMatrix of q’th matrix power of recursion matrix mat

Parameters:

  • mat (Python function for matrix) – The Python function must have signature (alf, bet, i, j, gn=1)

  • q (int) – Matrix power

  • alf, bet (numbers) – Jacobi parameters

  • M, N (int) – Shape of returned matrix

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.psi(alf, bet, n, k)[source]

Normalization factor for

\[\partial^k P^{(\alpha, \beta)}_n = \psi^{(k,\alpha,\beta)}_{n} P^{(\alpha+k,\beta+k)}_{n-k}, \quad n \ge k, \quad (*)\]

where \(\partial^k\) represents the \(k\)’th derivative

Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n, k (int) – Parameters in (*)

shenfun.jacobi.recursions.qn(alf, bet, n, gn=<function cn>)[source]

Specialized Jacobi polynomial

\[Q_n^{(\alpha, \beta)} = g_n^{(\alpha, \beta)} P^{(\alpha, \beta)}_n\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

  • gn (scaling function, optional) – Chebyshev of first and second kind use cn and un, respectively. Legendre uses gn=1.

shenfun.jacobi.recursions.un(alf, bet, n)[source]

Scaling function

\[u_n^{(\alpha, \beta)} = \frac{n+1}{P^{(\alpha,\beta)}_n(1)}\]
Parameters:

  • alf, bet (numbers) – Jacobi parameters

  • n (int) – Index

Note

Used by Chebyshev polynomials of second kind

Module contents