Skip to content

Chebyshev Derivative Operators

ChebyshevDerivative1D

Bases: Module

1D Chebyshev derivative operator using the precomputed differentiation matrix.

Mathematical Formulation:

For a function u(x) sampled at Chebyshev nodes xⱼ on [-L, L]:

du/dx |ⱼ = Σₖ D_{jk} u_k

where D is the (N+1)×(N+1) [or N×N for Gauss] differentiation matrix precomputed in ChebyshevGrid1D.

Higher-order derivatives are obtained by matrix powers

d²u/dx² = D @ (D @ u) = D² @ u

grid : ChebyshevGrid1D
    The 1D Chebyshev grid containing the differentiation matrix.
Source code in spectraldiffx/_src/chebyshev/operators.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
class ChebyshevDerivative1D(eqx.Module):
    """
    1D Chebyshev derivative operator using the precomputed differentiation matrix.

    Mathematical Formulation:
    -------------------------
    For a function u(x) sampled at Chebyshev nodes xⱼ on [-L, L]:

        du/dx |ⱼ = Σₖ D_{jk} u_k

    where D is the (N+1)×(N+1) [or N×N for Gauss] differentiation matrix
    precomputed in ChebyshevGrid1D.

    Higher-order derivatives are obtained by matrix powers:
        d²u/dx² = D @ (D @ u) = D² @ u

    Attributes:
    -----------
        grid : ChebyshevGrid1D
            The 1D Chebyshev grid containing the differentiation matrix.
    """

    grid: ChebyshevGrid1D

    def __call__(self, u: Array, order: int = 1) -> Float[Array, "N1"]:
        """
        Compute the n-th derivative of a field using the differentiation matrix.

        Operation:
            du/dx    = D @ u
            d²u/dx²  = D @ D @ u

        Parameters:
        -----------
        u : Array [N+1] for GL, [N] for Gauss
            Physical-space field values at Chebyshev nodes.
        order : int
            Derivative order (1 = first, 2 = second, ...). Default is 1.

        Returns:
        --------
        dnu_dxn : Array [N+1] or [N]
            n-th derivative at Chebyshev nodes.
        """
        D = self.grid.D
        result = u
        for _ in range(order):
            result = D @ result
        return result

    def gradient(self, u: Array) -> Float[Array, "N1"]:
        """
        Compute the first derivative du/dx.

        Parameters:
        -----------
        u : Array [N1]
            Physical-space field at Chebyshev nodes.

        Returns:
        --------
        du_dx : Array [N1]
        """
        return self(u, order=1)

    def laplacian(self, u: Array) -> Float[Array, "N1"]:
        """
        Compute the second derivative d²u/dx².

        Parameters:
        -----------
        u : Array [N1]
            Physical-space field at Chebyshev nodes.

        Returns:
        --------
        d2u_dx2 : Array [N1]
        """
        return self(u, order=2)

Functions

__call__(u, order=1)

Compute the n-th derivative of a field using the differentiation matrix.

Operation

du/dx = D @ u d²u/dx² = D @ D @ u

Parameters:

u : Array [N+1] for GL, [N] for Gauss Physical-space field values at Chebyshev nodes. order : int Derivative order (1 = first, 2 = second, ...). Default is 1.

Returns:

dnu_dxn : Array [N+1] or [N] n-th derivative at Chebyshev nodes.

Source code in spectraldiffx/_src/chebyshev/operators.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def __call__(self, u: Array, order: int = 1) -> Float[Array, "N1"]:
    """
    Compute the n-th derivative of a field using the differentiation matrix.

    Operation:
        du/dx    = D @ u
        d²u/dx²  = D @ D @ u

    Parameters:
    -----------
    u : Array [N+1] for GL, [N] for Gauss
        Physical-space field values at Chebyshev nodes.
    order : int
        Derivative order (1 = first, 2 = second, ...). Default is 1.

    Returns:
    --------
    dnu_dxn : Array [N+1] or [N]
        n-th derivative at Chebyshev nodes.
    """
    D = self.grid.D
    result = u
    for _ in range(order):
        result = D @ result
    return result

