Chebyshev Transforms, Calculus & Quadrature

Transforms

ChebyshevTransform1D

Bases: Module

Forward/inverse 1D Chebyshev transform.

Thin wrapper around :meth:ChebyshevGrid1D.transform that exposes :meth:to_spectral / :meth:from_spectral methods. Useful when you want to pass the transform around as a first-class object (e.g. into a PDE residual, or to a higher-order library).

Attributes:

Name	Type	Description
grid	ChebyshevGrid1D	Underlying grid carrying the node convention (GL or Gauss).

Examples:

>>> import jax.numpy as jnp
>>> grid = ChebyshevGrid1D.from_N_L(N=16, L=1.0)
>>> cheb = ChebyshevTransform1D(grid=grid)
>>> u = jnp.sin(jnp.pi * grid.x)
>>> a = cheb.to_spectral(u)
>>> u_roundtrip = cheb.from_spectral(a)  # ≈ u

Source code in spectraldiffx/_src/chebyshev/transforms.py

class ChebyshevTransform1D(eqx.Module):
    """Forward/inverse 1D Chebyshev transform.

    Thin wrapper around :meth:`ChebyshevGrid1D.transform` that exposes
    :meth:`to_spectral` / :meth:`from_spectral` methods.  Useful when you
    want to pass the transform around as a first-class object (e.g. into
    a PDE residual, or to a higher-order library).

    Attributes
    ----------
    grid : ChebyshevGrid1D
        Underlying grid carrying the node convention (GL or Gauss).

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid1D.from_N_L(N=16, L=1.0)
    >>> cheb = ChebyshevTransform1D(grid=grid)
    >>> u = jnp.sin(jnp.pi * grid.x)
    >>> a = cheb.to_spectral(u)
    >>> u_roundtrip = cheb.from_spectral(a)  # ≈ u
    """

    grid: ChebyshevGrid1D

    def to_spectral(self, u: Num[Array, "Npts"]) -> Num[Array, "Npts"]:
        """Forward Chebyshev transform: nodal values → coefficients aₖ."""
        return self.grid.transform(u, inverse=False)

    def from_spectral(self, a: Num[Array, "Npts"]) -> Num[Array, "Npts"]:
        """Inverse Chebyshev transform: coefficients aₖ → nodal values."""
        return self.grid.transform(a, inverse=True)

    def __call__(
        self, u: Num[Array, "Npts"], inverse: bool = False
    ) -> Num[Array, "Npts"]:
        """Alias for ``self.grid.transform(u, inverse=inverse)``."""
        return self.grid.transform(u, inverse=inverse)

Functions

to_spectral(u)

Forward Chebyshev transform: nodal values → coefficients aₖ.

Source code in spectraldiffx/_src/chebyshev/transforms.py

def to_spectral(self, u: Num[Array, "Npts"]) -> Num[Array, "Npts"]:
    """Forward Chebyshev transform: nodal values → coefficients aₖ."""
    return self.grid.transform(u, inverse=False)

from_spectral(a)

Inverse Chebyshev transform: coefficients aₖ → nodal values.

Source code in spectraldiffx/_src/chebyshev/transforms.py

def from_spectral(self, a: Num[Array, "Npts"]) -> Num[Array, "Npts"]:
    """Inverse Chebyshev transform: coefficients aₖ → nodal values."""
    return self.grid.transform(a, inverse=True)

__call__(u, inverse=False)

Alias for self.grid.transform(u, inverse=inverse).

Source code in spectraldiffx/_src/chebyshev/transforms.py

def __call__(
    self, u: Num[Array, "Npts"], inverse: bool = False
) -> Num[Array, "Npts"]:
    """Alias for ``self.grid.transform(u, inverse=inverse)``."""
    return self.grid.transform(u, inverse=inverse)

ChebyshevTransform2D

Bases: Module

Forward/inverse 2D Chebyshev transform (tensor product of 1D).

