shenfun.ultraspherical package

Submodules

shenfun.ultraspherical.bases module

Module for function spaces of ultraspherical type.

The ultraspherical polynomial \(Q^{(\alpha)}_k\) is here defined as

\[Q^{(\alpha)}_k = \frac{1}{P^{(\alpha,\alpha)}_k(1)} P^{(\alpha,\alpha)}_k\]

where \(P^{(\alpha,\alpha)}_k\) is the regular Jacobi polynomial with two equal parameters. The scaling with \((P^{(\alpha,\alpha)}_k(1))^{-1}\) is not standard, but it leads to the boundary values

\[{Q}^{(\alpha)}_k(\pm 1) = (\pm 1)^{k}\]

The Chebyshev (first and second kind) and Legendre polynomials can be defined as

\[\begin{split}T_k(x) &= Q^{(-1/2)}_k(x) \\ U_k(x) &= (k+1)Q^{(1/2)}_k(x) \\ L_k(x) &= Q^{(0)}_k(x)\end{split}\]
class shenfun.ultraspherical.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.ultraspherical.bases.CompactBiharmonic(N, quad='QG', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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{h_n}{b^{(2)}_{n+2,n}} \frac{(1-x^2)^2}{h^{(2,\alpha)}_{k+2}} \frac{d^2Q^{(\alpha)}_{k+2}}{dx^2},\]

where

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

This is Phi2 scaled such that the main diagonal of the stencil matrix is unity.

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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.CompactDirichlet(N, quad='QG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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 &= Q^{(\alpha)}_k - Q^{(\alpha)}_{k+2} \quad k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \tfrac{1}{2}(Q^{(\alpha)}_0 - Q^{(\alpha)}_1) \phi_{N-1} &= \tfrac{1}{2}(Q^{(\alpha)}_0 + Q^{(\alpha)}_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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.CompactNeumann(N, quad='QG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for Neumann boundary conditions

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

\[\]
n*(-2*a - n - 1)/(2*a*n + 4*a + n**2 + 5*n + 6)

phi_k &= Q^{(alpha)}_k - frac{k(2alpha+k+1)}{2 alpha k + 4k + k^2 + 5k + 6}Q^{(alpha)}_{k+2} quad k=0, 1, ldots, N-3, \ phi_{N-2} &= tfrac{1}{2}Q^{(alpha)}_1 - frac{(alpha+1)}{2(2alpha+3) Q^{(alpha)}_2} phi_{N-1} &= tfrac{1}{2}Q^{(alpha)}_1 + frac{(alpha+1)}{2(2alpha+3) Q^{(alpha)}_2}

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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Generic(N, quad='QG', bc={}, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, alpha=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 - QG - 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.

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

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

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

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

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

Bases: CompositeBase

Function space with single Dirichlet on left edge

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

\[\begin{split}\phi_k &= Q^{(\alpha)}_{k} + Q^{(\alpha)}_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= Q^{(\alpha)}_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) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • QG - Jacobi-Gauss

  • bc (2-tuple of (None, number), optional) – The number is the boundary condition value

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Orthogonal(N, quad='QG', alpha=0, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)[source]

Bases: JacobiBase

Function space for regular (orthogonal) ultraspherical polynomials

The orthogonal basis is

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

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

Parameters:

  • N (int) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • QG - Jacobi-Gauss

  • alpha (number, optional) – Parameter of the ultraspherical 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

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

derivative_jacobiQ(x, alpha, 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_bc_space()[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
static jacobiQ(x, alpha, N)[source]
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)

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]
sympy_stencil(i=i, j=j)[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.ultraspherical.bases.Phi1(N, quad='QG', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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)}_{k+1}} \frac{dQ^{(\alpha)}_{k+1}}{dx}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \tfrac{1}{2}(Q^{(\alpha)}_0 - Q^{(\alpha)}_1) \phi_{N-1} &= \tfrac{1}{2}(Q^{(\alpha)}_0 + Q^{(\alpha)}_1)\end{split}\]

where

\[\begin{split}h^{(1,\alpha)}_k&=\int_{-1}^1 \left(\frac{dQ^{(\alpha)}_{k}}{dx}\right)^2 (1-x^2)^{\alpha+1}dx, \\\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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Phi2(N, quad='QG', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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)}_{k+2}} \frac{d^2Q^{(\alpha)}_{k+2}}{dx^2},\]

where

\[\begin{split}h^{(2,\alpha)}_k&=\int_{-1}^1 \left(\frac{d^2Q^{(\alpha)}_{k}}{dx^2}\right)^2 (1-x^2)^{\alpha+2}dx, \\ &= \frac{2^{2 \alpha + 1} \cdot \left(2 \alpha + n + 1\right) \left(2 \alpha + n + 2\right) \Gamma^{2}\left(\alpha + n + 1\right)}{\left(2 \alpha + 2 n + 1\right) \Gamma\left(n - 1\right) \Gamma\left(2 \alpha + n + 1\right)},\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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Phi3(N, quad='QG', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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)}_{k+3}} \frac{d^3Q^{(\alpha)}_{k+3}}{dx^3}, \\\end{split}\]

where

\[\begin{split}h^{(3,\alpha)}_k=\int_{-1}^1 \left(\frac{d^3Q^{(\alpha)}_{k}}{dx^3}\right)^2 (1-x^2)^{\alpha+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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Phi4(N, quad='QG', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, alpha=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)}_{k+4}} \frac{d^4Q^{(\alpha)}_{k+4}}{dx^4}, \\\end{split}\]

where

\[\begin{split}h^{(4,\alpha)}_k&=\int_{-1}^1 \left(\frac{d^4Q^{(\alpha)}_{k}}{dx^4}\right)^2 (1-x^2)^{\alpha+4}dx, \\ &=\frac{2^{2 \alpha + 1} \cdot \left(2 \alpha + n + 1\right) \left(2 \alpha + n + 2\right) \left(2 \alpha + n + 3\right) \left(2 \alpha + n + 4\right) \Gamma^{2}\left(\alpha + n + 1\right)}{\left(2 \alpha + 2 n + 1\right) \Gamma\left(n - 3\right) \Gamma\left(2 \alpha + 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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.Phi6(N, quad='QG', 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, alpha=0, coordinates=None, **kw)[source]

Bases: CompositeBase

Function space for 12th order equations

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

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

where

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

The 12 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^{(k)}(-1)&=a_k, u^{(k)}(1)=b_k, k=0, 1, \ldots, 5\end{split}\]

The last 12 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

    • QG - Jacobi-Gauss

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

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.bases.UpperDirichlet(N, quad='QG', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, alpha=0, **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 &= Q^{(\alpha)}_{k} - Q^{(\alpha)}_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= Q^{(\alpha)}_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) – Number of quadrature points

  • quad (str, optional) – Type of quadrature

    • QG - Jacobi-Gauss

  • bc (2-tuple of (None, number), optional) – The number is the boundary condition value

  • alpha (number, optional) – Parameter of the ultraspherical 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.ultraspherical.matrices module

class shenfun.ultraspherical.matrices.BQQmat(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}=(Q_j, Q_k)_w,\]

\(Q_k \in\) ultraspherical.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.

Module contents