psiop

psiop.py — Symbolic‑numerical toolkit for pseudo‑differential operators in 1D/2D

Overview

The psiop module provides a unified framework for manipulating pseudo‑differential operators (ΨDOs) in one and two spatial dimensions. It combines symbolic construction of operator symbols (using SymPy) with numerical evaluation and visualisation (using NumPy/SciPy/Matplotlib). The package is designed for researchers and students working in microlocal analysis, spectral theory, and the numerical analysis of PDEs.

Key features include:

  • Symbol creation either from an explicit expression (symbol mode) or by automatic extraction from a differential operator (auto mode).

  • Computation of asymptotic expansions of symbols at high frequencies.

  • Determination of the operator order (homogeneity degree) via symbolic/numerical heuristics.

  • Asymptotic composition of operators in Kohn‑Nirenberg or Weyl quantisation.

  • Construction of formal left/right inverses and the formal adjoint.

  • Symbol of the exponential exp(tP) via series expansion.

  • Semiclassical trace formula (symbolic or numerical).

  • Hamiltonian flow associated with the principal symbol (symplectic phase‑space dynamics).

  • Pseudospectrum computation and visualisation with adaptive grid refinement and parallelisation.

  • Rich visualisation suite: symbol amplitude/phase, cotangent‑fiber structure, characteristic set, group velocity field, micro‑support, Hamiltonian trajectories.

  • Interactive dashboard (ipywidgets) for real‑time exploration of the symbol.

Mathematical background

A pseudo‑differential operator P acting on functions of x ∈ ℝⁿ is formally defined by its symbol p(x,ξ), a function on phase space T*ℝⁿ. The action on a function u is given by the Kohn‑Nirenberg quantisation

(P u)(x) = (2π)^{-n} ∫_{ℝⁿ} e^{i x·ξ} p(x,ξ) û(ξ) dξ ,

where û is the Fourier transform of u. If the symbol does not depend on x, the operator reduces to a Fourier multiplier. For a general spatially varying symbol, the above representation provides a rigorous extension of differential operators.

The asymptotic behaviour of the symbol as |ξ| → ∞ determines many properties of the operator. The principal symbol pₘ is the homogeneous component of highest order m. When the symbol is non‑homogeneous, the module attempts to estimate the effective order by expanding in inverse powers of |ξ| (or in a radial variable for 2D).

Symbolic calculus allows one to compose operators asymptotically. For two symbols p and q, the symbol of the composition P ∘ Q is given by an asymptotic series

(p ∘ q)(x,ξ) ~ Σ_{α} (i)^{-|α|}/α! ∂_ξ^α p(x,ξ) ∂_x^α q(x,ξ)

in the Kohn‑Nirenberg convention (a similar expansion exists for the Weyl star product). Truncating this series yields approximate compositions valid for high frequencies or slowly varying symbols.

Formal inverses and adjoints are constructed via similar asymptotic series, assuming the principal symbol never vanishes.

The Hamiltonian flow generated by the principal symbol describes the propagation of singularities along bicharacteristics – a cornerstone of microlocal analysis.

The pseudospectrum σ_ε(P) is the set of λ ∈ ℂ for which ‖(P-λI)^{-1}‖ ≥ ε^{-1}. It captures the near‑spectral behaviour of non‑normal operators and is visualised via contour plots of the resolvent norm.

References

class psiop.PseudoDifferentialOperator(expr, vars_x, var_u=None, mode='symbol', quantization='weyl')[source]

Bases: object

Pseudo-differential operator with dynamic symbol evaluation on spatial grids. Supports both 1D and 2D operators, and can be defined explicitly (symbol mode) or extracted automatically from symbolic equations (auto mode).

Parameters:

  • expr (sympy expression) – Symbolic expression representing the pseudo-differential symbol.

  • vars_x (list of sympy symbols) – Spatial variables (e.g., [x] for 1D, [x, y] for 2D).

  • var_u (sympy function, optional) – Function u(x, t) used in auto mode to extract the operator symbol.

  • mode (str, {'symbol', 'auto'}) –

    • ‘symbol’: directly uses expr as the operator symbol.

    • ’auto’: computes the symbol automatically by applying expr to exp(i x ξ).

dim

Spatial dimension (1 or 2).

Type:

int

fft, ifft

Fast Fourier transform and inverse (scipy.fft or scipy.fft2).

Type:

callable

p_func

Evaluated symbol function ready for numerical use.

Type:

callable

Notes

  • In ‘symbol’ mode, expr should be expressed in terms of spatial variables and frequency variables (ξ, η).

  • In ‘auto’ mode, the symbol is derived by applying the differential expression to a complex exponential.

  • Frequency variables are internally named ‘xi’ and ‘eta’ for consistency.

  • Uses numpy for numerical evaluation and scipy.fft for FFT operations.

Examples

>>> # Example 1: 1D Laplacian operator (symbol mode)
>>> from sympy import symbols
>>> x, xi = symbols('x xi', real=True)
>>> op = PseudoDifferentialOperator(expr=xi**2, vars_x=[x], mode='symbol')

>>> # Example 2: 1D transport operator (auto mode)
>>> from sympy import Function
>>> u = Function('u')
>>> expr = u(x).diff(x)
>>> op = PseudoDifferentialOperator(expr=expr, vars_x=[x], var_u=u(x), mode='auto')
animate_singularity(xi0=5.0, eta0=0.0, x0=0.0, y0=0.0, tmax=4.0, n_frames=100, projection=None)[source]

Animate the propagation of a singularity under the Hamiltonian flow.

This method visualizes how a singularity (x₀, y₀, ξ₀, η₀) evolves in phase space according to the Hamiltonian dynamics induced by the principal symbol of the operator. The animation integrates the Hamiltonian equations of motion and supports various projections: position (x-y), frequency (ξ-η), or mixed phase space coordinates.