Attributes:

Name	Type	Description
grid	ChebyshevGrid2D	Underlying 2D grid.

Examples:

>>> import jax.numpy as jnp
>>> grid = ChebyshevGrid2D.from_N_L(Nx=16, Ny=16, Lx=1.0, Ly=1.0)
>>> cheb = ChebyshevTransform2D(grid=grid)
>>> X, Y = grid.X
>>> u = jnp.sin(jnp.pi * X) * jnp.cos(jnp.pi * Y)
>>> a = cheb.to_spectral(u)
>>> u_roundtrip = cheb.from_spectral(a)

Source code in spectraldiffx/_src/chebyshev/transforms.py

class ChebyshevTransform2D(eqx.Module):
    """Forward/inverse 2D Chebyshev transform (tensor product of 1D).

    Attributes
    ----------
    grid : ChebyshevGrid2D
        Underlying 2D grid.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid2D.from_N_L(Nx=16, Ny=16, Lx=1.0, Ly=1.0)
    >>> cheb = ChebyshevTransform2D(grid=grid)
    >>> X, Y = grid.X
    >>> u = jnp.sin(jnp.pi * X) * jnp.cos(jnp.pi * Y)
    >>> a = cheb.to_spectral(u)
    >>> u_roundtrip = cheb.from_spectral(a)
    """

    grid: ChebyshevGrid2D

    def to_spectral(self, u: Num[Array, "Nypts Nxpts"]) -> Num[Array, "Nypts Nxpts"]:
        """Forward 2D Chebyshev transform: nodal values → coefficients a_{jk}."""
        return self.grid.transform(u, inverse=False)

    def from_spectral(self, a: Num[Array, "Nypts Nxpts"]) -> Num[Array, "Nypts Nxpts"]:
        """Inverse 2D Chebyshev transform."""
        return self.grid.transform(a, inverse=True)

    def __call__(
        self,
        u: Num[Array, "Nypts Nxpts"],
        inverse: bool = False,
    ) -> Num[Array, "Nypts Nxpts"]:
        """Alias for ``self.grid.transform(u, inverse=inverse)``."""
        return self.grid.transform(u, inverse=inverse)

Functions

to_spectral(u)

Forward 2D Chebyshev transform: nodal values → coefficients a_{jk}.

Source code in spectraldiffx/_src/chebyshev/transforms.py

def to_spectral(self, u: Num[Array, "Nypts Nxpts"]) -> Num[Array, "Nypts Nxpts"]:
    """Forward 2D Chebyshev transform: nodal values → coefficients a_{jk}."""
    return self.grid.transform(u, inverse=False)

from_spectral(a)

Inverse 2D Chebyshev transform.

Source code in spectraldiffx/_src/chebyshev/transforms.py

def from_spectral(self, a: Num[Array, "Nypts Nxpts"]) -> Num[Array, "Nypts Nxpts"]:
    """Inverse 2D Chebyshev transform."""
    return self.grid.transform(a, inverse=True)

__call__(u, inverse=False)

Alias for self.grid.transform(u, inverse=inverse).

Source code in spectraldiffx/_src/chebyshev/transforms.py

def __call__(
    self,
    u: Num[Array, "Nypts Nxpts"],
    inverse: bool = False,
) -> Num[Array, "Nypts Nxpts"]:
    """Alias for ``self.grid.transform(u, inverse=inverse)``."""
    return self.grid.transform(u, inverse=inverse)

cheb_dealias_product(grid, a, b)

Compute the dealiased pointwise product a·b on a Chebyshev grid.

Implements a 2/3-style truncation in Chebyshev-coefficient space:

1. Forward-transform a and b to coefficient space.
2. Zero modes with index > 2N/3 on both inputs.
3. Inverse-transform, multiply pointwise.
4. Forward-transform the product, zero modes with index > 2N/3,
   and inverse-transform once more.

