rkstiff.derivatives

Spectral differentiation utilities for Fourier and Chebyshev grids

This module provides high-accuracy derivative operators using spectral methods suitable for stiff PDE solvers. Each function supports NumPy arrays and implements differentiation via spectral transforms or Chebyshev differentiation matrices.

Overview

Spectral differentiation is based on the principle that differentiation in physical space corresponds to multiplication by \((i k_x)^n\) in Fourier space, or repeated application of a differentiation matrix in Chebyshev space:

\[\frac{d^n u}{dx^n} = \mathcal{F}^{-1}\!\left[(i k_x)^n \mathcal{F}[u]\right], \quad \text{(Fourier methods)}\]

and

\[\frac{d^n u}{dx^n} = D^n u, \quad \text{(Chebyshev methods)},\]

where \(\mathcal{F}\) denotes the Fourier transform and \(D\) is the Chebyshev differentiation matrix.

Functions

• dx_rfft() — Derivatives using the real FFT (rFFT)

• dx_fft() — Derivatives using the complex FFT

• dx_cheb() — Derivatives using the Chebyshev differentiation matrix

Functions

dx_cheb(d_cheb_matrix, u[, n])	Compute the n-th derivative using a Chebyshev spectral differentiation matrix.
dx_fft(kx, u[, n])	Compute the n-th derivative of a complex-valued array using the FFT.
dx_rfft(kx, u[, n])	Compute the n-th derivative of a real-valued array using the real FFT (rFFT).

rkstiff.derivatives.dx_rfft(kx, u, n=1)

Compute the n-th derivative of a real-valued array using the real FFT (rFFT).

Parameters:

• kx (np.ndarray) – Wavenumbers of the spectral grid. For an input size \(N\), kx must have size \(N/2 + 1\). Typically generated via np.fft.rfftfreq(N, d=dx) * 2π.

• u (np.ndarray) – Real-valued input array of length \(N\). Must not be complex.

• n (int, optional) – Derivative order (non-negative integer). Default is 1.

Returns:

The n-th derivative of u (real-valued array of same shape).

Return type:

np.ndarray

Raises:

• TypeError – If n is not an integer or u is complex-valued.

• ValueError – If n is negative or kx has incompatible shape.

Notes

The derivative is computed using the spectral relation:

\[\frac{d^n u}{dx^n} = \mathcal{F}^{-1}\!\left[(i k_x)^n \mathcal{F}[u]\right],\]

where \(\mathcal{F}\) is the discrete Fourier transform (rFFT variant).

For a real-valued input array of size \(N\), the rFFT output has size \(N/2 + 1\). Therefore, kx must match that shape exactly.

Examples

>>> import numpy as np
>>> N = 128
>>> L = 2 * np.pi
>>> x = np.linspace(0, L, N, endpoint=False)
>>> kx = np.fft.rfftfreq(N, d=L/N) * 2 * np.pi
>>> u = np.sin(x)
>>> dudx = dx_rfft(kx, u)
>>> np.allclose(dudx, np.cos(x), atol=1e-10)
True

rkstiff.derivatives.dx_fft(kx, u, n=1)

Compute the n-th derivative of a complex-valued array using the FFT.

Parameters:

• kx (np.ndarray) – Wavenumbers of the spectral grid. Must match u.shape.

• u (np.ndarray) – Complex-valued input array.

• n (int, optional) – Derivative order (non-negative integer). Default is 1.

Returns:

The n-th derivative of u.

Return type:

np.ndarray

Notes

Fourier differentiation rule:

\[\frac{d^n u}{dx^n} = \mathcal{F}^{-1}\!\left[(i k_x)^n \mathcal{F}[u]\right].\]

This method is applicable to periodic complex fields or spectral-space data.

Examples

>>> import numpy as np
>>> N = 128
>>> L = 2 * np.pi
>>> x = np.linspace(0, L, N, endpoint=False)
>>> kx = np.fft.fftfreq(N, d=L/N) * 2 * np.pi
>>> u = np.exp(1j * x)
>>> dudx = dx_fft(kx, u)
>>> np.allclose(dudx, 1j * np.exp(1j * x), atol=1e-10)
True

rkstiff.derivatives.dx_cheb(d_cheb_matrix, u, n=1)

Compute the n-th derivative using a Chebyshev spectral differentiation matrix.

Parameters:

• d_cheb_matrix (np.ndarray) – Chebyshev differentiation matrix \(D \in \mathbb{R}^{N\times N}\), typically constructed with construct_x_dx_cheb.

• u (np.ndarray) – Function values sampled on the Chebyshev–Gauss–Lobatto grid. Shape (N,) or (N, M).

• n (int, optional) – Derivative order (non-negative integer). Default is 1.

Returns:

The n-th derivative of u, same shape as the input.

Return type:

np.ndarray

Notes

The Chebyshev differentiation operator is applied as:

\[\frac{du}{dx} = D u, \qquad \frac{d^n u}{dx^n} = D^n u.\]

For functions defined on a domain \([a,b]\), the differentiation matrix must be scaled by \(2/(b-a)\) to account for mapping from \([-1,1]\).

This method achieves spectral accuracy for smooth functions, but high derivative orders may amplify rounding errors.

Examples

>>> import numpy as np
>>> from rkstiff.grids import construct_x_dx_cheb
>>> x, D = construct_x_dx_cheb(n=20, a=-1, b=1)
>>> u = np.sin(x)
>>> dudx = dx_cheb(D, u)
>>> np.allclose(dudx, np.cos(x), atol=1e-10)
True
>>> d2udx2 = dx_cheb(D, u, n=2)
>>> np.allclose(d2udx2, -np.sin(x), atol=1e-8)
True