gradient(u)

Compute the first derivative du/dx.

Parameters:

u : Array [N1] Physical-space field at Chebyshev nodes.

Returns:

du_dx : Array [N1]

Source code in spectraldiffx/_src/chebyshev/operators.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def gradient(self, u: Array) -> Float[Array, "N1"]:
    """
    Compute the first derivative du/dx.

    Parameters:
    -----------
    u : Array [N1]
        Physical-space field at Chebyshev nodes.

    Returns:
    --------
    du_dx : Array [N1]
    """
    return self(u, order=1)

laplacian(u)

Compute the second derivative d²u/dx².

Parameters:

u : Array [N1] Physical-space field at Chebyshev nodes.

Returns:

d2u_dx2 : Array [N1]

Source code in spectraldiffx/_src/chebyshev/operators.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def laplacian(self, u: Array) -> Float[Array, "N1"]:
    """
    Compute the second derivative d²u/dx².

    Parameters:
    -----------
    u : Array [N1]
        Physical-space field at Chebyshev nodes.

    Returns:
    --------
    d2u_dx2 : Array [N1]
    """
    return self(u, order=2)

ChebyshevDerivative2D

Bases: Module

2D Chebyshev derivative operators on [-Lx, Lx] × [-Ly, Ly].

Mathematical Formulation:

For u(x, y) on a (Ny_pts, Nx_pts) grid:

∂u/∂x [j,i] = (u @ Dxᵀ)[j,i]   (differentiation along axis 1)
∂u/∂y [j,i] = (Dy @ u)[j,i]     (differentiation along axis 0)

Laplacian: ∇²u = ∂²u/∂x² + ∂²u/∂y²

grid : ChebyshevGrid2D
    The 2D Chebyshev grid containing Dx and Dy matrices.