This is the Chebyshev analogue of Orszag's 2/3 rule for Fourier grids. It is exact for products whose highest relevant Chebyshev mode is below the cut-off, and otherwise prevents aliasing into the retained modes at the cost of some high-mode truncation error.

Parameters:

Name	Type	Description	Default
grid	ChebyshevGrid1D or ChebyshevGrid2D	Grid providing the transform; must have dealias='2/3' to have any effect (otherwise this function falls back to the naive product a * b).	required
a	Num[Array, ...]	Nodal fields on the same grid. Shapes (Npts,) for 1D or (Nypts, Nxpts) for 2D.	required
b	Num[Array, ...]	Nodal fields on the same grid. Shapes (Npts,) for 1D or (Nypts, Nxpts) for 2D.	required

Returns:

Type	Description
Num[Array, ...]	Dealiased pointwise product on the same grid.

Examples:

>>> import jax.numpy as jnp
>>> grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0, dealias="2/3")
>>> x = grid.x
>>> u = jnp.sin(jnp.pi * x)
>>> v = jnp.cos(jnp.pi * x)
>>> uv = dealias_product(grid, u, v)  # ≈ ½ sin(2πx)

Source code in spectraldiffx/_src/chebyshev/transforms.py

def dealias_product(
    grid: ChebyshevGrid1D | ChebyshevGrid2D,
    a: Num[Array, "..."],
    b: Num[Array, "..."],
) -> Num[Array, "..."]:
    """Compute the dealiased pointwise product a·b on a Chebyshev grid.

    Implements a 2/3-style truncation in Chebyshev-coefficient space:

        1. Forward-transform a and b to coefficient space.
        2. Zero modes with index > 2N/3 on both inputs.
        3. Inverse-transform, multiply pointwise.
        4. Forward-transform the product, zero modes with index > 2N/3,
           and inverse-transform once more.

    This is the Chebyshev analogue of Orszag's 2/3 rule for Fourier grids.
    It is exact for products whose highest relevant Chebyshev mode is
    below the cut-off, and otherwise prevents aliasing into the retained
    modes at the cost of some high-mode truncation error.

    Parameters
    ----------
    grid : ChebyshevGrid1D or ChebyshevGrid2D
        Grid providing the transform; must have ``dealias='2/3'`` to have
        any effect (otherwise this function falls back to the naive
        product ``a * b``).
    a, b : Num[Array, ...]
        Nodal fields on the same grid.  Shapes ``(Npts,)`` for 1D or
        ``(Nypts, Nxpts)`` for 2D.

    Returns
    -------
    Num[Array, ...]
        Dealiased pointwise product on the same grid.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0, dealias="2/3")
    >>> x = grid.x
    >>> u = jnp.sin(jnp.pi * x)
    >>> v = jnp.cos(jnp.pi * x)
    >>> uv = dealias_product(grid, u, v)  # ≈ ½ sin(2πx)
    """
    if grid.dealias != "2/3":
        # No dealiasing requested → return the naive product.
        return a * b

    if isinstance(grid, ChebyshevGrid1D):
        n_modes = grid.N + 1 if grid.node_type == "gauss-lobatto" else grid.N
        cutoff = int(2 * grid.N / 3)

        a_hat = grid.transform(a)
        b_hat = grid.transform(b)
        a_hat = _apply_dealias_1d(a_hat, n_modes, cutoff)
        b_hat = _apply_dealias_1d(b_hat, n_modes, cutoff)
        a_f = grid.transform(a_hat, inverse=True)
        b_f = grid.transform(b_hat, inverse=True)

        ab = a_f * b_f
        ab_hat = _apply_dealias_1d(grid.transform(ab), n_modes, cutoff)
        return grid.transform(ab_hat, inverse=True)

    if isinstance(grid, ChebyshevGrid2D):
        nx_modes = grid.Nx + 1 if grid.node_type == "gauss-lobatto" else grid.Nx
        ny_modes = grid.Ny + 1 if grid.node_type == "gauss-lobatto" else grid.Ny
        cutoff_x = int(2 * grid.Nx / 3)
        cutoff_y = int(2 * grid.Ny / 3)
        mask_x = (jnp.arange(nx_modes) <= cutoff_x).astype(a.dtype)
        mask_y = (jnp.arange(ny_modes) <= cutoff_y).astype(a.dtype)
        mask = mask_y[:, None] * mask_x[None, :]

        a_hat = grid.transform(a) * mask
        b_hat = grid.transform(b) * mask
        a_f = grid.transform(a_hat, inverse=True)
        b_f = grid.transform(b_hat, inverse=True)

        ab_hat = grid.transform(a_f * b_f) * mask
        return grid.transform(ab_hat, inverse=True)

    raise TypeError(
        f"grid must be a ChebyshevGrid1D or ChebyshevGrid2D, got {type(grid).__name__}."
    )

