Spectral Poisson / Helmholtz Solvers

Direct spectral solvers for 2-D elliptic PDEs on rectangular domains, using DST-I (Dirichlet), DCT-II (Neumann), or FFT (periodic) spectral methods.

DST-I (Dirichlet BCs)

`finitevolx.solve_poisson_dst(rhs, dx, dy, *, approximation='fd2')`

Solve ∇²ψ = f with homogeneous Dirichlet BCs using DST-I.

Convenience wrapper around :func:solve_helmholtz_dst with lambda_=0.

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side at interior grid points.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Solution ψ, same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_poisson_dst(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve ∇²ψ = f with homogeneous Dirichlet BCs using DST-I.

    Convenience wrapper around :func:`solve_helmholtz_dst` with ``lambda_=0``.

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side at interior grid points.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Solution ψ, same shape as *rhs*.
    """
    return solve_helmholtz_dst(rhs, dx, dy, lambda_=0.0, approximation=approximation)

`finitevolx.solve_helmholtz_dst(rhs, dx, dy, lambda_=0.0, *, approximation='fd2')`

Solve (∇² − λ)ψ = f with homogeneous Dirichlet BCs using DST-I.

The input rhs lives on the interior grid (boundary values are implicitly zero: ψ = 0 on all four edges).

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side at interior grid points.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`lambda_`	`float`	Helmholtz parameter λ. Default: 0.0 (Poisson).	`0.0`
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Solution ψ at interior grid points, same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_helmholtz_dst(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    lambda_: float = 0.0,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve (∇² − λ)ψ = f with homogeneous Dirichlet BCs using DST-I.

    The input *rhs* lives on the **interior** grid (boundary values are
    implicitly zero: ψ = 0 on all four edges).

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side at interior grid points.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    lambda_ : float
        Helmholtz parameter λ.  Default: 0.0 (Poisson).
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Solution ψ at interior grid points, same shape as *rhs*.
    """
    Ny, Nx = rhs.shape
    rhs_hat = dstn(rhs, type=1, axes=[0, 1])
    eigx = _eig_1d(
        dst1_eigenvalues, dst1_eigenvalues_ps, Nx, dx, (Nx + 1) * dx, approximation
    )
    eigy = _eig_1d(
        dst1_eigenvalues, dst1_eigenvalues_ps, Ny, dy, (Ny + 1) * dy, approximation
    )
    eig2d = eigy[:, None] + eigx[None, :] - lambda_
    psi_hat = rhs_hat / eig2d
    return idstn(psi_hat, type=1, axes=[0, 1])

DCT-II (Neumann BCs)

`finitevolx.solve_poisson_dct(rhs, dx, dy, *, approximation='fd2')`

Solve ∇²ψ = f with homogeneous Neumann BCs using DCT-II.

The Poisson problem has a one-dimensional null space (constant solutions). This function fixes the gauge by forcing the domain-mean of ψ to zero.

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Zero-mean solution ψ, same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_poisson_dct(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve ∇²ψ = f with homogeneous Neumann BCs using DCT-II.

    The Poisson problem has a one-dimensional null space (constant solutions).
    This function fixes the gauge by forcing the domain-mean of ψ to zero.

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Zero-mean solution ψ, same shape as *rhs*.
    """
    return solve_helmholtz_dct(rhs, dx, dy, lambda_=0.0, approximation=approximation)

`finitevolx.solve_helmholtz_dct(rhs, dx, dy, lambda_=0.0, *, approximation='fd2')`

Solve (∇² − λ)ψ = f with homogeneous Neumann BCs using DCT-II.

When lambda_ == 0 the (0,0) mode is singular (Λ[0,0] = 0, corresponding to the constant null mode of the Neumann Laplacian). This is handled by setting ψ̂[0,0] = 0 (zero-mean gauge).

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`lambda_`	`float`	Helmholtz parameter λ. Default: 0.0 (Poisson).	`0.0`
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Solution ψ, same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_helmholtz_dct(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    lambda_: float = 0.0,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve (∇² − λ)ψ = f with homogeneous Neumann BCs using DCT-II.

    When ``lambda_ == 0`` the (0,0) mode is singular (Λ[0,0] = 0,
    corresponding to the constant null mode of the Neumann Laplacian).
    This is handled by setting ψ̂[0,0] = 0 (zero-mean gauge).

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    lambda_ : float
        Helmholtz parameter λ.  Default: 0.0 (Poisson).
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Solution ψ, same shape as *rhs*.
    """
    Ny, Nx = rhs.shape
    rhs_hat = dctn(rhs, type=2, axes=[0, 1])
    eigx = _eig_1d(
        dct2_eigenvalues, dct2_eigenvalues_ps, Nx, dx, Nx * dx, approximation
    )
    eigy = _eig_1d(
        dct2_eigenvalues, dct2_eigenvalues_ps, Ny, dy, Ny * dy, approximation
    )
    eig2d = eigy[:, None] + eigx[None, :] - lambda_
    is_null = eig2d[0, 0] == 0.0
    eig2d_safe = eig2d.at[0, 0].set(jnp.where(is_null, 1.0, eig2d[0, 0]))
    psi_hat = rhs_hat / eig2d_safe
    psi_hat = psi_hat.at[0, 0].set(jnp.where(is_null, 0.0, psi_hat[0, 0]))
    return idctn(psi_hat, type=2, axes=[0, 1])

FFT (Periodic BCs)

`finitevolx.solve_poisson_fft(rhs, dx, dy, *, approximation='fd2')`

Solve ∇²ψ = f with periodic BCs using the 2-D FFT.

Convenience wrapper around :func:solve_helmholtz_fft with lambda_=0.

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side on the periodic domain.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Zero-mean solution ψ (real), same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_poisson_fft(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve ∇²ψ = f with periodic BCs using the 2-D FFT.

    Convenience wrapper around :func:`solve_helmholtz_fft` with ``lambda_=0``.

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side on the periodic domain.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Zero-mean solution ψ (real), same shape as *rhs*.
    """
    return solve_helmholtz_fft(rhs, dx, dy, lambda_=0.0, approximation=approximation)

`finitevolx.solve_helmholtz_fft(rhs, dx, dy, lambda_=0.0, *, approximation='fd2')`

Solve (∇² − λ)ψ = f with periodic BCs using the 2-D FFT.

Spectral algorithm:

Forward 2-D FFT: f̂ = FFT2(f) [Ny, Nx]
Eigenvalue matrix: Λ[j,i] = Λ_j^y + Λ_i^x − λ [Ny, Nx] where Λ^x, Λ^y are FD2 or PS eigenvalues (see approximation).
Spectral division: ψ̂[j,i] = f̂[j,i] / Λ[j,i] [Ny, Nx]
Inverse 2-D FFT: ψ = Re(IFFT2(ψ̂)) [Ny, Nx]

When lambda_ == 0 the (0,0) mode is singular (Λ[0,0] = 0). This is handled by setting ψ̂[0,0] = 0 (zero-mean gauge).

Parameters:

Name	Type	Description	Default
`rhs`	`Float[Array, 'Ny Nx']`	Right-hand side on the periodic domain.	required
`dx`	`float`	Grid spacing in x.	required
`dy`	`float`	Grid spacing in y.	required
`lambda_`	`float`	Helmholtz parameter λ. Default: 0.0 (Poisson).	`0.0`
`approximation`	`('fd2', 'spectral')`	Eigenvalue type. Default: `"fd2"`.	`"fd2"`

Returns:

Type	Description
`Float[Array, 'Ny Nx']`	Solution ψ (real), same shape as rhs.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/solvers.py

def solve_helmholtz_fft(
    rhs: Float[Array, "Ny Nx"],
    dx: float,
    dy: float,
    lambda_: float = 0.0,
    *,
    approximation: Approximation = "fd2",
) -> Float[Array, "Ny Nx"]:
    """Solve (∇² − λ)ψ = f with periodic BCs using the 2-D FFT.

    Spectral algorithm:

    1. Forward 2-D FFT:  f̂ = FFT2(f)                        [Ny, Nx]
    2. Eigenvalue matrix:
           Λ[j,i] = Λ_j^y + Λ_i^x − λ                      [Ny, Nx]
       where Λ^x, Λ^y are FD2 or PS eigenvalues (see *approximation*).
    3. Spectral division:  ψ̂[j,i] = f̂[j,i] / Λ[j,i]       [Ny, Nx]
    4. Inverse 2-D FFT:  ψ = Re(IFFT2(ψ̂))                  [Ny, Nx]

    When ``lambda_ == 0`` the (0,0) mode is singular (Λ[0,0] = 0).
    This is handled by setting ψ̂[0,0] = 0 (zero-mean gauge).

    Parameters
    ----------
    rhs : Float[Array, "Ny Nx"]
        Right-hand side on the periodic domain.
    dx : float
        Grid spacing in x.
    dy : float
        Grid spacing in y.
    lambda_ : float
        Helmholtz parameter λ.  Default: 0.0 (Poisson).
    approximation : {"fd2", "spectral"}
        Eigenvalue type.  Default: ``"fd2"``.

    Returns
    -------
    Float[Array, "Ny Nx"]
        Solution ψ (real), same shape as *rhs*.
    """
    Ny, Nx = rhs.shape
    rhs_hat = jnp.fft.fft2(rhs)
    eigx = _eig_1d(fft_eigenvalues, fft_eigenvalues_ps, Nx, dx, Nx * dx, approximation)
    eigy = _eig_1d(fft_eigenvalues, fft_eigenvalues_ps, Ny, dy, Ny * dy, approximation)
    eig2d = eigy[:, None] + eigx[None, :] - lambda_
    is_null = eig2d[0, 0] == 0.0
    eig2d_safe = eig2d.at[0, 0].set(jnp.where(is_null, 1.0, eig2d[0, 0]))
    psi_hat = rhs_hat / eig2d_safe
    psi_hat = psi_hat.at[0, 0].set(
        jnp.where(is_null, jnp.zeros_like(psi_hat[0, 0]), psi_hat[0, 0])
    )
    return jnp.real(jnp.fft.ifft2(psi_hat))

Eigenvalue Helpers

`finitevolx.dst1_eigenvalues(N, dx)`

1-D Laplacian eigenvalues for the DST-I basis (homogeneous Dirichlet BCs).

Returns the eigenvalues of the second-order finite-difference Laplacian

(L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

on N interior points with v_0 = v_{N+1} = 0. The DST-I diagonalises this operator, yielding:

λ_k = −4/dx² · sin²(π(k+1) / (2(N+1))),   k = 0, …, N−1

All eigenvalues are strictly negative (the discrete Dirichlet Laplacian is negative definite).

Parameters:

Name	Type	Description	Default
`N`	`int`	Number of interior grid points (excluding the two zero-boundary points).	required
`dx`	`float`	Uniform grid spacing.	required

Returns:

Type	Description
`Float[Array, ' N']`	1-D array of N eigenvalues, ordered k = 0, …, N−1. All < 0.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/eigenvalues.py

def dst1_eigenvalues(N: int, dx: float) -> Float[Array, " N"]:
    """1-D Laplacian eigenvalues for the DST-I basis (homogeneous Dirichlet BCs).

    Returns the eigenvalues of the second-order finite-difference Laplacian

        (L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

    on *N* interior points with v_0 = v_{N+1} = 0.  The DST-I diagonalises
    this operator, yielding:

        λ_k = −4/dx² · sin²(π(k+1) / (2(N+1))),   k = 0, …, N−1

    All eigenvalues are strictly negative (the discrete Dirichlet Laplacian
    is negative definite).

    Parameters
    ----------
    N : int
        Number of interior grid points (excluding the two zero-boundary points).
    dx : float
        Uniform grid spacing.

    Returns
    -------
    Float[Array, " N"]
        1-D array of *N* eigenvalues, ordered k = 0, …, N−1.  All < 0.
    """
    k = jnp.arange(N)
    return -4.0 / dx**2 * jnp.sin(jnp.pi * (k + 1) / (2 * (N + 1))) ** 2

`finitevolx.dct2_eigenvalues(N, dx)`

1-D Laplacian eigenvalues for the DCT-II basis (homogeneous Neumann BCs).

Returns the eigenvalues of the second-order finite-difference Laplacian

(L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

on N points with ∂v/∂n = 0 at both boundaries (ghost-point Neumann). The DCT-II diagonalises this operator, yielding:

λ_k = −4/dx² · sin²(πk / (2N)),   k = 0, …, N−1

The k=0 eigenvalue is exactly zero (the constant mode), corresponding to the one-dimensional null space of the Neumann Laplacian.

Parameters:

Name	Type	Description	Default
`N`	`int`	Number of grid points.	required
`dx`	`float`	Uniform grid spacing.	required

Returns:

Type	Description
`Float[Array, ' N']`	1-D array of N eigenvalues, ordered k = 0, …, N−1. λ_0 = 0; all others < 0.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/eigenvalues.py

def dct2_eigenvalues(N: int, dx: float) -> Float[Array, " N"]:
    """1-D Laplacian eigenvalues for the DCT-II basis (homogeneous Neumann BCs).

    Returns the eigenvalues of the second-order finite-difference Laplacian

        (L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

    on *N* points with ∂v/∂n = 0 at both boundaries (ghost-point Neumann).
    The DCT-II diagonalises this operator, yielding:

        λ_k = −4/dx² · sin²(πk / (2N)),   k = 0, …, N−1

    The k=0 eigenvalue is exactly zero (the constant mode), corresponding
    to the one-dimensional null space of the Neumann Laplacian.

    Parameters
    ----------
    N : int
        Number of grid points.
    dx : float
        Uniform grid spacing.

    Returns
    -------
    Float[Array, " N"]
        1-D array of *N* eigenvalues, ordered k = 0, …, N−1.
        λ_0 = 0; all others < 0.
    """
    k = jnp.arange(N)
    return -4.0 / dx**2 * jnp.sin(jnp.pi * k / (2 * N)) ** 2

`finitevolx.fft_eigenvalues(N, dx)`

1-D Laplacian eigenvalues for the FFT basis (periodic BCs).

Returns the eigenvalues of the second-order finite-difference Laplacian

(L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

on N points with periodic boundary conditions v_{N} = v_0. The DFT diagonalises this operator, yielding:

λ_k = −4/dx² · sin²(πk / N),   k = 0, …, N−1

The k=0 eigenvalue is exactly zero (the constant / zero-wavenumber mode), corresponding to the one-dimensional null space of the periodic Laplacian.

Parameters:

Name	Type	Description	Default
`N`	`int`	Number of grid points in one period.	required
`dx`	`float`	Uniform grid spacing.	required

Returns:

Type	Description
`Float[Array, ' N']`	1-D array of N eigenvalues in FFT order (k = 0, …, N−1). λ_0 = 0; all others ≤ 0.

Source code in .venv/lib/python3.12/site-packages/spectraldiffx/_src/fourier/eigenvalues.py

def fft_eigenvalues(N: int, dx: float) -> Float[Array, " N"]:
    """1-D Laplacian eigenvalues for the FFT basis (periodic BCs).

    Returns the eigenvalues of the second-order finite-difference Laplacian

        (L·v)_i = (v_{i-1} − 2v_i + v_{i+1}) / dx²

    on *N* points with periodic boundary conditions v_{N} = v_0.  The DFT
    diagonalises this operator, yielding:

        λ_k = −4/dx² · sin²(πk / N),   k = 0, …, N−1

    The k=0 eigenvalue is exactly zero (the constant / zero-wavenumber mode),
    corresponding to the one-dimensional null space of the periodic Laplacian.

    Parameters
    ----------
    N : int
        Number of grid points in one period.
    dx : float
        Uniform grid spacing.

    Returns
    -------
    Float[Array, " N"]
        1-D array of *N* eigenvalues in FFT order (k = 0, …, N−1).
        λ_0 = 0; all others ≤ 0.
    """
    k = jnp.arange(N)
    return -4.0 / dx**2 * jnp.sin(jnp.pi * k / N) ** 2