Source code in spectraldiffx/_src/chebyshev/operators.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
class ChebyshevDerivative2D(eqx.Module):
    """
    2D Chebyshev derivative operators on [-Lx, Lx] × [-Ly, Ly].

    Mathematical Formulation:
    -------------------------
    For u(x, y) on a (Ny_pts, Nx_pts) grid:

        ∂u/∂x [j,i] = (u @ Dxᵀ)[j,i]   (differentiation along axis 1)
        ∂u/∂y [j,i] = (Dy @ u)[j,i]     (differentiation along axis 0)

    Laplacian: ∇²u = ∂²u/∂x² + ∂²u/∂y²

    Attributes:
    -----------
        grid : ChebyshevGrid2D
            The 2D Chebyshev grid containing Dx and Dy matrices.
    """

    grid: ChebyshevGrid2D

    def gradient(
        self, u: Array
    ) -> tuple[Float[Array, "Ny1 Nx1"], Float[Array, "Ny1 Nx1"]]:
        """
        Compute partial derivatives (∂u/∂x, ∂u/∂y).

        Parameters:
        -----------
        u : Array [Ny_pts, Nx_pts]
            2D field at Chebyshev nodes.

        Returns:
        --------
        (du_dx, du_dy) : tuple of Arrays [Ny_pts, Nx_pts]
        """
        Dx = self.grid.Dx
        Dy = self.grid.Dy
        du_dx = u @ Dx.T  # apply Dx along x-axis (axis 1)
        du_dy = Dy @ u  # apply Dy along y-axis (axis 0)
        return du_dx, du_dy

    def laplacian(self, u: Array) -> Float[Array, "Ny1 Nx1"]:
        """
        Compute the 2D Laplacian ∇²u = ∂²u/∂x² + ∂²u/∂y².

        Uses precomputed second-derivative matrices Dx2 = Dx@Dx and Dy2 = Dy@Dy
        (stored on the grid) to avoid recomputing the O(N³) matrix products
        on every call.

        Parameters:
        -----------
        u : Array [Ny_pts, Nx_pts]

        Returns:
        --------
        lap_u : Array [Ny_pts, Nx_pts]
        """
        d2u_dx2 = u @ self.grid.Dx2.T
        d2u_dy2 = self.grid.Dy2 @ u
        return d2u_dx2 + d2u_dy2

    def divergence(self, vx: Array, vy: Array) -> Float[Array, "Ny1 Nx1"]:
        """
        Compute the divergence ∇·V = ∂vx/∂x + ∂vy/∂y.

        Parameters:
        -----------
        vx : Array [Ny_pts, Nx_pts]
            x-component of the vector field.
        vy : Array [Ny_pts, Nx_pts]
            y-component of the vector field.

        Returns:
        --------
        div : Array [Ny_pts, Nx_pts]
        """
        Dx = self.grid.Dx
        Dy = self.grid.Dy
        dvx_dx = vx @ Dx.T
        dvy_dy = Dy @ vy
        return dvx_dx + dvy_dy

    def curl(self, vx: Array, vy: Array) -> Float[Array, "Ny1 Nx1"]:
        """
        Compute the 2D scalar curl ζ = ∂vy/∂x - ∂vx/∂y.

        Parameters:
        -----------
        vx : Array [Ny_pts, Nx_pts]
        vy : Array [Ny_pts, Nx_pts]

        Returns:
        --------
        curl : Array [Ny_pts, Nx_pts]
        """
        Dx = self.grid.Dx
        Dy = self.grid.Dy
        dvy_dx = vy @ Dx.T
        dvx_dy = Dy @ vx
        return dvy_dx - dvx_dy

    def advection_scalar(
        self, vx: Array, vy: Array, q: Array
    ) -> Float[Array, "Ny1 Nx1"]:
        """
        Compute scalar advection (V·∇)q = vx·∂q/∂x + vy·∂q/∂y.

        Parameters:
        -----------
        vx : Array [Ny_pts, Nx_pts]
            x-velocity.
        vy : Array [Ny_pts, Nx_pts]
            y-velocity.
        q : Array [Ny_pts, Nx_pts]
            Scalar tracer field.

        Returns:
        --------
        adv : Array [Ny_pts, Nx_pts]
        """
        dq_dx, dq_dy = self.gradient(q)
        return vx * dq_dx + vy * dq_dy

Functions

gradient(u)

Compute partial derivatives (∂u/∂x, ∂u/∂y).

Parameters:

u : Array [Ny_pts, Nx_pts] 2D field at Chebyshev nodes.

Returns:

(du_dx, du_dy) : tuple of Arrays [Ny_pts, Nx_pts]

Source code in spectraldiffx/_src/chebyshev/operators.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def gradient(
    self, u: Array
) -> tuple[Float[Array, "Ny1 Nx1"], Float[Array, "Ny1 Nx1"]]:
    """
    Compute partial derivatives (∂u/∂x, ∂u/∂y).

    Parameters:
    -----------
    u : Array [Ny_pts, Nx_pts]
        2D field at Chebyshev nodes.

    Returns:
    --------
    (du_dx, du_dy) : tuple of Arrays [Ny_pts, Nx_pts]
    """
    Dx = self.grid.Dx
    Dy = self.grid.Dy
    du_dx = u @ Dx.T  # apply Dx along x-axis (axis 1)
    du_dy = Dy @ u  # apply Dy along y-axis (axis 0)
    return du_dx, du_dy

laplacian(u)

Compute the 2D Laplacian ∇²u = ∂²u/∂x² + ∂²u/∂y².

Uses precomputed second-derivative matrices Dx2 = Dx@Dx and Dy2 = Dy@Dy (stored on the grid) to avoid recomputing the O(N³) matrix products on every call.

Parameters:

u : Array [Ny_pts, Nx_pts]

Returns:

lap_u : Array [Ny_pts, Nx_pts]

