rkstiff.transforms

Transform utilities for axisymmetric and spectral-space computations

This module implements the Discrete Hankel Transform (DHT) for functions with cylindrical symmetry. The DHT provides a spectral decomposition for radial domains analogous to the Fourier transform for Cartesian coordinates, and is particularly useful in axially symmetric PDEs, such as paraxial wave equations, diffusion problems, or nonlinear Schrödinger-type systems.

Overview

The Hankel transform pair of order zero is defined as:

\[\begin{split}G(k) &= \int_0^{\infty} f(r) J_0(k r)\, r \, dr, \\ f(r) &= \int_0^{\infty} G(k) J_0(k r)\, k \, dk,\end{split}\]

where \(J_0\) is the Bessel function of the first kind of order zero.

The Discrete Hankel Transform (DHT) discretizes these integrals on a finite interval \(r \in [0, r_{\max}]\) using the zeros of \(J_0\). This quasi-discrete formulation is spectrally accurate for smooth, square-integrable functions on bounded domains.

Contents

HankelTransform — Implements forward and inverse DHT.

Classes

HankelTransform(nr[, rmax])

Discrete Hankel Transform (DHT) for axially symmetric functions.

class rkstiff.transforms.HankelTransform(nr, rmax=1.0)[source]

Bases: object

Discrete Hankel Transform (DHT) for axially symmetric functions.

This class implements the quasi-discrete Hankel transform (QDHT) following Guizar-Sicairos and Gutiérrez-Vega (2004), suitable for efficiently transforming between real-space and spectral-space representations of radially symmetric functions.

The continuous transform pair of order zero is:

\[\begin{split}G(k) &= \int_0^\infty f(r) J_0(k r)\, r \, dr, \\ f(r) &= \int_0^\infty G(k) J_0(k r)\, k \, dk,\end{split}\]

where \(J_0\) is the Bessel function of the first kind of order 0.

The discrete version used here approximates the above integrals by evaluating \(f(r)\) and \(G(k)\) on the grids:

\[r_i = \frac{j_i}{j_{N+1}}\,r_{\max}, \qquad k_i = \frac{j_i}{r_{\max}},\]

where \(j_i\) are the zeros of \(J_0\).

Parameters:

nr (int) – Number of radial sample points (≥ 4).
rmax (float, optional) – Maximum radial domain size. Default is 1.0.

r

Radial grid \([r_1, \ldots, r_N]\).

Type:: np.ndarray

kr

Spectral wavenumber grid corresponding to \(k_i\).

Type:: np.ndarray

_Y

Discrete Hankel transform matrix.

Type:: np.ndarray

_jN

Normalization constant (the \((N+1)\)-th zero of \(J_0\)).

Type:: float

References

Guizar-Sicairos, M. & Gutiérrez-Vega, J. C. Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields. J. Opt. Soc. Am. A 21, 53–58 (2004).

__init__(nr, rmax=1.0)[source]

Initialize a discrete Hankel transform on a radial domain.

hankel_matrix()[source]

Return the discrete Hankel transform matrix.

Returns:: Transform matrix \(Y_{ij} = 2 J_0(j_i j_j / j_{N+1}) / (j_{N+1} [J_1(j_j)]^2)\).
Return type:: np.ndarray, shape (nr, nr)

bessel_zeros()[source]

Return the zeros of the Bessel function \(J_0\).

Returns:: Array of the first \(N\) zeros of \(J_0\).
Return type:: np.ndarray

property nr: int: Number of radial grid points.

property rmax: float: Maximum radius of the domain.

ht(f)[source]

Compute the forward discrete Hankel transform.

Transforms a function \(f(r)\) from physical (radial) space to spectral (wavenumber) space according to:

\[G(k_i) \approx S \sum_{j=1}^{N} Y_{ij} f(r_j),\]

where \(S = r_{\max}^2 / j_{N+1}\) is the normalization factor.

Parameters:: f (np.ndarray) – Function values sampled at radial grid points r.
Returns:: Spectral coefficients sampled at corresponding wavenumbers kr.
Return type:: np.ndarray

Examples

>>> ht = HankelTransform(nr=64, rmax=5.0)
>>> f = np.exp(-ht.r**2)
>>> G = ht.ht(f)

iht(g)[source]

Compute the inverse discrete Hankel transform.

Transforms a spectral-space function \(G(k)\) back to real-space \(f(r)\) using:

\[f(r_i) \approx \frac{1}{S} \sum_{j=1}^{N} Y_{ij} G(k_j),\]

with the same transform matrix \(Y_{ij}\).

Parameters:: g (np.ndarray) – Spectral coefficients at wavenumber points kr.
Returns:: Real-space function evaluated at the radial grid r.
Return type:: np.ndarray

Examples

>>> ht = HankelTransform(nr=64, rmax=5.0)
>>> G = np.exp(-ht.kr**2)
>>> f = ht.iht(G)