Parameters:

  • xi0 (float) – Initial frequency components (ξ₀, η₀).

  • eta0 (float) – Initial frequency components (ξ₀, η₀).

  • x0 (float) – Initial spatial coordinates (x₀, y₀).

  • y0 (float) – Initial spatial coordinates (x₀, y₀).

  • tmax (float) – Total time of integration (final animation time).

  • n_frames (int) – Number of frames in the resulting animation.

  • projection (str or None) –

    Type of projection to display:

    • ’position’ : x vs y (or x alone in 1D)

    • ’frequency’: ξ vs η (or ξ alone in 1D)

    • ’phase’ : mixed coordinates like x vs ξ or x vs η

    If None, defaults to ‘phase’ in 1D and ‘position’ in 2D.

Returns:

Animation object that can be displayed interactively in Jupyter notebooks or saved as a video.

Return type:

matplotlib.animation.FuncAnimation

Notes

  • In 1D, only one spatial and one frequency variable are used.

  • Complex-valued Hamiltonian fields are truncated to their real parts for integration.

  • Trajectories are shown with both instantaneous position (dot) and full path (dashed line).

apply(u, x_grid, kx, boundary_condition='periodic', y_grid=None, ky=None, dealiasing_mask=None, freq_window='gaussian', clamp=1000000.0, space_window=False)[source]

Apply the pseudo-differential operator to the input field u.

This method dispatches the application of the pseudo-differential operator based on:

  • Whether the symbol is spatially dependent (x/y)

  • The boundary condition in use (periodic or dirichlet)

Supported operations:

  • Constant-coefficient symbols: applied via Fourier multiplication.

  • Spatially varying symbols: applied via Kohn–Nirenberg quantization.

  • Dirichlet boundary conditions: handled with non-periodic convolution-like quantization.

Dispatch Logic:

if not self.is_spatial: u ↦ Op(p)(D) ⋅ u = 𝓕⁻¹[ p(ξ) ⋅ 𝓕(u) ]

elif periodic: u ↦ Op(p)(x,D) ⋅ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ based of FFT (quicker)

elif dirichlet: u ↦ Op(p)(x,D) ⋅ u ≈ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ (slower)

Parameters:

  • u (ndarray) – Function to which the operator is applied

  • x_grid (ndarray) – Spatial grid in x direction

  • kx (ndarray) – Frequency grid in x direction

  • boundary_condition (str) – ‘periodic’ or ‘dirichlet’

  • y_grid (ndarray, optional) – Spatial grid in y direction (for 2D)

  • ky (ndarray, optional) – Frequency grid in y direction (for 2D)

  • dealiasing_mask (ndarray, optional) – Dealiasing mask

  • freq_window (str) – Frequency windowing (‘gaussian’ or ‘hann’)

  • clamp (float) – Clamp symbol values to [-clamp, clamp]

  • space_window (bool) – Apply spatial windowing

Returns:

Result of applying the operator

Return type:

ndarray

asymptotic_expansion(order=3)[source]

Compute the asymptotic expansion of the symbol as |ξ| → ∞ (high-frequency regime).

This method expands the pseudo-differential symbol in inverse powers of the frequency variable(s), either in 1D or 2D. It handles both polynomial and exponential symbols by performing a series expansion in 1/|ξ| up to the specified order.

The expansion is performed directly in Cartesian coordinates for 1D symbols. For 2D symbols, the method uses polar coordinates (ρ, θ) to perform the expansion at infinity in ρ, then converts the result back to Cartesian coordinates.

Parameters:

order (int, optional) – Maximum order of the asymptotic expansion. Default is 3.

Returns:

  • sympy.Expr – The asymptotic expansion of the symbol up to the given order, expressed in Cartesian coordinates. If expansion fails, returns the original unexpanded symbol.

  • Notes

  • - In 1D (expansion is performed directly in terms of ξ.)

  • - In 2D (the symbol is first rewritten in polar coordinates (ρ,θ), expanded asymptotically) – in ρ → ∞, then converted back to Cartesian coordinates (ξ,η).

  • - Handles special case when the symbol is an exponential function by expanding its argument.

    • Symbolic normalization is applied early (via simplify) for 2D expressions to improve convergence.

  • - Robust to failures (catches exceptions and issues warnings instead of raising errors.)

    • Final expression is simplified using powdenest and expand for improved readability.

clear_cache()[source]

Clear cached symbol evaluations.

commutator_symbolic(other, order=1, mode='kn', sign_convention=None)[source]

Compute the symbolic commutator [A, B] = A∘B − B∘A of two pseudo-differential operators using formal asymptotic expansion of their composition symbols.

This method computes the asymptotic expansion of the commutator’s symbol up to a given order, based on the symbolic calculus of pseudo-differential operators in the Kohn–Nirenberg quantization. The result is a purely symbolic sympy expression that captures the leading-order noncommutativity of the operators.

Parameters:

  • other (PseudoDifferentialOperator) – The pseudo-differential operator B to commute with this operator A.

  • order (int, default=1) – Maximum order of the asymptotic expansion. - order=1 yields the leading term proportional to the Poisson bracket {p, q}. - Higher orders include correction terms involving higher mixed derivatives.

Returns:

Symbolic expression for the asymptotic expansion of the commutator symbol σ([A,B]) = σ(A∘B − B∘A).

Return type:

sympy.Expr

compose_asymptotic(other, order=1, mode='kn', sign_convention=None)[source]

Compose two pseudo-differential operators using an asymptotic expansion in the chosen quantization scheme (Kohn–Nirenberg or Weyl).

Parameters:

  • other (PseudoDifferentialOperator) – The operator to compose with this one.

  • order (int, default=1) – Maximum order of the asymptotic expansion.

  • mode ({'kn', 'weyl'}, default='kn') – Quantization mode: - ‘kn’ : Kohn–Nirenberg quantization (left-quantized) - ‘weyl’ : Weyl symmetric quantization

  • sign_convention ({'standard', 'inverse'}, optional) – Controls the phase factor convention for the KN case: - ‘standard’ → (i)^(-n), gives [x, ξ] = +i (physics convention) - ‘inverse’ → (i)^(+n), gives [x, ξ] = -i (mathematical adjoint convention) If None, defaults to ‘standard’.

Returns:

Symbolic expression for the composed symbol up to the given order.

Return type:

sympy.Expr

