shenfun.fourier package

Submodules

shenfun.fourier.bases module

Module for function spaces using Fourier exponentials.

A basis function \(\phi_k\) is given as

\[\phi_k(x) = \exp(ikx)\]

and an expansion is given as

(1)\[u(x) = \sum_{k=-N/2}^{N/2-1} \hat{u}_k \exp(ikx)\]

However, since \(\exp(ikx) = \exp(i(k \pm N)x)\) this expansion can also be written as an interpolator

(2)\[u(x) = \sum_{k=-N/2}^{N/2} \frac{\hat{u}_k}{c_k} \exp(ikx)\]

where \(c_{N/2} = c_{-N/2} = 2\), whereas \(c_k = 1\) for \(k=-N/2+1, ..., N/2-1\). Furthermore, \(\hat{u}_{N/2} = \hat{u}_{-N/2}\).

The interpolator form is used for computing odd derivatives. Otherwise, it makes no difference and therefore (1) is used in transforms, since this is the form expected by fftw. The inner product is defined as

\[(u, v) = \frac{1}{L} \int_{0}^{L} u \overline{v} dx\]

where \(\overline{v}\) is the complex conjugate of \(v\), and \(L\) is the length of the (periodic) domain.

class shenfun.fourier.bases.C2C(N, padding_factor=1, domain=(0, 2 * pi), dealias_direct=False, coordinates=None, **kw)[source]

Bases: FourierBase

Fourier function space for complex to complex transforms

A basis function \(\phi_k\) is given as

\[\phi_k(x) = \exp(ikx), \quad k =-N/2, -N/2+1, \ldots, N/2-1\]

An expansion is given as

\[u(x) = \sum_{k=-N/2}^{N/2-1} \hat{u}_k \exp(ikx).\]
Parameters:

  • N (int) – Number of quadrature points. Should be even for efficiency, but this is not required.

  • padding_factor (float, optional) – Factor for padding backward transforms. padding_factor=1.5 corresponds to a 3/2-rule for dealiasing.

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

  • dealias_direct (bool, optional) – True for dealiasing using 2/3-rule. Must be used with padding_factor = 1.

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

convolve(u, v, uv=None, fast=True)[source]

Convolution of u and v.

Parameters:

  • u (array)

  • v (array)

  • uv (array, optional)

  • fast (bool, optional) – Whether to use fast transforms in computing convolution

Note

Note that this method is only valid for 1D data, and that for multidimensional arrays one should use corresponding method in the TensorProductSpace class.

count_trailing_zeros(u, reltol=1e-12, abstol=1e-15)[source]
get_mask_nyquist(bcast=True)[source]
shape(forward_output=True)[source]

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()[source]
slice()[source]

Return index set of current space

wavenumbers(bcast=True, scaled=False, eliminate_highest_freq=False)[source]

Return the wavenumbermesh

Parameters:

bcast (bool) – Whether or not to broadcast to dimensions() if an instance of this basis belongs to a TensorProductSpace

class shenfun.fourier.bases.R2C(N, padding_factor=1.0, domain=(0, 2 * pi), dealias_direct=False, coordinates=None, **kw)[source]

Bases: FourierBase

Fourier function space for real to complex transforms

A basis function \(\phi_k\) is given as

\[\phi_k(x) = \exp(ikx), \quad k =-N/2, -N/2+1, \ldots, N/2-1\]

An expansion is given as

\[u(x) = \sum_{k=-N/2}^{N/2-1} \hat{u}_k \exp(ikx).\]

Howewer, due to Hermitian symmetry \(\hat{u}_k = \overline{\hat{u}_{-k}}\), which is taken advantage of in this function space. Specifically, Fourier transforms make use of real-to-complex and complex-to-real algorithms, see FFTW and rfftn of mpi4py-fft.

Parameters:

  • N (int) – Number of quadrature points. Should be even for efficiency, but this is not required.

  • padding_factor (float, optional) – Factor for padding backward transforms. padding_factor=1.5 corresponds to a 3/2-rule for dealiasing.

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

  • dealias_direct (bool, optional) – True for dealiasing using 2/3-rule. Must be used with padding_factor = 1.

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

convolve(u, v, uv=None, fast=True)[source]

Convolution of u and v.

Parameters:

  • u (array)

  • v (array)

  • uv (array, optional)

  • fast (bool, optional) – Whether to use fast transforms in computing convolution

Note

Note that this method is only valid for 1D data, and that for multidimensional arrays one should use corresponding method in the TensorProductSpace class.

get_mask_nyquist(bcast=True)[source]

Return None or an array with zeros for Nyquist coefficients and one otherwise

Parameters:

bcast (boolean, optional) – If True then broadcast returned mask array to dimensions of the TensorProductSpace this base belongs to.

shape(forward_output=True)[source]

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()[source]
slice()[source]

Return index set of current space

wavenumbers(bcast=True, scaled=False, eliminate_highest_freq=False)[source]

Return the wavenumbermesh

Parameters:

bcast (bool) – Whether or not to broadcast to dimensions() if an instance of this basis belongs to a TensorProductSpace

shenfun.fourier.matrices module

Module for handling Fourier diagonal matrices

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

Bases: SpectralMatrix

Stiffness matrix for \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where

\[a_{kj}=(\exp(i k x), \cos^2(x) \exp(i l x)'')\]

where 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.

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

Bases: SpectralMatrix

Stiffness matrix for \(A=(a_{kj}) \in \mathbb{R}^{M \times N}\), where

\[a_{kj}=(\exp(i k x), \cos(x) \exp(i l x)'')\]

where 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.

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

Bases: SpectralMatrix

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

\[b_{kj}=(\exp(i k x), \cos^2(x) \exp(i l x))\]

where 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.

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

Bases: SpectralMatrix

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

\[b_{kj}=(\exp(i k x), \cos(x) \exp(i l x))\]

where 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.

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

Bases: SpectralMatrix

Derivative matrix for \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where

\[c_{kj}=(\exp(i k x), \sin(2x)/2 \exp(i l x)')\]

where 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.

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

Bases: SpectralMatrix

Derivative matrix for \(C=(c_{kj}) \in \mathbb{R}^{M \times N}\), where

\[c_{kj}=(\exp(i k x), \sin(x) \exp(i l x)')\]

where 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.

class shenfun.fourier.matrices.FourierMatDict[source]

Bases: SpectralMatDict

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

Bases: SpectralMatrix

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.

solve(b, u=None, axis=0, constraints=())[source]

Solve matrix system Au = b

where A is the current matrix (self)

Parameters:

  • b (array) – Array of right hand side on entry and solution on exit unless u is provided.

  • u (array, optional) – Output array

  • axis (int, optional) – The axis over which to solve for if b and u are multi- dimensional

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

Note

Vectors may be one- or multidimensional.

Module contents

shenfun.fourier.energy_fourier(u, T)[source]

Compute the energy of u using Parceval’s theorem

\[\int abs(u)^2 dx = N*\sum abs(u_hat)^2\]
Parameters:

  • u (Array) – The Fourier coefficients

  • T (TensorProductSpace)

See https://en.wikipedia.org/wiki/Parseval’s_theorem