Coefficient-space calculus

chebyshev_derivative_coeffs(a, L=1.0, order=1)

Differentiate a Chebyshev series in coefficient space (last axis).

For u(x) = Σₖ aₖ Tₖ(x/L), the derivative u'(x) = Σₖ a'ₖ Tₖ(x/L) has

cₖ a'ₖ = a'ₖ₊₂ + 2(k+1) aₖ₊₁,     a'_N = a'_{N+1} = 0

(c₀ = 2, cₖ = 1 otherwise), scaled by 1/L for the map x = L·ξ. Unrolling the recurrence gives a closed form: a'ₖ is a sum over the modes j > k of opposite parity,

a'ₖ = (2 / (cₖ L)) Σ_{j>k, j+k odd} j·aⱼ

which we evaluate with two reverse cumulative sums (one per parity), so the whole derivative costs O(N) and parallelises on accelerators — no sequential scan is needed.

Parameters:

Name	Type	Description	Default
a	Num[Array, '... Nmodes']	Chebyshev coefficients along the last axis (e.g. from :meth:ChebyshevGrid1D.transform).	required
L	float	Domain half-length. Default 1.0.	1.0
order	int	Derivative order (≥ 0).	1

Returns:

Type	Description
Num[Array, '... Nmodes']	Coefficients of the order-th derivative (same shape as a; the top order modes are zero).

Examples:

d/dx T₃(x) = 3 U₂(x) = 3 (T₀ + 2 T₂):

>>> import jax.numpy as jnp
>>> a = jnp.array([0.0, 0.0, 0.0, 1.0, 0.0])
>>> chebyshev_derivative_coeffs(a)  # ≈ [3, 0, 6, 0, 0]

Source code in spectraldiffx/_src/chebyshev/transforms.py