Notes

  • In 1D (Kohn–Nirenberg):

    (p ∘ q)(x, ξ) ~ Σₙ (1/n!) (i sgn)^n ∂_ξⁿ p(x, ξ) ∂_xⁿ q(x, ξ)

  • In 1D (Weyl):

    (p # q)(x, ξ) = exp[(i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q)] p(x, ξ) q(x, ξ) truncated at given order.

Examples

X = a*x, Y = b*ξ X_op.compose_asymptotic(Y_op, order=3, mode=’weyl’)

evaluate(X, Y, KX, KY, cache=True)[source]

Evaluate the pseudo-differential operator’s symbol on a grid of spatial and frequency coordinates.

The method dynamically selects between 1D and 2D evaluation based on the spatial dimension. If caching is enabled and a cached symbol exists, it returns the cached result to avoid recomputation.

Parameters:

  • X (ndarray) – Spatial grid coordinates. In 1D, Y is ignored.

  • Y (ndarray) – Spatial grid coordinates. In 1D, Y is ignored.

  • KX (ndarray) – Frequency grid coordinates. In 1D, KY is ignored.

  • KY (ndarray) – Frequency grid coordinates. In 1D, KY is ignored.

  • cache (bool, default=True) – If True, stores the computed symbol for reuse in subsequent calls to avoid redundant computation.

Returns:

Evaluated symbol values over the input grid. Shape matches the input spatial/frequency grids.

Return type:

ndarray

Raises:

NotImplementedError – If the spatial dimension is not 1D or 2D.

exponential_symbol(t=1.0, order=1, mode='kn', sign_convention=None)[source]

Compute the symbol of exp(tP) using asymptotic expansion methods.

This method calculates the exponential of a pseudo-differential operator using either a direct power series expansion or a Magnus expansion, depending on the structure of the symbol. The result is valid up to the specified asymptotic order.

Parameters:

  • t (float or sympy.Symbol, default=1.0) – Time or evolution parameter. Common uses: - t = -i*τ for Schrödinger evolution: exp(-iτH) - t = τ for heat/diffusion: exp(τΔ) - t for general propagators

  • order (int, default=3) – Maximum order of the asymptotic expansion. Higher orders include more composition terms, improving accuracy for small t or when non-commutativity effects are significant.

Returns:

Symbolic expression for the exponential operator symbol, computed as an asymptotic series up to the specified order.

Return type:

sympy.Expr

Notes

  • For commutative symbols (e.g., pure multiplication operators), the exponential is exact: exp(tP) = exp(t*p(x,ξ)).

  • For general non-commutative operators, the method uses the BCH-type expansion via iterated composition: exp(tP) ~ I + tP + (t²/2!)P∘P + (t³/3!)P∘P∘P + …

  • Each power P^n is computed via compose_asymptotic, which accounts for the non-commutativity through derivative terms.

  • The expansion is valid for |t| small enough or when the symbol has appropriate decay/growth properties.

  • In quantum mechanics (Schrödinger): U(t) = exp(-itH/ℏ) represents the time evolution operator.

  • In parabolic PDEs (heat equation): exp(tΔ) is the heat kernel.

formal_adjoint()[source]

Compute the formal adjoint symbol P* of the pseudo-differential operator.

The adjoint is defined such that for any test functions u and v, ⟨P u, v⟩ = ⟨u, P* v⟩ holds in the distributional sense. This is obtained by taking the complex conjugate of the symbol and expanding it asymptotically at infinity to ensure proper behavior under integration by parts.

Returns:

  • sympy.Expr – The adjoint symbol P*(x, ξ) in 1D or P*(x, y, ξ, η) in 2D.

  • Notes

  • - In 1D, the expansion is performed in powers of 1/|ξ|.

  • - In 2D, the expansion is radial in |ξ| = sqrt(ξ² + η²).

  • - This method ensures symbolic simplifications for readability and efficiency.

group_velocity_field(xlim=(-2, 2), klim=(-10, 10), density=30)[source]

Plot the group velocity field ∇_ξ p(x, ξ) for 1D pseudo-differential operators.

The group velocity represents the speed at which waves of different frequencies propagate in a dispersive medium. It is defined as the gradient of the symbol p(x, ξ) with respect to the frequency variable ξ.

Parameters:

  • xlim (tuple of float) – Spatial domain limits (x-axis).

  • klim (tuple of float) – Frequency domain limits (ξ-axis).

  • density (int) – Number of grid points per axis used for visualization.

Raises:

NotImplementedError – If called on a 2D operator, since this visualization is only implemented for 1D.

Notes

  • This method visualizes the vector field (∂p/∂ξ) in phase space.

  • Used for analyzing wave propagation properties and dispersion relations.

  • Requires symbolic expression self.expr depending on x and ξ.

interactive_symbol_analysis(xlim=(-2, 2), ylim=(-2, 2), xi_range=(0.1, 5), eta_range=(-5, 5), density=50)[source]

Launch an interactive dashboard for symbol exploration using ipywidgets.

This function provides a user-friendly interface to visualize various aspects of the pseudo-differential operator’s symbol. It supports multiple visualization modes in both 1D and 2D, including group velocity fields, micro-support estimates, symplectic vector fields, symbol amplitude/phase, cotangent fiber structure, characteristic sets and Hamiltonian flows.

Parameters:

  • pseudo_op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed interactively.

  • xlim (tuple of float) – Spatial domain limits along x and y axes respectively.

  • ylim (tuple of float) – Spatial domain limits along x and y axes respectively.

  • xi_range (tuple) – Frequency domain limits along ξ and η axes respectively.

  • eta_range (tuple) – Frequency domain limits along ξ and η axes respectively.

  • density (int) – Number of points per axis used to construct the evaluation grid. Controls resolution.

Notes

  • In 1D mode, sliders control the fixed frequency (ξ₀) and spatial position (x₀).

  • In 2D mode, additional sliders control the second frequency component (η₀) and second spatial coordinate (y₀).

  • Visualization updates dynamically as parameters are adjusted via sliders or dropdown menus.

  • Supported visualization modes:

    ‘Symbol Amplitude’ : |p(x,ξ)| or |p(x,y,ξ,η)| ‘Symbol Phase’ : arg(p(x,ξ)) or similar in 2D ‘Micro-Support (1/|p|)’ : Reciprocal of symbol magnitude ‘Cotangent Fiber’ : Structure of symbol over frequency space at fixed x ‘Characteristic Set’ : Zero set approximation {p ≈ 0} ‘Characteristic Gradient’ : |∇p(x, ξ)| or |∇p(x₀, y₀, ξ, η)| ‘Group Velocity Field’ : ∇_ξ p(x,ξ) or ∇_{ξ,η} p(x,y,ξ,η) ‘Symplectic Vector Field’ : (∇_ξ p, -∇_x p) or similar in 2D ‘Hamiltonian Flow’ : Trajectories generated by the Hamiltonian vector field

Raises:

  • NotImplementedError – If the spatial dimension is not 1D or 2D.

  • Prints –

  • ------ –

  • Interactive matplotlib figures with dynamic updates based on widget inputs. –

is_elliptic_numerically(x_grid, xi_grid, threshold=1e-08)[source]

Check if the pseudo-differential symbol p(x, ξ) is elliptic over a given grid.

A symbol is considered elliptic if its magnitude |p(x, ξ)| remains bounded away from zero across all points in the spatial-frequency domain. This method evaluates the symbol on a grid of spatial and frequency coordinates and checks whether its minimum absolute value exceeds a specified threshold.

Resampling is applied to large grids to prevent excessive memory usage, particularly in 2D.

Parameters:

  • x_grid (ndarray) – Spatial grid: either a 1D array (x) or a tuple of two 1D arrays (x, y).

  • xi_grid (ndarray) – Frequency grid: either a 1D array (ξ) or a tuple of two 1D arrays (ξ, η).

  • threshold (float, optional) – Minimum acceptable value for |p(x, ξ)|. If the smallest evaluated symbol value falls below this, the symbol is not considered elliptic.

Returns:

True if the symbol is elliptic on the resampled grid, False otherwise.

Return type:

bool

is_homogeneous(tol=1e-10)[source]

Check whether the symbol is homogeneous in the frequency variables.

Returns:

Tuple (is_homogeneous, degree) where: - is_homogeneous: True if the symbol satisfies p(λξ, λη) = λ^m * p(ξ, η) - degree: the detected degree m if homogeneous, or None

Return type:

(bool, Rational or float or None)

is_self_adjoint(tol=1e-10)[source]

Check whether the pseudo-differential operator is formally self-adjoint (Hermitian).

A self-adjoint operator satisfies P = P*, where P* is the formal adjoint of P. This property is essential for ensuring real-valued eigenvalues and stable evolution in quantum mechanics and symmetric wave propagation.

Parameters:

tol (float) – Tolerance for symbolic comparison between P and P*. Small numerical differences below this threshold are considered equal.

Returns:

  • bool – True if the symbol p(x, ξ) equals its formal adjoint p*(x, ξ) within the given tolerance, indicating that the operator is self-adjoint.

  • Notes

  • - The formal adjoint is computed via conjugation and asymptotic expansion at infinity in ξ.

  • - Symbolic simplification is used to verify equality, ensuring robustness against superficial – expression differences.

left_inverse_asymptotic(order=1)[source]

Construct a formal left inverse L such that the composition L ∘ P equals the identity operator up to terms of order ξ^{-order}. This expansion is performed asymptotically at infinity in the frequency variable(s).

The left inverse is built iteratively using symbolic differentiation and the method of asymptotic expansions for pseudo-differential operators. It ensures that:

L(P(x,ξ),x,D) ∘ P(x,D) = Id + smoothing operator of order -order

Parameters:

order (int, optional) – Maximum number of terms in the asymptotic expansion (default is 1). Higher values yield more accurate inverses at the cost of increased computational complexity.

Returns:

Symbolic expression representing the principal symbol of the formal left inverse operator L(x,ξ). This expression depends on spatial variables and frequencies, and includes correction terms up to the specified order.

Return type:

sympy.Expr

Notes

  • In 1D: Uses recursive application of the Leibniz formula for symbols.

  • In 2D: Generalizes to multi-indices for mixed derivatives in (x,y) and (ξ,η).

  • Each term involves combinations of derivatives of the original symbol p(x,ξ) and previously computed terms of the inverse.

  • Coefficients include powers of 1j (i) and factorial normalization for derivative terms.

plot_hamiltonian_flow(x0=0.0, xi0=5.0, y0=0.0, eta0=0.0, tmax=1.0, n_steps=100, show_field=True)[source]

Integrate and plot the Hamiltonian trajectories of the symbol in phase space.

This method numerically integrates the Hamiltonian vector field derived from the operator’s symbol to visualize how singularities propagate under the flow. It supports both 1D and 2D problems.

Parameters:

  • x0 (float) – Initial position and frequency (momentum) in 1D.

  • xi0 (float) – Initial position and frequency (momentum) in 1D.

  • y0 (float, optional) – Initial position and frequency in 2D; defaults to zero.

  • eta0 (float, optional) – Initial position and frequency in 2D; defaults to zero.

  • tmax (float) – Final integration time for the ODE solver.

  • n_steps (int) – Number of time steps used in the integration.

Notes

  • The Hamiltonian vector field is obtained from the symplectic flow of the symbol.

  • If the field is complex-valued, only its real part is used for integration.

  • In 1D, the trajectory is plotted in (x, ξ) phase space.

  • In 2D, the spatial trajectory (x(t), y(t)) is shown along with instantaneous momentum vectors (ξ(t), η(t)) using a quiver plot.

Raises:

  • NotImplementedError – If the spatial dimension is not 1D or 2D.

  • Displays –

  • -------- –

  • matplotlib plot – Phase space trajectory(ies) showing the evolution of position and momentum under the Hamiltonian dynamics.

plot_symplectic_vector_field(xlim=(-2, 2), klim=(-5, 5), density=30)[source]

Visualize the symplectic vector field (Hamiltonian vector field) associated with the operator’s symbol.

The plotted vector field corresponds to (∂_ξ p, -∂_x p), where p(x, ξ) is the principal symbol of the pseudo-differential operator. This field governs the bicharacteristic flow in phase space.

Parameters:

  • xlim (tuple of float) – Range for spatial variable x, as (x_min, x_max).

  • klim (tuple of float) – Range for frequency variable ξ, as (ξ_min, ξ_max).

  • density (int) – Number of grid points per axis for the visualization grid.

Raises:

NotImplementedError – If called on a 2D operator (currently only 1D implementation available).

Notes

  • Only supports one-dimensional operators.

  • Uses symbolic differentiation to compute ∂_ξ p and ∂_x p.

  • Numerical evaluation is done via lambdify with NumPy backend.

  • Visualization uses matplotlib quiver plot to show vector directions.

principal_symbol(order=1)[source]

Compute the leading homogeneous component of the pseudo-differential symbol.

This method extracts the principal part of the symbol, which is the dominant term under high-frequency asymptotics (|ξ| → ∞). The expansion is performed in polar coordinates for 2D symbols to maintain rotational symmetry, then converted back to Cartesian form.

Parameters:

order (int) – Order of the asymptotic expansion in powers of 1/ρ, where ρ = |ξ| in 1D or ρ = sqrt(ξ² + η²) in 2D. Only the leading-order term is returned.

Returns:

  • sympy.Expr – The principal symbol component, homogeneous of degree m - order, where m is the original symbol’s order.

  • Notes

  • - In 1D, uses direct series expansion in ξ.

  • - In 2D, expands in radial variable ρ while preserving angular dependence.

  • - Useful for microlocal analysis and constructing parametrices.

pseudospectrum_analysis(x_grid, lambda_real_range, lambda_imag_range, epsilon_levels=[0.1, 0.01, 0.001, 0.0001], resolution=100, method='spectral', L=None, N=None, use_sparse=False, parallel=True, n_workers=4, adaptive=False, adaptive_threshold=0.5, auto_range=True, plot=True)[source]

Compute and visualize the pseudospectrum of the operator.

Optimizations: - Uses apply() method instead of manual loops - Parallel computation of resolvent norms - Sparse matrix support for large N - Optional adaptive grid refinement

Parameters:

  • x_grid (array) – Spatial grid for quantization

  • lambda_real_range (tuple) – (min, max) for real part of λ

  • lambda_imag_range (tuple) – (min, max) for imaginary part of λ

  • epsilon_levels (list) – Levels for ε-pseudospectrum contours

  • resolution (int) – Grid resolution for λ sampling

  • method (str) – ‘spectral’ or ‘finite_difference’

  • L (float, optional) – Domain half-length for spectral method

  • N (int, optional) – Number of grid points

  • use_sparse (bool) – Use sparse matrices for large N

  • parallel (bool) – Enable parallel computation

  • n_workers (int) – Number of parallel workers

  • adaptive (bool) – Use adaptive grid refinement

  • adaptive_threshold (float) – Threshold for adaptive refinement

Returns:

Dictionary with pseudospectrum data and operator matrix

Return type:

dict

right_inverse_asymptotic(order=1)[source]

Construct a formal right inverse R of the pseudo-differential operator P such that the composition P ∘ R equals the identity plus a smoothing operator of order -order.

This method computes an asymptotic expansion for the right inverse using recursive corrections based on derivatives of the symbol p(x, ξ) and lower-order terms of R.

Parameters:

order (int) – Number of terms to include in the asymptotic expansion. Higher values improve approximation at the cost of complexity and computational effort.

Returns:

The symbolic expression representing the formal right inverse R(x, ξ), which satisfies: P ∘ R = Id + O(⟨ξ⟩^{-order}), where ⟨ξ⟩ = (1 + |ξ|²)^{1/2}.

Return type:

sympy.Expr

Notes

  • In 1D: The recursion involves spatial derivatives of R and derivatives of p with respect to ξ.

  • In 2D: The multi-index generalization is used with mixed derivatives in ξ and η.

  • The construction relies on the non-vanishing of the principal symbol p to ensure invertibility.

  • Each term in the expansion corresponds to higher-order corrections involving commutators between the operator P and the current approximation of R.

symbol_order(max_order=10, tol=0.001)[source]

Estimate the homogeneity order of the pseudo-differential symbol in high-frequency asymptotics.

This method attempts to determine the leading-order behavior of the symbol p(x, ξ) or p(x, y, ξ, η) as |ξ| → ∞ (in 1D) or |(ξ, η)| → ∞ (in 2D). The returned value represents the asymptotic growth or decay rate, which is essential for understanding the regularity and mapping properties of the corresponding operator.

The function uses symbolic preprocessing to ensure proper factorization of frequency variables, especially in sqrt and power expressions, to avoid erroneous order detection (e.g., due to hidden scaling).

Parameters:

  • max_order (int, optional) – Maximum number of terms to consider in the series expansion. Default is 10.

  • tol (float, optional) – Tolerance threshold for evaluating the coefficient magnitude. If the coefficient is too small, the detected order may be discarded. Default is 1e-3.

Returns:

  • If the symbol is homogeneous, returns its exact homogeneity degree as a float.

  • Otherwise, estimates the dominant asymptotic order from leading terms in the expansion.

  • Returns None if no valid order could be determined.

Return type:

float or None

Notes

  • In 1D:
    Two strategies are used:

    1. Expand directly in xi at infinity.

    2. Substitute xi = 1/z and expand around z = 0.

  • In 2D:

    • Transform the symbol into polar coordinates: (xi, eta) = rho*(cos(theta), sin(theta)).

    • Expand in rho at infinity, then extract the leading term’s power.

    • An alternative substitution using 1/z is also tried if the first method fails.

  • Preprocessing steps:

    • Sqrt expressions involving frequencies are rewritten to isolate the leading variable.

    • Power expressions are factored explicitly to ensure correct symbolic scaling.

  • If the symbol is not homogeneous, a warning is issued, and the result should be interpreted with care.

  • For non-homogeneous symbols, only the principal asymptotic term is considered.

Raises:

NotImplementedError – If the spatial dimension is neither 1 nor 2.

symplectic_flow()[source]

Compute the Hamiltonian vector field associated with the principal symbol.

This method derives the canonical equations of motion for the phase space variables (x, ξ) in 1D or (x, y, ξ, η) in 2D, based on the Hamiltonian formalism. These describe how position and frequency variables evolve under the flow generated by the symbol.

Returns:

A dictionary containing the components of the Hamiltonian vector field: - In 1D: keys are ‘dx/dt’ and ‘dxi/dt’, corresponding to dx/dt = ∂p/∂ξ and dξ/dt = -∂p/∂x. - In 2D: keys are ‘dx/dt’, ‘dy/dt’, ‘dxi/dt’, and ‘deta/dt’, with similar definitions:

dx/dt = ∂p/∂ξ, dy/dt = ∂p/∂η, dξ/dt = -∂p/∂x, dη/dt = -∂p/∂y.

Return type:

dict

Notes

  • The Hamiltonian here is the principal symbol p(x, ξ) itself.

  • This flow preserves the symplectic structure of phase space.

trace_formula(volume_element=None, numerical=False, x_bounds=None, xi_bounds=None)[source]

Compute the semiclassical trace of the pseudo-differential operator.

The trace formula relates the quantum trace of an operator to a phase-space integral of its symbol, providing a fundamental link between classical and quantum mechanics. This implementation supports both symbolic and numerical integration.

Parameters:

  • volume_element (sympy.Expr, optional) – Custom volume element for the phase space integration. If None, uses the standard Liouville measure dx dξ/(2π)^d.

  • numerical (bool, default=False) – If True, perform numerical integration over specified bounds. If False, attempt symbolic integration (may fail for complex symbols).

  • x_bounds (tuple of tuples, optional) – Spatial integration bounds. For 1D: ((x_min, x_max),) For 2D: ((x_min, x_max), (y_min, y_max)) Required if numerical=True.

  • xi_bounds (tuple of tuples, optional) – Frequency integration bounds. For 1D: ((xi_min, xi_max),) For 2D: ((xi_min, xi_max), (eta_min, eta_max)) Required if numerical=True.

Returns:

The trace of the operator. Returns a symbolic expression if numerical=False, or a float if numerical=True.

Return type:

sympy.Expr or float

Notes

  • The semiclassical trace formula states: Tr(P) = (2π)^{-d} ∫∫ p(x,ξ) dx dξ where d is the spatial dimension and p(x,ξ) is the operator symbol.

  • For 1D: Tr(P) = (1/2π) ∫_{-∞}^{∞} ∫_{-∞}^{∞} p(x,ξ) dx dξ

  • For 2D: Tr(P) = (1/4π²) ∫∫∫∫ p(x,y,ξ,η) dx dy dξ dη

  • This formula is exact for trace-class operators and provides an asymptotic approximation for general pseudo-differential operators.

  • Physical interpretation: the trace counts the “number of states” weighted by the observable p(x,ξ).

  • For projection operators (χ_Ω with χ² = χ), the trace gives the dimension of the range, related to the phase space volume of Ω.

  • The factor (2π)^{-d} comes from the quantum normalization of coherent states / Weyl quantization.

visualize_characteristic_gradient(x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0)[source]

Visualize the norm of the gradient of the symbol in phase space.

This method computes the magnitude of the gradient |∇p| of a pseudo-differential symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. The resulting colormap reveals regions where the symbol varies rapidly or remains nearly stationary, which is particularly useful for analyzing characteristic sets.

Parameters:

  • x_grid (numpy.ndarray) – 1D array of spatial coordinates for the x-direction.

  • xi_grid (numpy.ndarray) – 1D array of frequency coordinates (ξ).

  • y_grid (numpy.ndarray, optional) – 1D array of spatial coordinates for the y-direction (used in 2D mode). Default is None.

  • eta_grid (numpy.ndarray, optional) – 1D array of frequency coordinates (η) for the 2D case. Default is None.

  • x0 (float, optional) – Fixed x-coordinate for evaluating the symbol in 2D. Default is 0.0.

  • y0 (float, optional) – Fixed y-coordinate for evaluating the symbol in 2D. Default is 0.0.

Returns:

Displays a 2D colormap of |∇p| over the relevant phase-space domain.

Return type:

None

Notes

  • In 1D, the full gradient ∇p = (∂ₓp, ∂ξp) is computed over the (x, ξ) grid.

  • In 2D, the gradient ∇p = (∂ξp, ∂ηp) is computed at a fixed spatial point (x₀, y₀) over the (ξ, η) grid.

  • Numerical differentiation is performed using np.gradient.

  • High values of |∇p| indicate rapid variation of the symbol, while low values typically suggest characteristic regions.

visualize_characteristic_set(x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0, levels=[0.1])[source]

Visualize the characteristic set of the pseudo-differential symbol, defined as the approximate zero set p(x, ξ) ≈ 0.

In microlocal analysis, the characteristic set is the locus of points in phase space (x, ξ) where the symbol p(x, ξ) vanishes, playing a key role in understanding propagation of singularities.

Parameters:

  • x_grid (ndarray) – Spatial grid values (1D array) for plotting in 1D or evaluation point in 2D.

  • xi_grid (ndarray) – Frequency variable grid values (1D array) used to construct the frequency domain.

  • x0 (float, optional) – Fixed spatial coordinate in 2D case for evaluating the symbol at a specific x position.

  • y0 (float, optional) – Fixed spatial coordinate in 2D case for evaluating the symbol at a specific y position.

Notes

  • For 1D, this method plots the contour of |p(x, ξ)| = ε with ε = 1e-5 over the (x, ξ) plane.

  • For 2D, it evaluates the symbol at fixed (x₀, y₀) and plots the characteristic set in the (ξ, η) frequency plane.

  • This visualization helps identify directions of degeneracy or hypoellipticity of the operator.

Raises:

  • NotImplementedError – If called on a solver with dimensionality other than 1D or 2D.

  • Displays –

  • ------ –

  • A matplotlib contour plot showing either: –

    • The characteristic curve in the (x, ξ) phase plane (1D), - The characteristic surface slice in the (ξ, η) frequency plane at (x₀, y₀) (2D).

visualize_fiber(x_grid, xi_grid, x0=0.0, y0=0.0)[source]

Plot the cotangent fiber structure at a fixed spatial point (x₀[, y₀]).

This visualization shows how the symbol p(x, ξ) behaves on the cotangent fiber above a fixed spatial point. In microlocal analysis, this provides insight into the frequency content of the operator at that location.

Parameters:

  • x_grid (ndarray) – Spatial grid values (1D) for evaluation in 1D case.

  • xi_grid (ndarray) – Frequency grid values (1D) for evaluation in both 1D and 2D cases.

  • x0 (float, optional) – Fixed x-coordinate of the base point in space (1D or 2D).

  • y0 (float, optional) – Fixed y-coordinate of the base point in space (2D only).

Notes

  • In 1D: Displays |p(x, ξ)| over the (x, ξ) phase plane near the fixed point.

  • In 2D: Fixes (x₀, y₀) and evaluates p(x₀, y₀, ξ, η), showing the fiber over that point.

  • The color map represents the magnitude of the symbol, highlighting regions where it vanishes or becomes singular.

Raises:

NotImplementedError – If called in 2D with missing or improperly formatted grids.

visualize_micro_support(xlim=(-2, 2), klim=(-10, 10), threshold=0.001, density=300)[source]

Visualize the micro-support of the operator by plotting the inverse of the symbol magnitude 1 / |p(x, ξ)|.

The micro-support provides insight into the singularities of a pseudo-differential operator in phase space (x, ξ). Regions where |p(x, ξ)| is small correspond to large values in 1/|p(x, ξ)|, highlighting areas of significant operator influence or singularity.

Parameters:

  • xlim (tuple) – Spatial domain limits (x_min, x_max).

  • klim (tuple) – Frequency domain limits (ξ_min, ξ_max).

  • threshold (float) – Threshold below which |p(x, ξ)| is considered effectively zero; used for numerical stability.

  • density (int) – Number of grid points along each axis for visualization resolution.

Raises:

NotImplementedError – If called on a solver with dimension greater than 1 (only 1D visualization is supported).

Notes

  • This method evaluates the symbol p(x, ξ) over a grid and plots its reciprocal to emphasize regions where the symbol is near zero.

  • A small constant (1e-10) is added to the denominator to avoid division by zero.

  • The resulting plot helps identify characteristic sets.

visualize_phase(x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0)[source]

Plot the phase (argument) of the pseudodifferential operator’s symbol p(x, ξ) or p(x, y, ξ, η).

This visualization helps in understanding the oscillatory behavior and regularity properties of the operator in phase space. The phase is displayed modulo 2π using a cyclic colormap (‘twilight’) to emphasize its periodic nature.

Parameters:

  • x_grid (ndarray) – 1D array of spatial coordinates (x).

  • xi_grid (ndarray) – 1D array of frequency coordinates (ξ).

  • y_grid (ndarray, optional) – 2D spatial grid for y-coordinate (in 2D problems). Default is None.

  • eta_grid (ndarray, optional) – 2D frequency grid for η (in 2D problems). Not used directly but kept for API consistency.

  • xi0 (float, optional) – Fixed value of ξ for slicing in 2D visualization. Default is 0.0.

  • eta0 (float, optional) – Fixed value of η for slicing in 2D visualization. Default is 0.0.

  • Notes

  • 1D (- In)

  • 2D (- In)

  • π. (- Uses plt.pcolormesh with 'twilight' colormap to represent angles from -π to)

  • Raises

  • NotImplementedError (-)

visualize_symbol_amplitude(x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0)[source]

Display the modulus |p(x, ξ)| or |p(x, y, ξ₀, η₀)| as a color map.

This method visualizes the amplitude of the pseudodifferential operator’s symbol in either 1D or 2D spatial configuration. In 2D, the frequency variables are fixed to specified values (ξ₀, η₀) for visualization purposes.

Parameters:

  • x_grid (ndarray) – Spatial grids over which to evaluate the symbol. y_grid is optional and used only in 2D.

  • y_grid (ndarray) – Spatial grids over which to evaluate the symbol. y_grid is optional and used only in 2D.

  • xi_grid (ndarray) – Frequency grids. In 2D, these define the domain over which the symbol is evaluated, but the visualization fixes ξ = ξ₀ and η = η₀.

  • eta_grid (ndarray) – Frequency grids. In 2D, these define the domain over which the symbol is evaluated, but the visualization fixes ξ = ξ₀ and η = η₀.

  • xi0 (float, optional) – Fixed frequency values for slicing in 2D visualization. Defaults to zero.

  • eta0 (float, optional) – Fixed frequency values for slicing in 2D visualization. Defaults to zero.

Notes

  • In 1D: Visualizes |p(x, ξ)| over the (x, ξ) grid.

  • In 2D: Visualizes |p(x, y, ξ₀, η₀)| at fixed frequencies ξ₀ and η₀.

  • The color intensity represents the magnitude of the symbol, highlighting regions where the symbol is large or small.

psiop.freq_window_2d(P, KXb, KYb, kx_shift, ky_shift, mode)[source]

Apply a 2-D frequency window in-place and return P.

psiop.invalidate_kn_cache()[source]

Clear the phase-matrix cache.

Call this after any solver.setup() call that changes Lx, Nx, or the frequency grid, so the next apply re-builds the matrices for the new grid.

psiop.kohn_nirenberg_fft(u_vals, symbol_func, x_grid, kx, fft_func, ifft_func, dim=1, y_grid=None, ky=None, freq_window='gaussian', clamp=1000000.0, space_window=False)[source]

Numerically stable Kohn–Nirenberg quantization of a pseudo-differential operator.

Applies the pseudo-differential operator Op(p) to the function f via the Kohn–Nirenberg quantization:

[Op(p)f](x) = (1/(2π)^d) ∫ p(x, ξ) e^{ix·ξ} ℱ[f](ξ) dξ

where p(x, ξ) is a symbol that may depend on both spatial variables x and frequency variables ξ.

This method supports both 1D and 2D cases and includes optional smoothing techniques to improve numerical stability.

Parameters:

  • u_vals (ndarray) – Spatial samples of the input function f(x) or f(x, y), defined on a uniform grid.

  • symbol_func (callable) – A function representing the full symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. Must accept NumPy-compatible array inputs and return a complex-valued array.

  • freq_window (str) – Type of frequency-domain window to apply: - ‘gaussian’: smooth decay near high frequencies - ‘hann’: cosine-based tapering with hard cutoff - None: no frequency window applied

  • clamp (float) – Upper bound on the absolute value of the symbol. Prevents numerical blow-up from large values.

  • space_window (bool) – Whether to apply a spatial Gaussian window to suppress edge effects in physical space.

  • u_vals – Input function values

  • symbol_func – Symbol function p(x, ξ) or p(x, y, ξ, η)

  • x_grid (ndarray) – Spatial grid in x direction

  • kx (ndarray) – Frequency grid in x direction

  • fft_func (callable) – FFT function (fft or fft2)

  • ifft_func (callable) – Inverse FFT function

  • dim (int) – Dimension (1 or 2)

  • y_grid (ndarray, optional) – Spatial grid in y direction (for 2D)

  • ky (ndarray, optional) – Frequency grid in y direction (for 2D)

  • freq_window – Windowing function (‘gaussian’ or ‘hann’)

  • clamp – Clamp symbol values to [-clamp, clamp]

  • space_window – Apply spatial windowing

Returns:

  • ndarray – Result of applying the operator

  • 2-D optimizations

  • —————–

  • * The 4-D tensor (Nx, Ny, Nkx, Nky) is never fully materialized. – The x-axis is sliced into row-blocks; for each block only a (block_rows, Ny, Nkx, Nky) array is needed.

  • * Row-blocks are processed in parallel with ThreadPoolExecutor.

  • * Phase kernel exp(i*x*kx) is factored as a pair of 2-D arrays, – exploiting exp(i*(x*kx + y*ky)) = exp(i*x*kx) * exp(i*y*ky).

psiop.kohn_nirenberg_nonperiodic(u_vals, x_grid, xi_grid, symbol_func, freq_window: str = 'gaussian', clamp: float = 1000000.0, space_window: bool = False, _cache: dict = {})[source]

Apply a pseudo-differential operator using Kohn–Nirenberg quantization on a non‑periodic domain.

The operator is defined by

(Op(p) u)(x) = (1/(2π)) ∫ p(x, ξ) û(ξ) e^{i x ξ} dξ,

where û(ξ) is the non‑periodic Fourier transform of u(x) on the given grids.

For 1D inputs the function caches several objects that depend only on the spatial and frequency grids: the forward Fourier kernel, the inverse kernel, and pre‑computed frequency windows and spatial taper. These are reused automatically on subsequent calls with the same grids, reducing setup overhead. The cache can be cleared with invalidate_kn_cache().

For 2D inputs the computation is performed block‑wise along the first spatial dimension to control memory usage, and parallelised with a thread pool.

Parameters:

  • u_vals (ndarray) – Function values at the spatial grid points. One‑dimensional array for 1D, two‑dimensional array for 2D (first dimension corresponds to x₁, second to x₂).

  • x_grid (ndarray or tuple of ndarray) – Spatial grid(s). For 1D: a 1D array of coordinates. For 2D: a tuple (x1_grid, x2_grid) of two 1D arrays.

  • xi_grid (ndarray or tuple of ndarray) – Frequency grid(s). For 1D: a 1D array. For 2D: a tuple (xi1_grid, xi2_grid).

  • symbol_func (callable) –

    The symbol p(x, ξ) of the operator. * In 1D: called as symbol_func(x, xi) where x and xi are broadcastable

    arrays (e.g., x[:, None] and xi[None, :]). The function should support broadcasting to produce an array of shape (len(x), len(xi)).

    • In 2D: called as symbol_func(X1, X2, XI1, XI2) where the arguments are broadcastable arrays; the output shape should be compatible with (Nx1, Nx2, Nxi1, Nxi2).

    Must return an array of complex numbers (or castable to complex128).

  • freq_window ({'gaussian', 'hann'}, default='gaussian') – Frequency‑space window applied to the symbol before the inverse transform. * ‘gaussian’ : exp(‑(ξ/σ_w)⁴) with σ_w = 0.8·max(|ξ|). * ‘hann’ : Hanning window over the full frequency range.

  • clamp (float, default=1e6) – Absolute value beyond which symbol entries are clipped. Helps avoid numerical overflow.

  • space_window (bool, default=False) – If True, multiply the symbol by a Gaussian spatial taper centred in the middle of the spatial domain, with width equal to half the domain length.

  • _cache (dict, optional) – Internal cache dictionary. Do not use directly; exposed only for testing and debugging.

Returns:

Result of applying the pseudo‑differential operator to u_vals. Same shape as u_vals.

Return type:

ndarray

Notes

Caching (1D only)
The cache stores for a given grid pair (x, xi):

  • phase_ft : exp(‑i ξ x) – shape (Nxi, Nx)

  • exp_matrix : exp( i x ξ) – shape (Nx, Nxi)

  • window_gauss : 1D Gaussian window – shape (Nxi,)

  • window_hann : 1D Hann window – shape (Nxi,)

  • spatial_taper: 1D Gaussian taper – shape (Nx,) (if space_window ever used)

The cache key is derived from the first and last points of both grids, so the cache is automatically invalidated when either grid changes in shape or endpoints. A warning is issued the first time a new grid pair is encountered.

2D implementation

To avoid forming a full 4D array, the computation is split into blocks along the first spatial dimension (x₁). Each block is processed independently by a worker thread. The number of workers is taken from solver.FFT_WORKERS. Memory usage scales as O(block_size × Nx₂ × Nxi₁ × Nxi₂).

Examples

1D example: apply the derivative operator i·ξ.

>>> import numpy as np
>>> from kn_operator import kohn_nirenberg_nonperiodic, invalidate_kn_cache
>>> x = np.linspace(-10, 10, 200)
>>> xi = np.fft.fftshift(np.fft.fftfreq(len(x), d=x[1]-x[0])) * 2*np.pi
>>> u = np.exp(-x**2)
>>> def symbol(x, xi): return 1j * xi   # broadcasting works
>>> du = kohn_nirenberg_nonperiodic(u, x, xi, symbol)
>>> # Compare with analytical derivative: -2*x * exp(-x**2)