shenfun.chebyshevu package

Submodules

shenfun.chebyshevu.bases module

Module for defining function spaces in the Chebyshev family

class shenfun.chebyshevu.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.chebyshevu.bases.Compact3(N, quad='GU', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space for 6’th order equation

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

\[\phi_k &= U_k - \frac{3(n+1)}{n+5}U_{k+2} + \frac{3(n+1)(n+2)}{(n+5)(n+6)}U_{k+4} - \frac{(n+1)(n+2)(n+3)}{(n+5)(n+6)(n+7)}U_{k+6}\]

This is the same basis as Phi3, only scaled such that the main diagonal of the stencil matrix contains ones.

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 only used if there are nonzero boundary conditions.

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (6-tuple of floats, optional) – 3 boundary conditions at, respectively, 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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.CompactBiharmonic(N, quad='GU', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, scaled=False, **kw)

Bases: CompositeBase

Function space for biharmonic equation.

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

\[\begin{split}\phi_k &= \frac{h_k}{b^{(2)}_{k+2,k}}\frac{(1-x^2)^2 U''_{k+2}}{h^{(2)}_{k+2}} \\ &= U_k - \frac{2k+2}{k+4}U_{k+2} + \frac{(k+1)(k+2)}{(k+4)(k+5)}U_{k+4}\end{split}\]

where \(h^{(2)}_k = \frac{\pi (k+3)(k+2)k(k-1)}{2}\), \(h_k=\pi/2\) and \(b^{(2)}_{k+2, k}= \frac{1}{4(k+1)(k+2)}\). The 4 boundary functions are

\[\begin{split}\phi_{N-4} &= \tfrac{1}{2}U_0-\tfrac{5}{6}U_1+\tfrac{1}{32}U_3, \\ \phi_{N-3} &= \tfrac{3}{16}U_0-\tfrac{1}{16}U_1-\tfrac{1}{16}U_2+\tfrac{1}{32}U_3, \\ \phi_{N-2} &= \tfrac{1}{2}U_0+\tfrac{5}{16}U_1-\tfrac{1}{32}U_3), \\ \phi_{N-1} &= -\tfrac{3}{16}U_0-\tfrac{1}{16}U_1+\tfrac{1}{16}U_2+\tfrac{1}{32}U_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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (4-tuple of floats, optional) – 2 boundary conditions at, respectively, 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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.CompactDirichlet(N, quad='GU', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, scaled=False, **kw)

Bases: CompositeBase

Function space for Dirichlet boundary conditions.

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

\[\begin{split}\phi_k &= U_k - \frac{k+1}{k+3}U_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \tfrac{1}{2}U_0 - \tfrac{1}{4}U_1, \\ \phi_{N-1} &= \tfrac{1}{2}U_0 + \tfrac{1}{4}U_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 the first N-2 basis functions are

\[\begin{split}\phi_k &= \frac{U_k}{k+1} - \frac{U_{k+2}}{k+3}, \, k=0, 1, \ldots, N-3, \\\end{split}\]

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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.

• scaled (boolean, optional) – Whether or not to use scaled basis function.

• 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

This basis function is a scaled version of Phi1.

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.CompactNeumann(N, quad='GU', bc=(0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space for Neumann boundary conditions

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

\[\begin{split}\phi_k = {U_k}-\frac{k(k+1)}{(k+3)(k+4)} U_{k+2}, \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \tfrac{1}{16}(4U_1-U_2), \\ \phi_{N-1} &= \tfrac{1}{16}(4U_1+U_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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Generic(N, quad='GU', bc={}, domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.LowerDirichlet(N, quad='GU', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space with single Dirichlet on left edge

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

\[\begin{split}\phi_k &= U_{k} + \frac{k+1}{k+2} U_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= U_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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (2-tuple of (number, None), optional) – Boundary condition at 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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Orthogonal(N, quad='GU', domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: JacobiBase

Function space for Chebyshev series of second kind

The orthogonal basis is

\[U_k, \quad k = 0, 1, \ldots, N-1,\]

where \(U_k\) is the \(k\)’th Chebyshev polynomial of the second kind.

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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)

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

basis_function(i=0, x=x)

Return basis function i

Parameters:

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

• x (sympy Symbol, optional)

eval(x, u, output_array=None)

Evaluate Function u at position x

Parameters:

• x (float or array of floats)

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

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

Returns:

output_array

Return type:

array

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

Evaluate basis function i at points x

Parameters:

• x (float or array of floats)

• i (int, optional) – Basis function number

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

Returns:

output_array

Return type:

array

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

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

Parameters:

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

• i (int, optional) – Basis function number

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

• output_array (array, optional) – return array

Returns:

output_array

Return type:

array

evaluate_basis_derivative_all(x=None, k=0, argument=0)

Return k’th derivative of basis evaluated at x or all quadrature points as a Vandermonde matrix.

Parameters:

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

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

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

Returns:

Vandermonde matrix

Return type:

array

static family()

get_bc_space()

get_orthogonal(**kwargs)

Return orthogonal space (otherwise as self)

Returns:

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

Return type:

SpectralBase

property is_orthogonal

l2_norm_sq(i=None)

Return square of l2-norm

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

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

Parameters:

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

orthogonal_basis_function(i=0, x=x)

Return the orthogonal basis function i

Parameters:

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

• x (sympy Symbol, optional)

plan(shape, axis, dtype, options)

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)

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

stencil_matrix(N=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().

weight(x=x)

Weight of inner product space

Parameters:

x (coordinate)

class shenfun.chebyshevu.bases.Phi1(N, quad='GU', bc=(0.0, 0.0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, scaled=False, **kw)

Bases: CompositeBase

Function space for Dirichlet boundary conditions

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

\[\begin{split}\phi_k &= \frac{1}{\pi}\left(\frac{U_k}{k+1}-\frac{U_{k+2}}{k+3}\right) = \frac{(1-x^2)U'_{k+1}}{h^{(1)}_{k+1}} \, k=0, 1, \ldots, N-3, \\ \phi_{N-2} &= \tfrac{1}{2}U_0 - \tfrac{1}{4}U_1, \\ \phi_{N-1} &= \tfrac{1}{2}U_0 + \tfrac{1}{4}U_1,\end{split}\]

where \(h^{(1)}_n = \frac{\pi n(n+2)}{2}\). We have

\[\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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

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

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

• scaled (boolean, optional) – Whether or not to scale basis function n by 1/(n+2)

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

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

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

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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Phi2(N, quad='GU', bc=(0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, scaled=False, **kw)

Bases: CompositeBase

Function space for biharmonic equation.

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

\[\begin{split}\phi_k &= \frac{(1-x^2)^2 U''_{k+2}}{h^{(2)}_{k+2}} \\ &= \frac{1}{2\pi(k+1)(k+2)}\left(U_k- \frac{2(k+1)}{k+4}U_{k+2} + \frac{(k+1)(k+2)}{(k+4)(k+5)}U_{k+4} \right)\end{split}\]

where \(h^{(2)}_n = \frac{\pi (n+3)(n+2)n(n-1)}{2}\). The 4 boundary functions are

\[\begin{split}\phi_{N-4} &= \tfrac{1}{2}U_0-\tfrac{5}{6}U_1+\tfrac{1}{32}U_3, \\ \phi_{N-3} &= \tfrac{3}{16}U_0-\tfrac{1}{16}U_1-\tfrac{1}{16}U_2+\tfrac{1}{32}U_3, \\ \phi_{N-2} &= \tfrac{1}{2}U_0+\tfrac{5}{16}U_1-\tfrac{1}{32}U_3), \\ \phi_{N-1} &= -\tfrac{3}{16}U_0-\tfrac{1}{16}U_1+\tfrac{1}{16}U_2+\tfrac{1}{32}U_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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (4-tuple of floats, optional) – 2 boundary conditions at, respectively, 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.

• scaled (boolean, optional) – Whether or not to scale basis function n by 1/(n+3)

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

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

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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Phi3(N, quad='GU', bc=(0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space for 6’th order equation

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

\[\begin{split}\phi_k &= \frac{(1-x^2)^3}{h^{(3)}_{k+3}} U^{(3)}_{k+3} \\ h^{(3)}_{k+3} &= \frac{\pi \Gamma (k+8)}{2(k+4)k!} = \int_{-1}^1 U^{(3)}_k U^{(3)}_k (1-x^2)^{3.5}} dx.\end{split}\]

where \(U^{(3)}_k\) is the 3rd derivative of \(U_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 only used if there are nonzero boundary conditions.

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (6-tuple of floats, optional) – 3 boundary conditions at, respectively, 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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Phi4(N, quad='GU', bc=(0, 0, 0, 0, 0, 0, 0, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space for 8th order equation

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

\[\begin{split}\phi_k &= \frac{(1-x^2)^4}{h^{(4)}_{k+4}} U^{(4)}_{k+4} \\ h^{(4)}_{k+4} &= \frac{\pi \Gamma (k+10)}{2(k+5)k!} = \int_{-1}^1 U^{(4)}_{k+4} U^{(4)}_{k+4} (1-x^2)^{4.5}} dx.\end{split}\]

where \(U^{(4)}_k\) is the 4th derivative of \(U_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, u''(1)=g, u'''(1)=h\end{split}\]

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

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.Phi6(N, quad='GU', 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)

Bases: CompositeBase

Function space for 12th order equation

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

\[\begin{split}\phi_k &= \frac{(1-x^2)^6}{h^{(6)}_{k+6}} U^{(6)}_{k+6} \\ h^{(6)}_k &= \frac{\pi (k+7)!}{2(k-6)!} = \int_{-1}^1 U^{(6)}_k U^{(6)}_k (1-x^2)^{6.5} dx,\end{split}\]

where \(U^{(6)}_k\) is the 6th derivative of \(U_k\). The 12 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^{(k)}(-1)&=a_k, u^{(k)}(1)=b_k, k=0, 1, \ldots, 5\end{split}\]

The last 12 basis functions are only used if there are nonzero boundary conditions.

Parameters:

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

• quad (str, optional) – Type of quadrature

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-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

static boundary_condition()

static short_name()

class shenfun.chebyshevu.bases.UpperDirichlet(N, quad='GU', bc=(None, 0), domain=(-1, 1), dtype=<class 'float'>, padding_factor=1, dealias_direct=False, coordinates=None, **kw)

Bases: CompositeBase

Function space with single Dirichlet on upper edge

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

\[\begin{split}\phi_k &= U_{k} - \frac{k+1}{k+2} U_{k+1}, \, k=0, 1, \ldots, N-2, \\ \phi_{N-1} &= U_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

  • GC - Chebyshev-Gauss

  • GU - Chebyshevu-Gauss

• bc (2-tuple of (None, number), optional) – Boundary condition at 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

static boundary_condition()

static short_name()

shenfun.chebyshevu.matrices module

This module contains specific inner product matrices for the different bases in the Chebyshev family of the second kind.

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 as long as they are all in the Chebyshev family, either first or second kind. The next letters in the matrix name uses the ‘short_name’ for all these different bases, see chebyshevu.bases.py.

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 BTTmat may be looked up as

>>> from shenfun.chebyshevu.matrices import mat
>>> from shenfun.chebyshevu.bases import Orthogonal as T
>>> B = mat[((T, 0), (T, 0))]

and an instance of the matrix can be created as

>>> B0 = T(10)
>>> BM = B((B0, 0), (B0, 0))
>>> import numpy as np
>>> d = {0: np.pi/2}
>>> [np.all(BM[k] == v) for k, v in d.items()]
[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]

To see that this is in fact the BUUmat:

>>> print(BM.__class__)
<class 'shenfun.chebyshevu.matrices.BUUmat'>

class shenfun.chebyshevu.matrices.AP1SDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

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\) chebyshevu.bases.Phi1, the trial function \(\psi_j \in\) chebyshev.bases.ShenDirichlet, and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

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

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

get_solver()

Return appropriate solver for self

Note

Fall back on generic Solve, which is using Scipy sparse matrices with splu/spsolve. This is still pretty fast.

class shenfun.chebyshevu.matrices.AP1SNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

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\) chebyshevu.bases.Phi1, the trial function \(\psi_j \in\) chebyshev.bases.ShenNeumann, and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

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

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

get_solver()

Return appropriate solver for self

Note

Fall back on generic Solve, which is using Scipy sparse matrices with splu/spsolve. This is still pretty fast.

class shenfun.chebyshevu.matrices.BP1SDmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

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

\[b_{kj}=(\psi_j, \phi_k)_w,\]

where the test function \(\phi_k \in\) chebyshevu.bases.Phi1, the trial \(\psi_j \in\) chebyshev.bases.ShenDirichlet, and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

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

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

class shenfun.chebyshevu.matrices.BP1SNmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

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

\[b_{kj}=(\psi_j, \phi_k)_w,\]

where the test function \(\phi_k \in\) chebyshevu.bases.Phi1, the trial function \(\psi_j \in\) chebyshev.bases.ShenNeumann, and test and trial spaces have dimensions of M and N, respectively.

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

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

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

class shenfun.chebyshevu.matrices.BUUmat(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)

Bases: SpectralMatrix

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

\[b_{kj}=(U_j, U_k)_w,\]

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

assemble(method)

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

• ‘exact’

• ‘quadrature’

Note

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

• legendre.matrix

• chebyshev.matrix

• chebyshevu.matrix

• ultraspherical.matrix

• fourier.matrix

• laguerre.matrix

• hermite.matrix

• jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

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

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

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

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

Module contents

Functionality for working with Chebyshev bases of second kind