Source code in spectraldiffx/_src/chebyshev/operators.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def laplacian(self, u: Array) -> Float[Array, "Ny1 Nx1"]:
    """
    Compute the 2D Laplacian ∇²u = ∂²u/∂x² + ∂²u/∂y².

    Uses precomputed second-derivative matrices Dx2 = Dx@Dx and Dy2 = Dy@Dy
    (stored on the grid) to avoid recomputing the O(N³) matrix products
    on every call.

    Parameters:
    -----------
    u : Array [Ny_pts, Nx_pts]

    Returns:
    --------
    lap_u : Array [Ny_pts, Nx_pts]
    """
    d2u_dx2 = u @ self.grid.Dx2.T
    d2u_dy2 = self.grid.Dy2 @ u
    return d2u_dx2 + d2u_dy2

divergence(vx, vy)

Compute the divergence ∇·V = ∂vx/∂x + ∂vy/∂y.

Parameters:

vx : Array [Ny_pts, Nx_pts] x-component of the vector field. vy : Array [Ny_pts, Nx_pts] y-component of the vector field.

Returns:

div : Array [Ny_pts, Nx_pts]

Source code in spectraldiffx/_src/chebyshev/operators.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
def divergence(self, vx: Array, vy: Array) -> Float[Array, "Ny1 Nx1"]:
    """
    Compute the divergence ∇·V = ∂vx/∂x + ∂vy/∂y.

    Parameters:
    -----------
    vx : Array [Ny_pts, Nx_pts]
        x-component of the vector field.
    vy : Array [Ny_pts, Nx_pts]
        y-component of the vector field.

    Returns:
    --------
    div : Array [Ny_pts, Nx_pts]
    """
    Dx = self.grid.Dx
    Dy = self.grid.Dy
    dvx_dx = vx @ Dx.T
    dvy_dy = Dy @ vy
    return dvx_dx + dvy_dy

curl(vx, vy)

Compute the 2D scalar curl ζ = ∂vy/∂x - ∂vx/∂y.

Parameters:

vx : Array [Ny_pts, Nx_pts] vy : Array [Ny_pts, Nx_pts]

Returns:

curl : Array [Ny_pts, Nx_pts]

Source code in spectraldiffx/_src/chebyshev/operators.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def curl(self, vx: Array, vy: Array) -> Float[Array, "Ny1 Nx1"]:
    """
    Compute the 2D scalar curl ζ = ∂vy/∂x - ∂vx/∂y.

    Parameters:
    -----------
    vx : Array [Ny_pts, Nx_pts]
    vy : Array [Ny_pts, Nx_pts]

    Returns:
    --------
    curl : Array [Ny_pts, Nx_pts]
    """
    Dx = self.grid.Dx
    Dy = self.grid.Dy
    dvy_dx = vy @ Dx.T
    dvx_dy = Dy @ vx
    return dvy_dx - dvx_dy

advection_scalar(vx, vy, q)

Compute scalar advection (V·∇)q = vx·∂q/∂x + vy·∂q/∂y.

Parameters:

vx : Array [Ny_pts, Nx_pts] x-velocity. vy : Array [Ny_pts, Nx_pts] y-velocity. q : Array [Ny_pts, Nx_pts] Scalar tracer field.

Returns:

adv : Array [Ny_pts, Nx_pts]

Source code in spectraldiffx/_src/chebyshev/operators.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def advection_scalar(
    self, vx: Array, vy: Array, q: Array
) -> Float[Array, "Ny1 Nx1"]:
    """
    Compute scalar advection (V·∇)q = vx·∂q/∂x + vy·∂q/∂y.

    Parameters:
    -----------
    vx : Array [Ny_pts, Nx_pts]
        x-velocity.
    vy : Array [Ny_pts, Nx_pts]
        y-velocity.
    q : Array [Ny_pts, Nx_pts]
        Scalar tracer field.

    Returns:
    --------
    adv : Array [Ny_pts, Nx_pts]
    """
    dq_dx, dq_dy = self.gradient(q)
    return vx * dq_dx + vy * dq_dy