def chebyshev_derivative_coeffs(
    a: Num[Array, "... Nmodes"],
    L: float = 1.0,
    order: int = 1,
) -> Num[Array, "... Nmodes"]:
    """Differentiate a Chebyshev series in coefficient space (last axis).

    For u(x) = Σₖ aₖ Tₖ(x/L), the derivative u'(x) = Σₖ a'ₖ Tₖ(x/L) has

        cₖ a'ₖ = a'ₖ₊₂ + 2(k+1) aₖ₊₁,     a'_N = a'_{N+1} = 0

    (c₀ = 2, cₖ = 1 otherwise), scaled by 1/L for the map x = L·ξ.
    Unrolling the recurrence gives a closed form: a'ₖ is a sum over the
    modes j > k of opposite parity,

        a'ₖ = (2 / (cₖ L)) Σ_{j>k, j+k odd} j·aⱼ

    which we evaluate with two reverse cumulative sums (one per parity),
    so the whole derivative costs O(N) and parallelises on accelerators —
    no sequential ``scan`` is needed.

    Parameters
    ----------
    a : Num[Array, "... Nmodes"]
        Chebyshev coefficients along the last axis (e.g. from
        :meth:`ChebyshevGrid1D.transform`).
    L : float
        Domain half-length.  Default 1.0.
    order : int
        Derivative order (≥ 0).

    Returns
    -------
    Num[Array, "... Nmodes"]
        Coefficients of the ``order``-th derivative (same shape as ``a``;
        the top ``order`` modes are zero).

    Examples
    --------
    d/dx T₃(x) = 3 U₂(x) = 3 (T₀ + 2 T₂):

    >>> import jax.numpy as jnp
    >>> a = jnp.array([0.0, 0.0, 0.0, 1.0, 0.0])
    >>> chebyshev_derivative_coeffs(a)  # ≈ [3, 0, 6, 0, 0]
    """
    if order < 0:
        raise ValueError(f"order must be >= 0, got {order}")
    n_modes = a.shape[-1]
    k = jnp.arange(n_modes)
    odd = (k % 2).astype(bool)
    c = jnp.where(k == 0, 2.0, 1.0)

    def rev_cumsum(v: Array) -> Array:
        return jnp.flip(jnp.cumsum(jnp.flip(v, axis=-1), axis=-1), axis=-1)

    out = a
    for _ in range(order):
        b = 2.0 * k * out  # 2j·aⱼ
        s_odd = rev_cumsum(jnp.where(odd, b, 0.0))  # Σ_{j≥k, j odd} 2j aⱼ
        s_even = rev_cumsum(jnp.where(odd, 0.0, b))  # Σ_{j≥k, j even} 2j aⱼ
        # Mode k collects the opposite-parity tail (j = k is excluded
        # automatically because it has the same parity as k).
        out = jnp.where(odd, s_even, s_odd) / (c * L)
    return out

chebyshev_antiderivative_coeffs(a, L=1.0)

Indefinite integral of a Chebyshev series, vanishing at x = −L.

Uses ∫T₀ = T₁, ∫T₁ = T₂/4 and, for k ≥ 2,

∫Tₖ dξ = T_{k+1} / (2(k+1)) − T_{k−1} / (2(k−1))

so that the antiderivative U(x) = Σₖ Bₖ Tₖ(x/L) has

Bₖ = L (cₖ₋₁ aₖ₋₁ − aₖ₊₁) / (2k),     k ≥ 1   (c₀ = 2, a_{N+1} = 0)

and B₀ is fixed by U(−L) = Σₖ Bₖ (−1)ᵏ = 0.

The exact antiderivative has degree N+1; its top mode B_{N+1} = L a_N / (2(N+1)) is dropped to keep the array shape, which is exact whenever a_N = 0 (e.g. for any dealiased or resolved field).

Parameters:

Name	Type	Description	Default
a	Num[Array, '... Nmodes']	Chebyshev coefficients along the last axis.	required
L	float	Domain half-length.	1.0

Returns:

Type	Description
Num[Array, '... Nmodes']	Coefficients of U(x) = ∫_{−L}^{x} u(s) ds.

Examples:

∫_{−1}^{x} 1 ds = x + 1 = T₀ + T₁:

>>> import jax.numpy as jnp
>>> chebyshev_antiderivative_coeffs(jnp.array([1.0, 0.0, 0.0]))  # ≈ [1, 1, 0]

Source code in spectraldiffx/_src/chebyshev/transforms.py

def chebyshev_antiderivative_coeffs(
    a: Num[Array, "... Nmodes"],
    L: float = 1.0,
) -> Num[Array, "... Nmodes"]:
    """Indefinite integral of a Chebyshev series, vanishing at x = −L.

    Uses ∫T₀ = T₁, ∫T₁ = T₂/4 and, for k ≥ 2,

        ∫Tₖ dξ = T_{k+1} / (2(k+1)) − T_{k−1} / (2(k−1))

    so that the antiderivative U(x) = Σₖ Bₖ Tₖ(x/L) has

        Bₖ = L (cₖ₋₁ aₖ₋₁ − aₖ₊₁) / (2k),     k ≥ 1   (c₀ = 2, a_{N+1} = 0)

    and B₀ is fixed by U(−L) = Σₖ Bₖ (−1)ᵏ = 0.

    The exact antiderivative has degree N+1; its top mode
    B_{N+1} = L a_N / (2(N+1)) is dropped to keep the array shape, which
    is exact whenever a_N = 0 (e.g. for any dealiased or resolved field).

    Parameters
    ----------
    a : Num[Array, "... Nmodes"]
        Chebyshev coefficients along the last axis.
    L : float
        Domain half-length.

    Returns
    -------
    Num[Array, "... Nmodes"]
        Coefficients of U(x) = ∫_{−L}^{x} u(s) ds.

    Examples
    --------
    ∫_{−1}^{x} 1 ds = x + 1 = T₀ + T₁:

    >>> import jax.numpy as jnp
    >>> chebyshev_antiderivative_coeffs(jnp.array([1.0, 0.0, 0.0]))  # ≈ [1, 1, 0]
    """
    n_modes = a.shape[-1]
    k = jnp.arange(n_modes)
    c = jnp.where(k == 0, 2.0, 1.0)
    zero = jnp.zeros_like(a[..., :1])
    a_prev = jnp.concatenate([zero, (c * a)[..., :-1]], axis=-1)  # cₖ₋₁ aₖ₋₁
    a_next = jnp.concatenate([a[..., 1:], zero], axis=-1)  # aₖ₊₁
    k_safe = jnp.where(k == 0, 1, k)
    B = jnp.where(k == 0, 0.0, L * (a_prev - a_next) / (2.0 * k_safe))
    sign = jnp.where(k % 2 == 0, 1.0, -1.0)
    B0 = -jnp.sum(B * sign, axis=-1)
    return B.at[..., 0].set(B0)

chebyshev_integral_coeffs(a, L=1.0)

Definite integral ∫_{−L}^{L} u(x) dx from Chebyshev coefficients.

Since ∫_{−1}^{1} Tₖ(ξ) dξ = 2 / (1 − k²) for even k and 0 for odd k,

∫_{−L}^{L} u dx = L Σ_{k even} 2 aₖ / (1 − k²)

On Gauss–Lobatto nodes this is identical to Clenshaw–Curtis quadrature; on Gauss nodes it is the corresponding Fejér-type rule.

Parameters:

Name	Type	Description	Default
a	Num[Array, '... Nmodes']	Chebyshev coefficients along the last axis.	required
L	float	Domain half-length.	1.0

Returns:

Type	Description
Num[Array, '...']	The integral, reducing the last axis.

Source code in spectraldiffx/_src/chebyshev/transforms.py

def chebyshev_integral_coeffs(
    a: Num[Array, "... Nmodes"],
    L: float = 1.0,
) -> Num[Array, "..."]:
    """Definite integral ∫_{−L}^{L} u(x) dx from Chebyshev coefficients.

    Since ∫_{−1}^{1} Tₖ(ξ) dξ = 2 / (1 − k²) for even k and 0 for odd k,

        ∫_{−L}^{L} u dx = L Σ_{k even} 2 aₖ / (1 − k²)

    On Gauss–Lobatto nodes this is identical to Clenshaw–Curtis
    quadrature; on Gauss nodes it is the corresponding Fejér-type rule.

    Parameters
    ----------
    a : Num[Array, "... Nmodes"]
        Chebyshev coefficients along the last axis.
    L : float
        Domain half-length.

    Returns
    -------
    Num[Array, "..."]
        The integral, reducing the last axis.
    """
    n_modes = a.shape[-1]
    k = jnp.arange(n_modes)
    even = k % 2 == 0
    k2m1 = jnp.where(even, 1.0 - k * k, 1.0)
    w = jnp.where(even, 2.0 / k2m1, 0.0)
    return L * jnp.sum(a * w, axis=-1)

Clenshaw–Curtis quadrature

clenshaw_curtis_weights(N, L=1.0)

Clenshaw–Curtis quadrature weights on Gauss–Lobatto nodes of [−L, L].

Parameters:

Name	Type	Description	Default
N	int	Chebyshev polynomial degree. The grid has N+1 Gauss–Lobatto nodes.	required
L	float	Domain half-length (default 1). The weights scale linearly with L.	1.0

Returns:

Type	Description
Float[Array, 'Npts']	Weights w such that ∫_{−L}^{L} f(x) dx ≈ Σⱼ w[j] · f(xⱼ). Exact for polynomials of degree ≤ N (and degree N+1 when N is even). For analytic integrands CC converges at the same asymptotic rate as Gauss quadrature despite the smaller degree of precision.

Examples:

Integrate exp(x) on [−1, 1] (exact value e − 1/e):

>>> import jax.numpy as jnp
>>> w = clenshaw_curtis_weights(N=32, L=1.0)
>>> x = jnp.cos(jnp.pi * jnp.arange(33) / 32)
>>> float(jnp.sum(w * jnp.exp(x)))

Source code in spectraldiffx/_src/chebyshev/quadrature.py

def clenshaw_curtis_weights(N: int, L: float = 1.0) -> Float[Array, "Npts"]:
    """Clenshaw–Curtis quadrature weights on Gauss–Lobatto nodes of [−L, L].

    Parameters
    ----------
    N : int
        Chebyshev polynomial degree.  The grid has N+1 Gauss–Lobatto nodes.
    L : float
        Domain half-length (default 1).  The weights scale linearly with L.

    Returns
    -------
    Float[Array, "Npts"]
        Weights ``w`` such that ``∫_{−L}^{L} f(x) dx ≈ Σⱼ w[j] · f(xⱼ)``.
        Exact for polynomials of degree ≤ N (and degree N+1 when N is
        even).  For analytic integrands CC converges at the same
        asymptotic rate as Gauss quadrature despite the smaller degree of
        precision.

    Examples
    --------
    Integrate exp(x) on [−1, 1] (exact value e − 1/e):

    >>> import jax.numpy as jnp
    >>> w = clenshaw_curtis_weights(N=32, L=1.0)
    >>> x = jnp.cos(jnp.pi * jnp.arange(33) / 32)
    >>> float(jnp.sum(w * jnp.exp(x)))  # doctest: +SKIP
    """
    return jnp.asarray(_cc_weights_numpy(N, L))

clenshaw_curtis_integrate_1d(grid, f)

Integrate a 1D nodal field over [−L, L] using Clenshaw–Curtis.

The grid must use Gauss–Lobatto nodes (Gauss nodes would require a different quadrature rule — Gauss–Chebyshev — which is not provided here because it is less accurate for smooth non-periodic f).

Parameters:

Name	Type	Description	Default
grid	ChebyshevGrid1D	Must have node_type == 'gauss-lobatto'.	required
f	Num[Array, 'Npts']	Nodal values of the integrand.	required

Returns:

Type	Description
Float[Array, '']	Scalar approximation to ∫_{−L}^{L} f(x) dx.

Examples:

>>> import jax.numpy as jnp
>>> grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0)
>>> f = jnp.exp(grid.x)
>>> I = clenshaw_curtis_integrate_1d(grid, f)  # ≈ e − 1/e

Source code in spectraldiffx/_src/chebyshev/quadrature.py

def clenshaw_curtis_integrate_1d(
    grid: ChebyshevGrid1D,
    f: Num[Array, "Npts"],
) -> Float[Array, ""]:
    """Integrate a 1D nodal field over [−L, L] using Clenshaw–Curtis.

    The grid must use Gauss–Lobatto nodes (Gauss nodes would require a
    different quadrature rule — Gauss–Chebyshev — which is not provided
    here because it is less accurate for smooth non-periodic f).

    Parameters
    ----------
    grid : ChebyshevGrid1D
        Must have ``node_type == 'gauss-lobatto'``.
    f : Num[Array, "Npts"]
        Nodal values of the integrand.

    Returns
    -------
    Float[Array, ""]
        Scalar approximation to ∫_{−L}^{L} f(x) dx.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0)
    >>> f = jnp.exp(grid.x)
    >>> I = clenshaw_curtis_integrate_1d(grid, f)  # ≈ e − 1/e
    """
    if grid.node_type != "gauss-lobatto":
        raise ValueError(
            "clenshaw_curtis_integrate_1d requires Gauss–Lobatto nodes; got "
            f"node_type={grid.node_type!r}."
        )
    w = clenshaw_curtis_weights(grid.N, grid.L)
    return jnp.sum(w * f)

clenshaw_curtis_integrate_2d(grid, f)

Integrate a 2D nodal field over [−Lx, Lx] × [−Ly, Ly].

Uses the tensor product of 1D Clenshaw–Curtis rules:

∫∫ f(x, y) dx dy ≈ Σⱼᵢ w_y[j] · w_x[i] · f[j, i]

Parameters:

Name	Type	Description	Default
grid	ChebyshevGrid2D	Must use Gauss–Lobatto nodes in both directions.	required
f	Num[Array, 'Nypts Nxpts']	Nodal values of the integrand on the (Nᵧ+1, Nₓ+1) grid.	required

Returns:

Type	Description
Float[Array, '']	Scalar approximation to ∫∫ f dA.

Examples:

>>> import jax.numpy as jnp
>>> grid = ChebyshevGrid2D.from_N_L(Nx=24, Ny=24, Lx=1.0, Ly=1.0)
>>> X, Y = grid.X
>>> f = jnp.exp(X + Y)
>>> I = clenshaw_curtis_integrate_2d(grid, f)  # ≈ (e − 1/e)²

Source code in spectraldiffx/_src/chebyshev/quadrature.py

def clenshaw_curtis_integrate_2d(
    grid: ChebyshevGrid2D,
    f: Num[Array, "Nypts Nxpts"],
) -> Float[Array, ""]:
    """Integrate a 2D nodal field over [−Lx, Lx] × [−Ly, Ly].

    Uses the tensor product of 1D Clenshaw–Curtis rules:

        ∫∫ f(x, y) dx dy ≈ Σⱼᵢ w_y[j] · w_x[i] · f[j, i]

    Parameters
    ----------
    grid : ChebyshevGrid2D
        Must use Gauss–Lobatto nodes in both directions.
    f : Num[Array, "Nypts Nxpts"]
        Nodal values of the integrand on the (Nᵧ+1, Nₓ+1) grid.

    Returns
    -------
    Float[Array, ""]
        Scalar approximation to ∫∫ f dA.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid2D.from_N_L(Nx=24, Ny=24, Lx=1.0, Ly=1.0)
    >>> X, Y = grid.X
    >>> f = jnp.exp(X + Y)
    >>> I = clenshaw_curtis_integrate_2d(grid, f)  # ≈ (e − 1/e)²
    """
    if grid.node_type != "gauss-lobatto":
        raise ValueError(
            "clenshaw_curtis_integrate_2d requires Gauss–Lobatto nodes; got "
            f"node_type={grid.node_type!r}."
        )
    wx = clenshaw_curtis_weights(grid.Nx, grid.Lx)  # (Nxpts,)
    wy = clenshaw_curtis_weights(grid.Ny, grid.Ly)  # (Nypts,)
    W = wy[:, None] * wx[None, :]  # (Nypts, Nxpts)
    return jnp.sum(W * f)