psiop package

The core symbolic-numerical toolkit for scalar pseudo-differential operators, Peetre decomposition, and asymptotic calculus.

psiop — Symbolic–numerical toolkit for pseudo‑differential operators in 1D/2D

Overview

The psiop package provides a unified framework for constructing, manipulating, and numerically applying pseudo‑differential operators (ΨDOs) in one and two spatial dimensions. It combines symbolic symbol calculus (SymPy) with numerical evaluation, quantization, diagnostics, and visualization (NumPy / SciPy / Matplotlib).

The package is intended for researchers and students working in microlocal analysis, spectral theory, semiclassical analysis, Hamiltonian / geometric optics, and the numerical analysis of PDEs.

Subpackage layout

psiop (this module): Core scalar PseudoDifferentialOperator, Peetre decomposition, asymptotic symbolic calculus, quantization conversion, operator diagnostics, and Hamiltonian-flow / visualization wrappers.
psiop.matpsiop: MatrixPseudoDifferentialOperator — N × N matrix-valued symbols acting on vector fields and matrix-valued fields, with noncommutative asymptotic composition, matrix commutators, and matrix exponential symbols.
psiop.psiop_apply: Standalone numerical kernels: Kohn–Nirenberg FFT / non-periodic quadrature, NUFFT joint-residual backend, AAA rational approximation, and Chebyshev / SVD low-rank factorization.
psiop.psiop_solver: Time-stepping solvers and propagators: build_propagator, PropagatorFamily, solve_first_order, solve_second_order, solve_matrix_field, solve_sylvester_field, solve_ricci_flow_conformal_2d, and grid-generation utilities.

Main object

PseudoDifferentialOperator: Scalar pseudo‑differential operator in 1D or 2D. The symbol may be given explicitly (mode='symbol') or extracted automatically from a differential expression acting on a test function (mode='auto').

Key features

Symbol creation (1D and 2D):: explicit symbol mode p(x, ξ) or p(x, y, ξ, η); automatic extraction from differential operators via plane‑wave testing.
Quantization conventions:: Kohn–Nirenberg (default) and Weyl quantization; asymptotic conversion between the two via the exponential cross-derivative series (exact and finite for polynomial symbols).
High‑frequency symbolic analysis:: principal symbol extraction; asymptotic expansion as |ξ| → ∞; homogeneity detection; symbolic / numerical order estimation.
Asymptotic symbolic calculus:: operator composition p ∘ q (KN and Weyl / Moyal); commutators [P, Q]; formal left and right inverses; formal adjoint P*; fractional powers P^α; exponential symbols exp(tP).
Peetre‑style symbolic decomposition:: local polynomial part in the frequency variables; separable terms a(x) q(ξ); genuinely joint, non‑separable residual; automatic backend selection for the joint residual (NUFFT, AAA, or low-rank) — see psiop.psiop_apply for the numerical kernels.
Operator diagnostics:: symbolic and numerical ellipticity tests; formal self‑adjointness checks; semiclassical trace formula (symbolic or numerical); pseudospectrum computation with resolvent-norm contours, eigenvalue overlays, sparse / dense solvers, and optional adaptive refinement.
Hamiltonian and microlocal dynamics:: Hamiltonian vector field of the principal symbol; symplectic flow visualization; bicharacteristic integration; singularity propagation animation (2D / 3D phase-space projections); characteristic set and gradient visualization.
Visualization:: Symbol amplitude / phase plots, cotangent fibers, micro-support, wavefront set, group velocity fields, and interactive ipywidgets dashboards — all delegated to the microlocal plotting module.

Mathematical background

Symbols and quantization

A pseudo‑differential operator P acting on functions of x ∈ ℝⁿ (n = 1 or 2) is defined by its symbol p(x, ξ) on phase space T*ℝⁿ = ℝⁿ_x × ℝⁿ_ξ. Symbols are classified by their growth in ξ: p ∈ S^m_{1,0} if

|∂_x^β ∂_ξ^α p(x, ξ)| ≤ C_{α,β} (1 + |ξ|)^(m - |α|)

for all multi‑indices α, β. The integer / real number m is the order of the operator.

Quantization maps

Kohn–Nirenberg (left / standard, default):

(P u)(x) = Op^KN(p) u(x) = (2π)⁻ⁿ ∫_{ℝⁿ} exp(i x·ξ) p(x, ξ) û(ξ) dξ

Weyl (symmetric):

Op^w(a) u(x) = (2π)⁻ⁿ ∫∫ exp(i(x-y)·ξ) a((x+y)/2, ξ) u(y) dy dξ

Conversion between the two is realized as an asymptotic series in mixed x–ξ derivatives:

a_KN(x, ξ) = exp(−i/2 · ∂_x ∂_ξ) a_Weyl(x, ξ)
∼ Σ_{k≥0} (−i/2)^k / k! · (∂_x ∂_ξ)^k a_Weyl(x, ξ)

In 2D the cross-derivative operator sums both coordinate pairs: ∂_x ∂_ξ → ∂_x ∂_ξ + ∂_y ∂_η. For polynomial symbols the series is exact and finite; for general S^m_{1,0} symbols it is interpreted asymptotically.

Asymptotic composition

For two symbols p and q, the Kohn–Nirenberg expansion is

(p ∘ q)(x, ξ) ∼ Σ_α i^(-|α|) / α! · ∂_ξ^α p · ∂_x^α q

A Weyl / Moyal star product is also available. The same machinery underlies commutators, formal inverses, adjoints, fractional powers, and exponential symbols.

Principal symbol and order

The principal symbol is the leading homogeneous component as |ξ| → ∞. If p is homogeneous of degree m, p(x, λξ) = λ^m p(x, ξ) for λ > 0.

Peetre decomposition

For numerical application the symbol is split into

p(x, ξ) = p_local(x, ξ) + p_sep(x, ξ) + p_joint(x, ξ)

where p_local is polynomial in ξ (differential part), p_sep = Σ_k a_k(x) q_k(ξ) is a sum of separable Fourier multipliers, and p_joint contains the genuinely entangled space–frequency residual — routed automatically to NUFFT, AAA, or low-rank backends (see psiop.psiop_apply).

Hamiltonian flow

The Hamiltonian vector field of the principal symbol governs singularity propagation along bicharacteristics:

dx/dt = ∂_ξ p, dξ/dt = −∂_x p (1D) dx/dt = ∂_ξ p, dy/dt = ∂_η p, dξ/dt = −∂_x p, dη/dt = −∂_y p (2D)

Pseudospectrum

The ε-pseudospectrum of P is

σ_ε(P) = { λ ∈ ℂ : ‖(P − λI)⁻¹‖ ≥ ε⁻¹ }

determined by the smallest singular value of the shifted matrix discretization: ‖(P − λI)⁻¹‖ = 1 / σ_min(P − λI). Particularly informative for non-normal operators.

References

class psiop.PseudoDifferentialOperator(expr, vars_x, var_u=None, mode='symbol', quantization='kohn-nirenberg', apply_backend='peetre', compute_peetre=False, peetre_options=None)[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. Thin delegate to the module-level animate_singularity engine (previously ~130 lines of duplicated Hamiltonian/ODE setup here, plus a near-identical copy further down the module). cf. package microlocal.py.

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, weyl_order=4, backend=None, apply_joint=True, joint_backend='direct', joint_degree=6, joint_tol=1e-05, joint_bounds=None, joint_max_rel_error=None)[source]

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

Dispatches based on: - Whether the symbol is spatially dependent (x / y). - The boundary condition (‘periodic’ or ‘dirichlet’). - The quantization convention stored in self.quantization.

Supported quantizations

‘kohn-nirenberg’ (default)

Standard KN quantization. The symbol is used as-is.

‘weyl’

Weyl quantization. The Weyl symbol is first converted to its KN equivalent via the asymptotic series

a_KN = exp(+i/2 * d_x d_xi) a_Weyl

truncated at weyl_order. The corrected KN symbol is then passed to the existing KN numerical backend — no new numerical kernel is needed.

Dispatch logic

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)

param u:: Function to which the operator is applied.
type u:: ndarray
param x_grid:: Spatial grid in x direction.
type x_grid:: ndarray
param kx:: Frequency grid in x direction.
type kx:: ndarray
param boundary_condition:: Boundary condition type.
type boundary_condition:: {‘periodic’, ‘dirichlet’}, default=’periodic’
param y_grid:: Spatial grid in y direction (2D only).
type y_grid:: ndarray, optional
param ky:: Frequency grid in y direction (2D only).
type ky:: ndarray, optional
param dealiasing_mask:: Dealiasing mask applied in Fourier space.
type dealiasing_mask:: ndarray, optional
param freq_window:: Frequency-domain smoothing window.
type freq_window:: {‘gaussian’, ‘hann’, None}, default=’gaussian’
param clamp:: Clip symbol values to [-clamp, clamp].
type clamp:: float, default=1e6
param space_window:: Apply a Gaussian spatial taper.
type space_window:: bool, default=False
param weyl_order:: Truncation order for the Weyl -> KN asymptotic correction. Ignored when self.quantization != 'weyl'. The series is exact and finite for polynomial symbols.
type weyl_order:: int, default=4
param backend:: Numerical application backend. If None (default), uses self.apply_backend. ‘peetre’ applies the operator via the local/separable/low-rank Peetre decomposition; ‘direct’ evaluates the Kohn-Nirenberg (or Weyl-corrected) integral pointwise.
type backend:: {‘peetre’, ‘direct’}, optional
param apply_joint:: When backend=’peetre’, whether to include the genuinely joint (non-separable) residual term of the Peetre decomposition. If False, only the local and separable terms are applied.
type apply_joint:: bool, default=True
param joint_backend:: Numerical representation used for the joint residual term when apply_joint=True (e.g. direct evaluation, low-rank/AAA, etc.).
type joint_backend:: str, default=”direct”
param joint_degree:: Polynomial/Chebyshev degree used when approximating the joint residual term.
type joint_degree:: int, default=6
param joint_tol:: Tolerance controlling the accuracy of the joint-residual approximation.
type joint_tol:: float, default=1e-5
param joint_bounds:: Bounding box for the joint residual’s domain of approximation. If None, inferred automatically.
type joint_bounds:: tuple, optional
param joint_max_rel_error:: If given, an upper bound on the acceptable relative error of the joint-residual approximation; the backend may refine or fall back to a more accurate method if this is exceeded.
type joint_max_rel_error:: float, optional
returns:: Result of applying the operator, same shape as u.
rtype:: ndarray
raises ValueError:: If boundary_condition is not ‘periodic’ or ‘dirichlet’.

apply_hybrid(u, x_grid, kx, y_grid=None, ky=None, **kwargs)[source]

Apply the operator via Peetre decomposition with per-term backend routing.

The symbol is split with peetre_decomposition into a separable part and a joint (non-separable) residual. The separable part is applied directly through apply_peetre (apply_joint=False). The joint residual is then expanded and split into its individual additive terms (sympy.Add.make_args on the expanded symbol); each nonzero term is wrapped in its own PseudoDifferentialOperator and applied with joint_backend=’auto’, so _auto_select_joint_backend picks a backend independently for that term — NUFFT if a decomposition exists for it, AAA if it’s a rational function (or has a symbolic polynomial denominator), otherwise low-rank.

Splitting before routing means terms that individually admit a fast decomposition get one, instead of the whole joint residual being treated as a single mixed term and falling back to O(N²) direct quadrature.

Parameters:

u (array_like) – Input field values, matching the shape/dtype convention used by apply/apply_peetre.
x_grid (array_like) – Spatial grid(s). y_grid is only used for 2D operators.
y_grid (array_like) – Spatial grid(s). y_grid is only used for 2D operators.
kx (array_like) – Frequency grid(s) conjugate to x_grid/y_grid. ky is only used for 2D operators.
ky (array_like) – Frequency grid(s) conjugate to x_grid/y_grid. ky is only used for 2D operators.
**kwargs – Forwarded to every apply_peetre call this method makes (both the separable-part call and each per-term joint call) — e.g. boundary_condition, freq_window, clamp, space_window, weyl_order, joint_degree, joint_tol, joint_bounds, joint_max_rel_error, joint_num_samples, joint_seed, use_cache. Do not pass apply_joint or joint_backend here: both are fixed internally by this method, and supplying them will raise TypeError: apply_peetre() got multiple values for keyword argument.

Returns:

Complex128 array, same shape as u: the separable contribution plus the sum of the per-term joint contributions.

Return type:

numpy.ndarray

Notes

Backend selection happens per additive term of the joint symbol, not once for the whole residual, so this can outperform calling apply_peetre directly with joint_backend=’auto’ whenever the joint part mixes terms that individually favor different backends. Compare with apply_hybrid_old, which instead loops over peetre_decomposition’s own joint_residual list without expanding it into additive terms.

apply_peetre(u, x_grid, kx, boundary_condition='periodic', y_grid=None, ky=None, dealiasing_mask=None, freq_window='gaussian', clamp=1000000.0, space_window=False, weyl_order=4, apply_joint=True, decomposition=None, use_cache=True, separable_local=False, joint_backend='direct', joint_degree=6, joint_tol=1e-05, joint_bounds=None, joint_max_rel_error=None, joint_num_samples=10000, joint_seed=42)[source]

Apply the operator by summing the pieces of its Peetre decomposition.

peetre_decomposition splits the symbol into three kinds of contributions, each applied by whichever method suits its structure:

Local terms (local_terms, or the legacy local/local_symbol keys): applied directly via PseudoDifferentialOperator.apply when a term isn’t separable, or as an amplitude-times-operator pair otherwise.
Separable terms a(x) * q(xi) (deco[‘separable’]): q is applied to u as its own operator, then the result is multiplied pointwise by a evaluated on x_grid/y_grid (falling back to applying the full product a * q as one operator if a can’t be lambdified).
Joint residual (deco[‘joint_symbol’], the non-separable remainder): applied via _apply_joint_residual using joint_backend (‘direct’, ‘lowrank’, ‘nufft’, or ‘aaa’), unless apply_joint=False, in which case it is dropped and a warning is issued — the returned result is then only a local+separable (asymptotic) approximation of the true operator.

For Weyl-quantized operators, the Weyl symbol is first converted to its Kohn-Nirenberg equivalent (weyl_to_kn_symbol, truncated at weyl_order) before the decomposition above is computed or used; all downstream sub-operators are built with quantization=’kohn-nirenberg’.

Parameters:

u (array_like) – Input field values on the grid(s).
x_grid (array_like) – Spatial grid(s). y_grid is required for 2D operators.
y_grid (array_like) – Spatial grid(s). y_grid is required for 2D operators.
kx (array_like) – Frequency grid(s) conjugate to x_grid/y_grid. ky is required for 2D operators.
ky (array_like) – Frequency grid(s) conjugate to x_grid/y_grid. ky is required for 2D operators.
boundary_condition (str, optional) – Boundary handling passed through to each sub-operator’s apply.
dealiasing_mask (optional) – Passed through to each sub-operator’s apply call.
freq_window (optional) – Passed through to each sub-operator’s apply call.
clamp (optional) – Passed through to each sub-operator’s apply call.
space_window (optional) – Passed through to each sub-operator’s apply call.
weyl_order (int, optional) – Truncation order used when converting a Weyl symbol to Kohn-Nirenberg.
apply_joint (bool, optional) – If False, skip the joint residual entirely (faster, but only an approximation) and warn that it was ignored.
decomposition (dict, optional) – A precomputed peetre_decomposition result to reuse instead of computing one. If omitted, it’s computed (and cached, per use_cache) internally.
use_cache (optional) – Forwarded to peetre_decomposition when decomposition is not supplied.
separable_local (optional) – Forwarded to peetre_decomposition when decomposition is not supplied.
joint_backend ({'direct', 'lowrank', 'nufft', 'aaa'}, optional) – Backend used to apply the joint residual.
joint_degree (optional) – Backend-specific parameters forwarded to _apply_joint_residual.
joint_tol (optional) – Backend-specific parameters forwarded to _apply_joint_residual.
joint_bounds (optional) – Backend-specific parameters forwarded to _apply_joint_residual.
joint_max_rel_error (optional) – Backend-specific parameters forwarded to _apply_joint_residual.
joint_num_samples (optional) – Backend-specific parameters forwarded to _apply_joint_residual.
joint_seed (optional) – Backend-specific parameters forwarded to _apply_joint_residual.

Returns:

Complex128 array, same shape as u, equal to the sum of the local, separable, and (if apply_joint) joint contributions.

Return type:

numpy.ndarray

Raises:

ValueError – If self.dim == 2 and y_grid or ky is not provided.

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 and Peetre decompositions.

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.
mode ({'kn', 'weyl'}, default='kn') – Quantization mode forwarded to compose_asymptotic for both A∘B and B∘A.
sign_convention ({'standard', 'inverse'}, optional) – Phase-factor convention forwarded to compose_asymptotic (controls the sign of [x, ξ]); see compose_asymptotic for details. If None, defaults to ‘standard’.

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, do_simplify=True)[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’.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

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.

decompose_symbol_peetre(*args, **kwargs)[source]: Alias for peetre_decomposition(), for compatibility with the standalone symbolic_decompose.py naming style.

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=2, mode='kn', sign_convention=None, do_simplify=True)[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=2) – Maximum order of the asymptotic expansion. Higher orders include more composition terms, improving accuracy for small t or when non-commutativity effects are significant.
mode ({'kn', 'weyl'}, default='kn') – Quantization mode forwarded to compose_asymptotic at each power P^n in the series.
sign_convention ({'standard', 'inverse'}, optional) – Phase-factor convention forwarded to compose_asymptotic; see compose_asymptotic for details. If None, defaults to ‘standard’.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

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.

fractional_power(alpha, order=1, method='symbolic', x_grid=None, L=None, N=None)[source]

Compute the symbol or discrete matrix representation of the fractional or complex power P^alpha of the pseudo-differential operator.

This method calculates the fractional power using either an optimized symbolic asymptotic approach or a numerical matrix approach.

For the symbolic method, the algorithm distinguishes between spatially independent symbols (pure Fourier multipliers) and spatially dependent symbols (heterogeneous media):

Fast Path (Pure Multipliers): If the symbol has no spatial dependence, it computes the exact algebraic power p(xi)^alpha. It automatically detects negative symbols (e.g., -xi**2) and extracts the appropriate complex phase factor to ensure the correct branch of the complex power is taken.

Slow Path (Spatially Dependent): If the symbol depends on space, the fractional power is currently restricted to square roots (alpha = 0.5). For order=0, it returns the principal square root. For order=1, it skips heavy iterative Newton-Raphson routines and directly computes the first-order microlocal spatial correction via the subprincipal error: delta_q = E / (2 * q_0).

For the numerical method, the operator is discretized into a matrix (via spectral or finite-difference methods) and the fractional matrix power is computed using SciPy.

Parameters:

alpha (float or complex) – The exponent to which the operator is raised.
order (int, default=1) – The asymptotic expansion order for the symbolic correction. - 0: Returns the principal symbol’s fractional power (q_0). - 1: Computes the first-order microlocal correction term (delta_q). Ignored for pure multipliers and the numerical method.
method (str, {'symbolic', 'numerical'}, default='symbolic') –
The computation method to use. - ‘symbolic’: Returns a SymPy expression for the fractional power

symbol using direct asymptotic calculus.
- ’numerical’: Returns a NumPy array representing the discrete fractional power matrix.
x_grid (ndarray, optional) – Spatial grid for the numerical method. If None, a default grid is generated based on L and N.
L (float, optional) – Domain half-length for the numerical spectral method.
N (int, optional) – Number of grid points for the numerical method.

Returns:

If method=’symbolic’, returns the SymPy expression of the fractional power symbol.
If method=’numerical’, returns the NumPy array representing the fractional power matrix.

Return type:

sympy.Expr or ndarray

Raises:

NotImplementedError – If method=’symbolic’ is used on a spatially dependent symbol with alpha != 0.5. The direct subprincipal asymptotic correction for heterogeneous media is currently optimized strictly for square roots to maintain clean, high-performance symbolic evaluations.
ValueError – If method is not ‘symbolic’ or ‘numerical’.

Notes

The symbolic method automatically handles negative overall symbols (e.g., the negative Laplacian) by factoring out the complex phase (e.g., I or I**(2*alpha)) and evaluating the absolute value.
For the numerical method, the underlying discrete operator is built using _build_operator_matrix, defaulting to a spectral method if x_grid, L, and N are not explicitly provided.
The fractional power of a pseudo-differential operator is a cornerstone in defining fractional PDEs (e.g., fractional Schrödinger equations, anomalous diffusion) and spectral zeta functions.

group_velocity_field(xlim=(-2, 2), klim=(-10, 10), density=30)[source]: Quiver plot of the group velocity field (1, dp/dxi). 1D only. cf. package microlocal.py.

interactive_symbol_analysis(xlim=(-2, 2), ylim=(-2, 2), xi_range=(0.1, 5), eta_range=(-5, 5), density=50)[source]: Launch an ipywidgets dashboard for symbol exploration. Same modes, same sliders, same defaults as before – rewritten as a mode-table dispatcher that delegates to the visualize_*/plot_* methods above instead of duplicating their 1D/2D branches inline (previously ~260 lines of near-duplicated if-elif chains). cf. package microlocal.py.

is_elliptic_numerically(x_grid, xi_grid, order=None, threshold=1e-06, xi_frac=0.5, n_random=100000, n_edge=32, seed=None)[source]

Check ellipticity using the standard asymptotic (high-frequency) definition.

A symbol p(x, ξ) of order m is elliptic if there exist constants C > 0 and R > 0 such that

|p(x, ξ)| ≥ C · (1 + |ξ|²)^(m/2) for all |ξ| ≥ R (and x in the
region of interest).

Ellipticity is a high-frequency condition: it makes no requirement on the symbol near ξ = 0, so this method only probes the region |ξ| ≥ R, where R is taken as a fraction (xi_frac) of the largest |ξ| present in xi_grid. This intentionally excludes the low-frequency region, so a symbol like the harmonic oscillator x² + ξ² — which has an isolated zero at (0, 0) but is the textbook example of a globally elliptic symbol — is correctly detected as elliptic instead of being penalised for vanishing at a single point that asymptotic ellipticity doesn’t care about.

Parameters:

x_grid (ndarray or tuple of ndarray) – Spatial grid(s). For 1D: a 1D array of x coordinates. For 2D: a tuple (x, y) of two 1D arrays.
xi_grid (ndarray or tuple of ndarray) – Frequency grid(s). For 1D: a 1D array of ξ coordinates. For 2D: a tuple (ξ, η) of two 1D arrays.
order (float, optional) – The homogeneity order m used in the (1 + |ξ|²)^(m/2) normalization. If None (default), it is auto-detected via self.symbol_order(). Raises ValueError if auto-detection fails and no order was supplied.
threshold (float, optional) – Minimum acceptable value of the normalized symbol |p(x,ξ)| / (1 + |ξ|²)^(m/2), i.e. the constant C in the ellipticity definition above. Default is 1e-6.
xi_frac (float, optional) – Fraction (in (0, 1)) of the largest |ξ| available in xi_grid used as the high-frequency cutoff R. Only points with |ξ| ≥ R are tested. Default is 0.5.
n_random (int, optional) – Number of random points sampled in the high-frequency region. Default 100_000.
n_edge (int, optional) – Number of spatial points per dimension used in the deterministic check performed at the largest available |ξ| (the edge of the supplied grid). Default 32.
seed (int, optional) – Seed for the random sampling, for reproducible results. Default None (non-deterministic across calls).

Returns:

True if the normalized symbol stays above threshold on every tested point in the high-frequency region, otherwise False.

Return type:

bool

Notes

This method combines two checks:

Random sampling in the shell |ξ| ≥ R (rejection-sampled from the hyperrectangle defined by xi_grid’s extremes).
Deterministic edge sampling — the symbol is evaluated at the largest available |ξ| (the corners/extremes of xi_grid) over a spatial grid of n_edge points per dimension, to guarantee the most informative high-frequency points are always checked, not just randomly sampled ones.

Because random sampling is probabilistic, a True result does not guarantee ellipticity in the strict mathematical sense — it only indicates that no tested point violated the condition. For rigorous analysis, combine this with symbolic methods (e.g., principal_symbol() and is_homogeneous()).

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.

kn_to_weyl_symbol(order=4)[source]

Convert the Kohn-Nirenberg symbol to its Weyl equivalent.

This is the inverse of weyl_to_kn_symbol:

a_Weyl(x, xi) = exp(+i/2 * d_x d_xi) a_KN(x, xi)
              ~ sum_{k=0}^{order}  (+i/2)^k / k!
                  * (d_x d_xi)^k a_KN(x, xi)

In 2D, d_x d_xi is replaced by d_x d_xi + d_y d_eta.

Parameters:: order (int, default=4) – Truncation order. Same remarks as for weyl_to_kn_symbol.
Returns:: Equivalent Weyl symbol (simplified SymPy expression).
Return type:: sympy.Expr

Examples

1D – multiplication-by-x*xi operator:

>>> op.kn_to_weyl_symbol(order=2)
x*xi + I/2      # opposite sign to weyl_to_kn_symbol

Notes

Typical use case: the operator is defined by its KN symbol (e.g. extracted in ‘auto’ mode), and you want the Weyl symbol to exploit its spectral properties – self-adjointness when a_Weyl is real, composition via the Moyal star product, etc.

Round-trip consistency check:

op_weyl = PseudoDifferentialOperator(op.kn_to_weyl_symbol(order=N), ...)
op_weyl.weyl_to_kn_symbol(order=N)  # should recover op.symbol up to order N

Submodules

psiop.matpsiop module

N x N matrix-valued pseudo-differential operators and coupled systems.

matpsiop — N × N matrix-valued pseudo-differential operators

Overview

The matpsiop submodule extends the scalar pseudo-differential framework to systems of coupled partial differential equations, Dirac-type equations, and matrix-valued fields (such as density matrices or matrix Green’s functions).

It provides the MatrixPseudoDifferentialOperator class, which wraps an N × N matrix of scalar symbols Pᵢⱼ(x, ξ) and orchestrates their numerical application and symbolic calculus. Because matrix multiplication is inherently non-commutative, this module implements specialized asymptotic expansions for composition, commutators, and exponentials that preserve matrix ordering.

Main object

MatrixPseudoDifferentialOperator: N × N matrix-valued pseudo-differential operator. Each entry Pᵢⱼ is internally wrapped as its own scalar PseudoDifferentialOperator. This design allows the matrix operator to seamlessly inherit all scalar numerical backends (FFT, Peetre decomposition, NUFFT, AAA, low-rank) entrywise, without requiring separate matrix-specific numerical kernels.

Key features

Vector and matrix field application:: Entrywise application to vector fields u = (u₁, …, uₙ); left and right matrix-multiplication actions on N × N matrix-valued fields (essential for Sylvester-type equations ∂ₜU = P U - U Q).
Symbol matrix evaluation and spectral analysis:: Pointwise numerical evaluation of the full (…, N, N) symbol on spatial/frequency grids; pointwise eigenvalues and eigenvectors of the symbol matrix (closed-form for N=2, used to build per-branch Hamiltonians for coupled systems).
Non-commutative asymptotic symbolic calculus:: Matrix composition P ∘ Q (exact for constant-coefficient entries, non-commutative even at 0th order); matrix commutators [P, Q]; formal left and right asymptotic inverses (requiring an invertible principal symbol); formal Hermitian adjoint P*; matrix exponential symbols exp(t · Op[P]) for propagators of coupled/vector-valued PDE systems.

Mathematical background

Matrix-valued symbols and quantization

A matrix pseudo-differential operator P acting on a vector field u(x) = (u₁(x), …, uₙ(x))ᵀ is defined by an N × N matrix of scalar symbols P(x, ξ). The Kohn–Nirenberg quantization is applied entrywise:

(P u)ᵢ(x) = (2π)⁻ᵈ Σⱼ ∫ exp(i x·ξ) Pᵢⱼ(x, ξ) ûⱼ(ξ) dξ

For matrix-valued fields U(x) (e.g., density matrices), the operator can act from the left (P U) or from the right (U Q), corresponding to the row-wise or column-wise application of the scalar entries.

Asymptotic matrix composition

The symbol of the composition P ∘ Q admits the Kohn–Nirenberg asymptotic expansion:

(P ∘ Q)ᵢₖ(x, ξ) ∼ Σⱼ Σ_α [ i^(-|α|) / α! ] ∂_ξ^α Pᵢⱼ(x, ξ) · ∂_x^α Qⱼₖ(x, ξ)

Unlike the scalar case, matrices do not commute. Therefore, the 0th-order term of the composition is the standard matrix product P(x, ξ) Q(x, ξ), and P ∘ Q ≠ Q ∘ P even for constant-coefficient symbols. For x-independent symbols, all higher-order derivative terms vanish identically, making the 0th-order matrix product the exact composition symbol.

Matrix exponential and propagators

The symbol of the exponential operator exp(t P) is computed via a truncated Taylor series using the matrix composition rule:

exp(t P) ∼ I + t P + (t²/2!) (P ∘ P) + …

This is the fundamental building block for time-stepping coupled systems ∂ₜu = P u. Because the composition is exact for constant-coefficient matrices, the propagator symbol for such systems reduces to the standard matrix exponential exp(t P(ξ)).

Formal adjoints and inverses

The formal Hermitian adjoint P* is obtained by taking the formal adjoint of each scalar entry (complex conjugation + asymptotic expansion in |ξ| → ∞), and then transposing the resulting matrix.

Formal asymptotic inverses R (such that P ∘ R ∼ I) require the principal symbol matrix to be invertible (det P ≠ 0). The recursion preserves matrix multiplication order, ensuring the inverse factor remains on the correct side to cancel P.

Numerical design notes

Entrywise numerical orchestration:: MatrixPseudoDifferentialOperator introduces no new numerical kernels for application. The apply() method simply orchestrates N² independent calls to the scalar PseudoDifferentialOperator.apply(). Consequently, it automatically inherits the Peetre decomposition, NUFFT, AAA, and low-rank backends for each entry.
Memory and performance:: Because each of the N² scalar operators is constructed and cached independently, memory usage scales linearly with N². For large systems, constant-coefficient entries benefit from the scalar fast-path FFT multiplier, making matrix application highly efficient.
Sylvester-type equations:: The apply_matrix_field (left action) and apply_matrix_field_right (right action) methods are designed to support operator splitting for equations of the form ∂ₜU = P U - U Q. Left and right actions commute as operations ((P U) Q = P (U Q)), enabling efficient Lie-Trotter or Strang splitting schemes even when the underlying scalar operators do not commute.

class psiop.matpsiop.MatrixPseudoDifferentialOperator(P_expr, vars_x, mode='symbol', quantization='kohn-nirenberg', apply_backend='peetre', compute_peetre=False, peetre_options=None)[source]

Bases: object

N x N matrix-valued pseudodifferential operator, built from a sympy Matrix of scalar symbols P_ij(x[, y], xi[, eta]).

Each entry P_ij is wrapped as its own scalar PseudoDifferentialOperator (same expr/mode/quantization/apply_backend conventions), so apply() reuses the existing periodic FFT / Peetre machinery entirely. Entries may depend on x (and y) for variable-coefficient systems, or be xi`(, eta`)-only for constant-coefficient ones – both are supported by the same class; nothing here forces one or the other.

Parameters:

P_expr (sympy.Matrix or nested list of sympy.Expr) – N x N matrix of scalar symbol expressions, in the same x[, y], xi[, eta] convention as PseudoDifferentialOperator.
vars_x (list of sympy symbols) – Spatial variables; length 1 or 2, as for PseudoDifferentialOperator.
mode – Forwarded to each entry’s PseudoDifferentialOperator.
quantization – Forwarded to each entry’s PseudoDifferentialOperator.
apply_backend – Forwarded to each entry’s PseudoDifferentialOperator.
compute_peetre – Forwarded to each entry’s PseudoDifferentialOperator.
peetre_options – Forwarded to each entry’s PseudoDifferentialOperator.

size

Matrix dimension N (the constructor accepts any square N, but eigen_symbol’s closed-form path is specific to N=2).

Type:: int

entries

entries[i][j] is the scalar operator for P_ij.

Type:: list of list of PseudoDifferentialOperator

apply(u, x_grid, kx, y_grid=None, ky=None, **apply_kwargs)[source]

Apply P(x, xi) to a vector field u = (u_1, …, u_N).

Parameters:

u (sequence of N ndarrays) – Vector field components, each sampled on the grid.
x_grid – As for PseudoDifferentialOperator.apply.
kx – As for PseudoDifferentialOperator.apply.
y_grid – As for PseudoDifferentialOperator.apply.
ky – As for PseudoDifferentialOperator.apply.
**apply_kwargs – Forwarded to every entry’s scalar apply() (e.g. boundary_condition, freq_window, …).

Returns:

(P u)_i = sum_j Op[P_ij](u_j).

Return type:

list of N ndarrays

apply_matrix_field(U, x_grid, kx, y_grid=None, ky=None, **apply_kwargs)[source]

Apply P(x, xi) to a matrix-valued field U(x) by left matrix multiplication on U’s own N x N structure:

(P U)_ik(x) = sum_j Op[P_ij] (U_jk) (x)

Unlike apply(), which propagates a single vector field u = (u_1, …, u_N), here U itself carries an extra N x N index pair (e.g. a density matrix or matrix Green’s function) that P acts on only from the left. Each column U[:, k] of U is an ordinary vector field, so this reduces to N independent calls to apply(), one per column, with the results reassembled into the matrix-shaped output; P is applied to U, never U to P.

Parameters:

U (sequence of N sequences of N ndarrays, or ndarray of shape (N, N, ...)) – Matrix-valued field; U[j][k] (equivalently U[j, k] for an ndarray) is the scalar (j, k) component field sampled on the grid, so that U plays the role of an N x N matrix at every grid point.
x_grid – As for apply.
kx – As for apply.
y_grid – As for apply.
ky – As for apply.
**apply_kwargs – Forwarded to every entry’s scalar apply() (e.g. boundary_condition, freq_window, …).

Returns:

out[i][k] holds (P U)_ik, indexed the same way as U.

Return type:

list of N lists of N ndarrays

Raises:

ValueError – If U is not an N x N array of fields, with N equal to self.size.

apply_matrix_field_right(U, x_grid, kx, y_grid=None, ky=None, **apply_kwargs)[source]

Apply this operator’s symbol Q(x, xi) to a matrix-valued field U(x) by right matrix multiplication on U’s own N x N structure:

(U Q)_ik(x) = sum_j Op[Q_jk] (U_ij) (x)

This is the mirror image of apply_matrix_field (which left-multiplies by P): here each row U[i, :] of U is an ordinary vector field acted on from the right by Q, which is equivalent to the left action of the transposed symbol matrix Q^T on that row – hence the index order Op[Q_jk], not Op[Q_kj], so the two methods are genuinely different unless Q is symmetric.

Together with apply_matrix_field, this is the numerical primitive needed to time-step Sylvester-type equations d_t U = P U - U Q, since left- and right-multiplication always commute as operations ((P U) Q == P (U Q)), even though the underlying scalar operators Op[P_ij] and Op[Q_jk] need not commute with each other when the symbols depend on x. See solve_sylvester_field for the corresponding time-stepper.

Parameters:

U (sequence of N sequences of N ndarrays, or ndarray of shape (N, N, ...)) – Matrix-valued field; U[i][j] (equivalently U[i, j] for an ndarray) is the scalar (i, j) component field sampled on the grid.
x_grid – As for apply.
kx – As for apply.
y_grid – As for apply.
ky – As for apply.
**apply_kwargs – Forwarded to every entry’s scalar apply().

Returns:

out[i][k] holds (U Q)_ik, indexed the same way as U.

Return type:

list of N lists of N ndarrays

Raises:

ValueError – If U is not an N x N array of fields, with N equal to self.size.

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

Symbol of the commutator [Op[self], Op[other]], generalizing PseudoDifferentialOperator.commutator_symbolic to matrices.

Unlike the scalar case (whose 0th-order term always vanishes, since scalars commute), the matrix commutator is generally nonzero already at 0th order: it’s the ordinary matrix commutator P(x,xi) Q(x,xi) - Q(x,xi) P(x,xi). Higher orders add the noncommutative analogue of the Poisson-bracket correction. For constant-coefficient self/other this is exact, and reduces exactly to the plain matrix commutator (see compose_asymptotic).

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

Compose two matrix-valued symbols via the same asymptotic Kohn-Nirenberg / Weyl expansion as PseudoDifferentialOperator.compose_asymptotic, generalized to (order-preserving) matrix multiplication: this is the symbol of Op[self] . Op[other], and self.compose_asymptotic(other) != other.compose_asymptotic(self) in general – unlike the scalar case, matrices don’t commute even at 0th order.

For x[, y]-independent (“constant-coefficient”) symbols this is exact at any order: every n>=1 term involves a spatial derivative of a xi-only expression, which vanishes identically, so the result reduces to the ordinary matrix product P(xi) Q(xi).

Parameters:

other (MatrixPseudoDifferentialOperator) – Same size and dim as self.
order – As for the scalar compose_asymptotic.
mode – As for the scalar compose_asymptotic.
do_simplify (sign_convention and) – As for the scalar compose_asymptotic.

Returns:

The composed symbol.

Return type:

sympy.Matrix, shape (size, size)

eigen_symbol(*args)[source]: Pointwise eigenvalues/eigenvectors of the symbol matrix P(x[, y], xi[, eta]) at given point(s). If called without arguments, computes symbolic eigenvalues/eigenvectors.

exponential_symbol(t=1.0, order=2, mode='kn', sign_convention=None, do_simplify=True)[source]

Symbol of exp(t Op[self]) for the matrix-valued operator, via the matrix analogue of PseudoDifferentialOperator.exponential_symbol.

Same truncated power series as the scalar case,

exp(tP) ~ I + t P + (t^2/2!) P^{.2} + (t^3/3!) P^{.3} + …

but “P^{.n}” means the symbol of Op[P] . Op[P] . … . Op[P] (n times), computed via the matrix compose_asymptotic – i.e. ordinary matrix multiplication P @ P order-corrected by the KN/Weyl derivative terms – since matrix symbols do not commute and Matrix.__mul__(P, P) alone is only the 0th-order (frozen-coefficient) approximation to that composition. Works for both dim == 1 and dim == 2 – compose_asymptotic already branches on dimension internally.

Parameters:

t (float or sympy.Symbol, default=1.0) – Evolution parameter, same conventions as the scalar version (e.g. t = -i*tau for exp(-i*tau*H), t = tau for exp(tau*Delta)).
order (int, default=2) – Truncation order, used both for the outer Taylor series and as the order passed to each compose_asymptotic call.
mode ({'kn', 'weyl'}, default='kn') – Quantization convention for the composition (2D Weyl is not implemented for matrix symbols – see compose_asymptotic).
sign_convention (optional) – Forwarded to compose_asymptotic.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Returns:

Truncated symbol of exp(t Op[self]).

Return type:

sympy.Matrix, shape (size, size)

Notes

For x[, y]-independent (“constant-coefficient”) self, compose_asymptotic is exact, so this reduces to the exact truncated matrix exponential series of P(xi[, eta]); compare against scipy.linalg.expm at sample points to sanity-check.
Non-commutativity means self and other’s roles in each compose_asymptotic call matter; here every factor is self, so ordering is moot, but see commutator_symbolic for the general two-operator case.

formal_adjoint(n_terms=6)[source]

Formal Hermitian adjoint symbol P* of the matrix operator.

Each entry gets the same scalar treatment as PseudoDifferentialOperator.formal_adjoint (conjugate + asymptotic expansion at infinity in |xi|); the resulting matrix is then transposed (not conjugate-transposed again – conjugation already happened entrywise) because (Op[P]u, v) = (u, Op[P]* v) swaps the row/column roles of the symbol, same as for a plain matrix adjoint.

left_inverse_asymptotic(order=1)[source]: Formal left inverse L such that Op[L] . Op[self] ~ Id up to O(<xi>^-order), matrix analogue of PseudoDifferentialOperator.left_inverse_asymptotic. Requires the symbol P(x, xi) to be invertible as a matrix.

right_inverse_asymptotic(order=1)[source]: Formal right inverse R such that Op[self] . Op[R] ~ Id up to O(<xi>^-order), matrix analogue of PseudoDifferentialOperator.right_inverse_asymptotic. Requires the symbol P(x, xi) to be invertible as a matrix.

symbol_matrix(*args)[source]

Numerically evaluate P(x[, y], xi[, eta]) at a point or broadcastable arrays, returning an ndarray of shape (…, N, N). If called without arguments, returns the symbolic sympy.Matrix.

Parameters:: *args – The point(s) to evaluate at, in the order each entry’s p_func expects: (x, xi) for 1D, (x, y, xi, eta) for 2D. Arguments may be broadcastable ndarrays (e.g. full grids), in which case the leading dimensions of the output match their broadcast shape.
Return type:: ndarray, shape (…, N, N) or sympy.Matrix if no args are provided.

psiop.psiop_apply module

Standalone numerical kernels for Kohn-Nirenberg quantization, NUFFT, AAA rational approximation, and Low-Rank SVD factorization.

psiop_apply.py — Numerical backends for pseudo‑differential operator application

Overview

The psiop_apply module provides the heavy‑lifting numerical kernels for applying pseudo‑differential operators (ΨDOs) to spatial fields. While the core psiop module handles symbolic calculus, asymptotic expansions, and the Peetre decomposition, this module is responsible for the actual evaluation of the resulting integrals on discrete spatial and frequency grids.

It implements both periodic (FFT‑based) and non‑periodic (direct quadrature) Kohn–Nirenberg quantization, alongside specialized factorization and approximation backends for the genuinely joint space‑frequency residuals that arise in the Peetre decomposition.

Main objects and workflows

Periodic and Non‑Periodic Kohn–Nirenberg Application

kohn_nirenberg_fft: Applies the operator on a periodic domain using FFTs. Features an automatic fast‑path for spatially independent symbols (pure Fourier multipliers) and a memory‑bounded, multi‑threaded slow path for spatially dependent symbols.

kohn_nirenberg_nonperiodic: Applies the operator on a non‑periodic (Dirichlet‑type) domain via direct discrete Fourier quadrature. Utilizes aggressive caching of phase matrices and windowing arrays to minimize redundant computations across repeated calls.

Low‑Rank Chebyshev/SVD Factorization

factorize_symbolic: Approximates a joint space‑frequency symbol p(x, ξ) as a sum of separable terms ∑ₖ aₖ(x) qₖ(ξ) using Chebyshev interpolation followed by Singular Value Decomposition (SVD) truncation. Includes Monte‑Carlo quality diagnostics.

NUFFT‑Based Joint‑Residual Backend

try_nufft_decomposition_* / apply_nufft_*: Targets joint residuals with genuinely oscillatory phases of the form exp(i · Λ(x) · M(ξ)). Extracts the phase and amplitude, and evaluates the resulting non‑uniform FFT (Type 3) via the optional finufft library, with a pure‑NumPy O(N·M) direct‑sum fallback.

AAA Rational Approximation Backend

try_aaa_decomposition_* / aaa_plan_to_callable_*: Targets joint residuals that are rational functions or exhibit explicit poles / algebraic decay. Uses a vector‑valued Adaptive Antoulas‑Algorithm (AAA) barycentric rational interpolation to build a compact, fast‑evaluating surrogate for the symbol.

Key features

Memory‑bounded execution:: The 1D and 2D slow paths avoid O(N²ᵈ) RAM allocation by evaluating the symbol on chunked space‑frequency sub‑grids (~256 MB max per block) and accumulating via optimized Einstein summation (np.einsum).
Multi‑threaded row‑blocking:: The 2D slow path distributes spatial row‑blocks across a ThreadPoolExecutor, achieving near‑linear speedup on multi‑core machines for spatially dependent symbols.
Automatic fast‑path detection:: Before executing the expensive quadrature, the periodic and non‑periodic kernels probe the symbol at a few test frequencies. If the symbol is spatially independent, the code bypasses the quadrature entirely and applies the symbol as a pure Fourier multiplier, reducing complexity from O(N²) to O(N log N) in 1D, and O(N⁴) to O(N² log N) in 2D.
Phase‑matrix caching:: Non‑periodic transforms precalculate and cache discrete Fourier transform phases, reconstruction phases, and window arrays. The cache keys are derived from grid shapes and endpoints, ensuring automatic invalidation when grid resolution changes.
Quality‑gated approximations:: The low‑rank, NUFFT, and AAA backends all compute relative L2 errors against the exact symbol. If the approximation error exceeds the requested tolerance, the backend gracefully falls back to the exact (but slower) direct Kohn–Nirenberg quadrature.

Mathematical background and numerical design

Kohn–Nirenberg quantization (Periodic)

The operator Op(p) is applied to a periodic function u via:

[Op(p) u](x) = (2π)⁻ᵈ ∫ p(x, ξ) e^{i x·ξ} ℱ[u](ξ) dξ

where ℱ[u] is the discrete Fourier transform. Fast‑path: If p(x, ξ) = p(ξ), the x‑dependence drops out, and the integral collapses to the pure multiplier:

[Op(p) u](x) = ℱ⁻¹[ p(ξ) · ℱ[u](ξ) ]

Slow‑path: For space‑dependent p(x, ξ), the integral is evaluated directly. To prevent memory exhaustion, the spatial domain is split into blocks of size B, and the frequency domain into chunks of size C. The quadrature is accumulated block‑by‑block:

result[i₀:i₁] = (Δξ / 2π) ∑_{k‑chunk} P_{blk} · ℱ[u]_{chunk} · e^{i x_{blk} · ξ_{chunk}}

Kohn–Nirenberg quantization (Non‑Periodic)

On a non‑periodic domain, the continuous Fourier transform is replaced by a direct discrete quadrature:

[Op(p) u](x) = (2π)⁻ᵈ ∫ p(x, ξ) e^{i x·ξ} [ ∫ e^{-i y·ξ} u(y) dy ] dξ

The inner integral (forward transform) and outer integral (reconstruction) are represented as dense matrix‑vector products using precomputed phase matrices Φ_{ft} = e^{-i ξ xᵀ} and Φ_{rec} = e^{i x ξᵀ}. These matrices are cached globally. The fast‑path logic is identical to the periodic case, bypassing the matrix multiplications when p is independent of x.

Peetre Joint Residual Factorization

When the Peetre decomposition yields a genuinely joint residual p_joint(x, ξ) that cannot be written as a(x)q(ξ), it is routed to one of three specialized backends based on its algebraic structure:

Low‑Rank (Chebyshev/SVD): For smooth, non‑oscillatory kernels. The symbol is interpolated on a tensor‑product Chebyshev grid, reshaped into a matrix C ∈ ℂ^{N_x × N_ξ}, and truncated via SVD:

C ≈ U_r Σ_r V_r^H ⇒ p_joint(x, ξ) ≈ ∑_{k=1}^r aₖ(x) qₖ(ξ)

The basis functions are explicit Chebyshev polynomials mapped to the physical bounding box.
NUFFT (Oscillatory): For residuals containing a bilinear phase exp(i Λ(x) M(ξ)). The symbol is factored as:

p_joint(x, ξ) = c(x) g(ξ) exp(i Λ(x) M(ξ))

The application is reformulated as a Type 3 Non‑Uniform FFT, evaluating the sum:

f(x) = ∑_{j} w_j exp(i (x · Λ(x) + μ_j · M(ξ_j)))

where w_j are the weighted Fourier coefficients of u. This achieves O(N log N) complexity via finufft, avoiding the polynomial basis convergence issues of the low‑rank method.
AAA (Rational / Poles): For residuals with explicit poles or algebraic decay (e.g., resolvent‑like structures). A vector‑valued AAA barycentric rational interpolant is constructed:

r(ξ) = ∑_{k} wₖ fₖ / (ξ - zₖ) / ∑_{k} wₖ / (ξ - zₖ)

where the support points zₖ and weights wₖ are selected adaptively to minimize the residual. The spatial dependence is handled by building a separate AAA fit for each Chebyshev node in x, followed by barycentric Lagrange interpolation in x.

Numerical stability

All application kernels enforce numerical stability through:

Magnitude clamping: Symbol values exceeding clamp (default 10⁶) are scaled down while preserving their complex phase.
Frequency windowing: Optional Gaussian or Hann tapers in the frequency domain to attenuate high‑frequency numerical artifacts.
Spatial tapering: Optional centered Gaussian tapers in the spatial domain to mitigate edge boundary artifacts in non‑periodic settings.

psiop.psiop_apply.aaa_plan_to_callable_1d(plan)[source]: Wrap an aaa_decomposition_1d plan as p(x, xi) -> ndarray, matching the symbol_func signature kohn_nirenberg_fft/nonperiodic expect.

psiop.psiop_apply.aaa_plan_to_callable_2d(plan)[source]: Wrap an aaa_decomposition_2d plan as p(x, y, xi, eta) -> ndarray, matching the symbol_func signature kohn_nirenberg_fft/nonperiodic expect for dim=2.

psiop.psiop_apply.apply_nufft_1d(u, plan, x_grid, kx, dx, dxi, eps=1e-12, freq_window='gaussian')[source]

Apply Op(p_joint) via the NUFFT tier, 1D case. plan is the output of try_nufft_decomposition_1d (a list of term dicts).

freq_window matches kohn_nirenberg_fft’s default – without applying it here too, results silently diverge from joint_backend=’direct’ even at freq_window=’gaussian’ defaults (found via end-to-end testing against the real dispatcher, not from the isolated unit tests, which never exercised the default windowing at all).

psiop.psiop_apply.apply_nufft_2d(u, kind, plan, x_grid, y_grid, kx, ky, dx, dy, dxi, deta, eps=1e-12, freq_window='gaussian')[source]: Apply Op(p_joint) via the NUFFT tier, 2D case. kind/plan are the output of try_nufft_decomposition_2d. See apply_nufft_1d docstring on why freq_window must be matched to the direct path’s default.

psiop.psiop_apply.evaluate_decomposition_quality(orig_expr, symbolic_pairs, x_syms, xi_syms, bounds, num_samples=10000, seed=42)[source]

Estimate the symbol-level approximation error of a separable/low-rank decomposition against the original expression, via Monte Carlo sampling at random off-grid points (so the error reflects genuine approximation quality rather than exact agreement at the fitting nodes):

orig_expr(x, xi) ≈ sum_k a_k(x) q_k(xi)

Parameters:

orig_expr (sympy.Expr) – Original joint symbol being approximated.
symbolic_pairs (list of tuple) – Candidate decomposition, as pairs (a_k(x), q_k(xi)) of sympy expressions.
x_syms (list of sympy symbols) – Spatial variables of orig_expr.
xi_syms (list of sympy symbols) – Frequency variables of orig_expr.
bounds (dict) – Mapping from each symbol in x_syms + xi_syms to a (min, max) sampling range.
num_samples (int, default 10000) – Number of random points drawn uniformly within bounds.
seed (int, default 42) – Seed for the random number generator, for reproducibility.

Returns:

Dictionary with keys ‘rel_l2_error’ (relative L2 error over the sampled points), ‘max_abs_error’ and ‘mean_abs_error’ (pointwise absolute-error statistics).

Return type:

dict

psiop.psiop_apply.factorize_symbolic(expr, x_syms, xi_syms, bounds, degree=6, tol=1e-05, num_samples=10000, seed=42, digits=5)[source]

Low-rank Chebyshev/SVD factorization of a joint symbol:

p(x, xi) ≈ sum_{k=1}^r a_k(x) q_k(xi)

The approximation is valid on the bounded rectangle given by bounds.

Parameters:

expr (sympy.Expr) – Symbol to factorize, usually the Peetre joint residual.
x_syms (list of sympy symbols) – Spatial variables.
xi_syms (list of sympy symbols) – Frequency variables.
bounds (dict) – Dictionary mapping each symbol to (min, max).
degree (int) – Chebyshev degree in each variable.
tol (float) – Relative singular-value cutoff and coefficient pruning threshold.
num_samples (int) – Number of Monte Carlo samples for quality diagnostics.
seed (int) – RNG seed.
digits (int) – Number of digits used when converting floating coefficients to SymPy.

Returns:

symbolic_pairs (list of tuple) – List of (a_k(x), q_k(xi)) SymPy expressions.
metrics (dict) – Symbol-level approximation diagnostics.

psiop.psiop_apply.invalidate_kn_cache() → None[source]: Clear the phase-matrix cache for non-periodic 1D operations.

psiop.psiop_apply.kohn_nirenberg_fft(u_vals: ndarray, symbol_func: Callable[[...], ndarray], x_grid: ndarray, kx: ndarray, fft_func: Callable, ifft_func: Callable, dim: int = 1, y_grid: ndarray | None = None, ky: ndarray | None = None, freq_window: str | None = 'gaussian', clamp: float = 1000000.0, space_window: bool = False, is_spatial: bool = False) → ndarray[source]

Numerically stable Kohn–Nirenberg quantization of a pseudo-differential operator on a periodic domain using FFTs.

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

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

where ℱ[u] is the discrete Fourier transform of u, and p(x, ξ) is a symbol that may depend on both spatial variables (x, y) and frequency variables (ξ, η).

This implementation supports 1D and 2D spatial dimensions, providing numerical stability through symbol magnitude clamping, optional frequency windowing (Gaussian/Hann), and optional spatial tapering.

Fast-Path Optimization (Spatial Independence)

When the symbol p is independent of spatial coordinates (and space_window=False), the operator simplifies to a pure Fourier multiplier: ifft(P * fft(u)).

The function detects spatial independence using a multi-point evaluation heuristic across non-zero frequencies (to avoid false positives at ξ = 0). When active, computational complexity drops from O(N^2) to O(N log N) in 1D, and from O(N^4) to O(N^2 log N) in 2D.

Memory-Bounded Slow Path (Spatial Dependence)

For space-dependent symbols, the implementation avoids O(N^2d) RAM allocation: - 1D: Slices spatial evaluation into memory-bounded chunks (~256 MB max). - 2D: Combines parallel multi-threaded row-blocking with 2D frequency chunking

and phased factorized inner products (np.einsum).

Parameters:

u_vals (ndarray) – Spatial samples of the input field u(x) or u(x, y).
symbol_func (callable) – Symbol evaluator p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. Must accept NumPy-broadcastable positional arguments. Return arrays are automatically broadcasted and cast to complex128 to handle scalar outputs safely (e.g., from sympy.lambdify).
x_grid (ndarray) – 1D array of spatial coordinates along the x-axis.
kx (ndarray) – 1D array of spatial frequencies in the x-direction.
fft_func (callable) – Forward Fourier transform function (e.g., scipy.fft.fft or fft2).
ifft_func (callable) – Inverse Fourier transform function (e.g., scipy.fft.ifft or ifft2).
dim ({1, 2}, default=1) – Spatial dimensionality of the domain.
y_grid (ndarray, optional) – 1D array of spatial coordinates along the y-axis (required if dim=2).
ky (ndarray, optional) – 1D array of spatial frequencies in the y-direction (required if dim=2).
freq_window ({'gaussian', 'hann', None}, default='gaussian') – Frequency-domain window/taper applied to attenuate high-frequency instabilities.
clamp (float, default=1e6) – Maximum allowed magnitude for the symbol entries. Magnitudes exceeding this value are clipped to prevent overflow.
space_window (bool, default=False) – If True, applies a centered Gaussian spatial taper to mitigate edge boundary artifacts. Note: Setting `space_window=True` forces execution through the slow path.
is_spatial (bool or None, default=None) – Explicit hint about whether the 2D symbol depends on (x1, x2). True forces the slow (space-dependent) path, False forces the fast (space-independent) path, None triggers the sampling heuristic.

Returns:

Resulting array of the same shape and type (complex128) as u_vals after applying the pseudo-differential operator.

Return type:

ndarray

Raises:

ValueError – If dim=2 and y_grid or ky are not supplied, or if dim is not 1 or 2.

psiop.psiop_apply.kohn_nirenberg_nonperiodic(u_vals: ndarray, x_grid: ndarray | Tuple[ndarray, ndarray], xi_grid: ndarray | Tuple[ndarray, ndarray], symbol_func: Callable[[...], ndarray], freq_window: str | None = 'gaussian', clamp: float = 1000000.0, space_window: bool = False, is_spatial: bool | None = None, _cache: Dict = {}) → ndarray[source]

Numerically stable Kohn–Nirenberg quantization of a pseudo-differential operator on a non-periodic domain using direct matrix/quadrature-based transforms.

Applies the pseudo-differential operator Op(p) to the input function u via the non-periodic Kohn–Nirenberg integral formula:

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

where ℱ_NP[u] is the direct discrete Fourier integral transform evaluated over an arbitrary non-periodic spatial grid x and frequency grid xi.

Supports 1D and 2D spatial dimensions, featuring a caching mechanism for both 1D and 2D phase/window matrices, symbol magnitude clamping, optional frequency windowing (Gaussian/Hann), and spatial tapering.

1D Cached Path

In 1D, precalculates and caches discrete Fourier transform phases (phase_ft), reconstruction phases (exp_matrix), and window arrays in _cache to accelerate repeated function evaluations on identical grids.

2D Cached Path

In 2D, precalculates and caches the analogous grid-only objects — forward phase matrices (phase1, phase2), reconstruction phase matrices (exp1, exp2), the frequency window (freq_win_2d), and (if space_window=True) the spatial taper arrays (sw_x1_full, sw_x2). The cache key includes freq_window and space_window since, unlike the 1D case, these options change the cached content itself, not just which arrays get used.

2D Fast-Path Optimization (Spatial Independence)

Symbol spatial-dependence is controlled by is_spatial:

is_spatial=False forces the fast path (symbol treated as x-independent).
is_spatial=True forces the slow path (symbol treated as x-dependent).
is_spatial=None (default) runs a multi-point sampling heuristic (skipped, and slow path forced, whenever space_window=True, since the taper itself introduces x-dependence).

When spatial independence holds, matrix multiplications perform global frequency filtering in O(N^3) complexity instead of the full O(N^4) space-dependent integration.

2D Memory-Bounded Slow Path (Spatial Dependence)

For space-dependent symbols in 2D, a multi-tiered execution strategy prevents RAM spikes: - Row-Based Parallelization: Slices x1 into spatial blocks distributed across

a worker thread pool (ThreadPoolExecutor).

Dual Frequency Chunking: Iterates across sub-blocks of xi1 and xi2 to guarantee intermediate tensor evaluation (sv_chunk, phase_chunk) remains strictly bounded (~256 MB max).
Tensor Contraction: Employs optimized Einstein summation (np.einsum) for localized quadrature integration.

Parameters:

u_vals (ndarray) – Spatial samples of the input field u(x) [1D] or u(x1, x2) [2D].
x_grid (ndarray or tuple of ndarray) – Spatial coordinate grid x (1D) or tuple (x1, x2) (2D).
xi_grid (ndarray or tuple of ndarray) – Frequency grid xi (1D) or tuple (xi1, xi2) (2D).
symbol_func (callable) – Symbol evaluator p(x, ξ) in 1D or p(x1, x2, ξ1, ξ2) in 2D. Must accept NumPy-broadcastable positional arguments. Returns are automatically broadcasted, type-cast to complex128, and reshaped/copied safely to handle scalar or reduced-dimension outputs (e.g., from sympy.lambdify).
freq_window ({'gaussian', 'hann', None}, default='gaussian') – Frequency-domain window/taper applied to attenuate high-frequency numerical artifacts.
clamp (float, default=1e6) – Maximum allowed magnitude for symbol values. Entries exceeding this threshold are clipped to prevent overflow.
space_window (bool, default=False) – If True, applies a centered Gaussian spatial taper to attenuate edge boundary artifacts. Note: Enabling this disables the 2D spatial-independence fast path.
is_spatial (bool or None, default=None) – Explicit hint about whether the 2D symbol depends on (x1, x2). True forces the slow (space-dependent) path, False forces the fast (space-independent) path, None triggers the sampling heuristic.
_cache (dict, optional) – Global or local cache dictionary storing reusable 1D and 2D phase/window matrices. Defaults to module-level _KN_CACHE.

Returns:

Resulting complex-valued array (complex128) of the same dimensionality and shape as u_vals after applying the non-periodic pseudo-differential operator.

Return type:

ndarray

Raises:

NotImplementedError – If u_vals.ndim is not 1 or 2.

psiop.psiop_apply.try_aaa_decomposition_1d(joint_expr, x_sym, xi_sym, x_bounds, xi_bounds, n_cheb=24, n_xi_samples=100, rtol=1e-08)[source]: 1D bivariate rational decomposition via vector-AAA. Returns a plan dict (with a fast numpy callable, see aaa_plan_to_callable_1d) or None if the quality gate (rel_l2_error > 10*rtol) isn’t met.

psiop.psiop_apply.try_aaa_decomposition_2d(joint_expr, x_sym, y_sym, xi_sym, eta_sym, x_bounds, y_bounds, xi_bounds, eta_bounds, n_cheb_x=10, n_cheb_y=10, n_xi_samples=30, n_eta_samples=30, rtol=1e-08)[source]: 2D decomposition via sequential vector-AAA (xi support points chosen at a representative eta slice – see module docstring caveat above; stage 2 compresses eta from the EXACT symbolic slice at each xi support point). Returns a plan dict or None if the quality gate fails.

psiop.psiop_apply.try_nufft_decomposition_1d(joint_expr, x_sym, xi_sym)[source]: 1D (phase space (x,xi)) NUFFT classifier. Returns a list of term plans, or None if any additive term doesn’t fit (falls back).

psiop.psiop_apply.try_nufft_decomposition_2d(joint_expr, x_sym, y_sym, xi_sym, eta_sym)[source]

2D (phase space (x,y,xi,eta)) NUFFT classifier. Tries, in order:

axis-separable: term factors as A(x,xi)*B(y,eta) (disjoint variable groups) – cheapest, two independent 1D passes.
single-joint-term: term’s phase is one product Lambda(x,y)*M(xi,eta) – needs a 3D NUFFT embedding.

A symbol whose terms need genuinely independent coupling on BOTH axes simultaneously (4D embedding) is not representable by either tier and returns None (finufft has no type-3 transform above 3D). Returns (‘axis_sep’, combo_plan) or (‘joint3d’, plans) or None.

psiop.psiop_solver module

Time-stepping solvers, exponential propagators, and grid generation utilities.

psiop_solver — Time-stepping solvers, propagators, and grid utilities

Overview

The psiop_solver submodule orchestrates the numerical application of scalar and matrix-valued pseudo-differential operators (from psiop and psiop.matpsiop) to evolve fields in time. It provides high-level time-stepping routines, amortized exponential propagators, operator-splitting schemes for coupled matrix systems, and specialized quasi-linear solvers, along with the underlying periodic spatial and frequency grid generation.

The submodule is designed to seamlessly bridge the symbolic asymptotic calculus of psiop with efficient, memory-bounded numerical integration, handling both simple linear evolution and complex non-commutative or quasi-linear systems.

Main objects and workflows

PropagatorFamily and build_propagator: Construct the one-step approximate propagator exp(dt · Op[s]) for a given symbol s. PropagatorFamily amortizes the cost by building the asymptotic exponential-symbol expansion once with dt left as a free symbolic parameter, allowing subsequent propagators for different step sizes to be generated via a single, cheap sympy.subs call rather than repeating the full compose_asymptotic recursion.
solve_first_order: Time-steps the first-order evolution equation ∂ₜu = Op(s)(u) by repeated application of the asymptotic propagator.
solve_second_order: Handles second-order-in-time equations ∂²ₜu = Op(S)(u) by reducing them to a first-order block companion system of double the dimension, then time-stepping via solve_first_order.
solve_matrix_field: Extends the propagator machinery to matrix-valued fields U(x), evolving ∂ₜU = Op(P)U where P acts on U from the left (e.g., density matrices or matrix Green’s functions).
solve_sylvester_field: Time-steps Sylvester-type equations ∂ₜU = Op(P)U − U·Op(Q) using Lie-Trotter or Strang operator splitting between the independent left-acting P and right-acting Q symbols.
solve_ricci_flow_conformal_2d: Integrates the 2D conformal Ricci flow ∂ₜφ = e⁻²ᵠΔφ. Because the coefficient depends on the evolving solution, this quasi-linear equation cannot be described by a fixed symbol. It uses an IMEX/Lie splitting: an explicit Euler correction for the deviation from the spatial average, followed by an exact stiff step for the averaged part applied directly via FFT.
Grid utilities: make_grid_1d and make_grid_2d construct uniform periodic spatial grids and their associated FFT-ordered angular frequency grids.

Key features

Amortized propagator construction:: PropagatorFamily computes the symbolic exp(t · Op[s]) once; build_propagator caches these families process-wide, making parameter sweeps, adaptive stepping, and dt-convergence studies extremely cheap.
Block companion reduction:: Automatic conversion of scalar or matrix second-order operators into first-order 2k × 2k block systems [[0, I], [S, 0]].
Operator splitting for non-commutative matrix fields:: Lie (𝒪(dt)) and Strang (𝒪(dt²)) splitting for Sylvester-type equations, leveraging the fact that left and right operator actions commute as operations even when the underlying symbols do not.
Quasi-linear IMEX integration:: Rothe-type linearization for the conformal Ricci flow, freezing the non-linear coefficient once per step and applying the stiff, constant-coefficient remainder exactly via Fourier multipliers.
Robust time-loop management:: Built-in snapshot cadence, automatic finite-value checking to prevent silent divergence, and unified handling of scalar vs. matrix-valued initial conditions.

Mathematical background and numerical design

Asymptotic exponential symbol

The propagator exp(dt · P) is constructed via the truncated asymptotic series:

exp(tP) ~ I + tP + (t²/2!)P∘P + (t³/3!)P∘P∘P + …

where ∘ denotes the asymptotic composition (Kohn–Nirenberg or Weyl). For constant-coefficient symbols, all derivative corrections in the composition vanish, making the series exact at any truncation order (reducing to the standard matrix/scalar exponential of the symbol).

Amortized propagator families

Evaluating the asymptotic series for exp(tP) requires repeated calls to compose_asymptotic, which involves heavy symbolic differentiation. PropagatorFamily performs this expensive recursion exactly once, leaving t as a symbolic parameter _dt_family. Requesting a propagator for a specific numeric dt then requires only a single sympy.subs operation, bypassing the recursion entirely. This is cached process-wide by build_propagator.

Second-order block companion reduction

The equation ∂²ₜu = Op(S)(u) is rewritten as a first-order system for the state vector [u, v]ᵀ where v = ∂ₜu:

∂ₜ [u] [ 0 I ] [u]
[v] = [ S 0 ] [v]

For a k × k matrix symbol S, the companion matrix is 2k × 2k. This allows the full machinery of solve_first_order to be reused without deriving a separate second-order integrator.

Sylvester-type operator splitting

For ∂ₜU = Op(P)U − U·Op(Q), the left action L_P(U) = Op(P)U and right action R_Q(U) = U·Op(Q) commute as operations: L_P(R_Q(U)) = R_Q(L_P(U)). If P and Q are constant-coefficient (Fourier multipliers), the exact solution over a step dt is U(t) = exp(tP) U(0) exp(-tQ). If they depend on x, the underlying scalar operators Op[P_ij] and Op[Q_jk] may not commute. The solver then falls back to: - Lie-Trotter: U^{n+1} = exp(dt·P) exp(-dt·Q) U^n (𝒪(dt) error) - Strang: U^{n+1} = exp(dt/2·P) exp(-dt·Q) exp(dt/2·P) U^n (𝒪(dt²))

Quasi-linear IMEX for conformal Ricci flow

The 2D Ricci flow in conformal gauge g = e²ᵠ(dx² + dy²) reduces to the scalar quasi-linear heat equation ∂ₜφ = e⁻²ᵠΔφ. Because the coefficient c(x) = e⁻²ᵠ depends on the solution, no fixed symbol describes the operator ahead of time. Each step uses an IMEX/Lie splitting: 1. Explicit correction: Compute Δφ via a plain Laplacian. Take one

explicit Euler sub-step using the deviation of c from its spatial average c₀: φ_explicit = φ + dt · (c − c₀)Δφ.

Stiff step: Propagate the spatially averaged, x-independent generator c₀·Δ. Since c₀·Δ is a pure Fourier multiplier, exp(dt · c₀ · Δ) has the closed form exp(-dt · c₀ · |k|²) on the grid. This is applied directly via FFT, avoiding the need to rebuild an asymptotic exponential-symbol propagator at every step. This makes the stiff step exact for the frozen coefficient and significantly cheaper.

Numerical design notes

Grid generation strictly follows FFT ordering (np.fft.fftfreq) to ensure compatibility with the fast-path FFT multipliers in psiop_apply.
Time-stepping loops include automatic np.isfinite checks to raise FloatingPointError immediately upon divergence, rather than silently returning NaN-filled arrays.
Snapshot saving is decoupled from the integration step via a save_every parameter, allowing high-resolution temporal integration with low-resolution storage.
Matrix-valued solvers rigorously distinguish between left-action (apply_matrix_field) and right-action (apply_matrix_field_right), ensuring correct index contraction for density matrices and Green’s functions.

References

class psiop.psiop_solver.PropagatorFamily(s_expr, vars_x, order=3, quantization='kohn-nirenberg', mode_composition='kn', apply_backend='peetre', do_simplify=True)[source]

Bases: object

A one-step propagator family exp(dt · Op[s]) for a fixed symbol s but variable dt, built so that changing dt is cheap.

build_propagator rebuilds the full truncated exponential series (via exponential_symbol, i.e. order rounds of symbolic differentiation and compose_asymptotic) from scratch every time it is called – expensive, and wasteful if the same symbol is re-solved at several different step sizes (a dt-convergence study, a parameter sweep, or an adaptive/embedded stepper that changes dt every few steps).

exponential_symbol already accepts t as a sympy Symbol, so this class builds it ONCE with dt left symbolic, and produces a concrete propagator for any numeric dt via a single cheap sympy.subs call – skipping the compose_asymptotic recursion entirely on every subsequent request. build_propagator below uses this internally (via a small process-wide cache), so existing call sites benefit automatically with no code changes required.

Parameters:

s_expr (sympy.Expr or sympy.Matrix) – Symbol expression (scalar or matrix-valued), same as build_propagator.
vars_x (list of sympy.Symbol) – Spatial variables.
order (int, default 3) – Truncation order for both the Taylor series in t and the asymptotic composition at each power.
quantization ({'kohn-nirenberg', 'weyl'}, default 'kohn-nirenberg') – Quantization convention for the resulting propagator.
mode_composition ({'kn', 'weyl'}, default 'kn') – Composition rule used inside exponential_symbol.
apply_backend ({'peetre', 'direct'}, default 'peetre') – Numerical backend attached to the propagator operator.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Examples

>>> family = PropagatorFamily(xi**2, [x], order=4)
>>> prop_1, _, _ = family.propagator_for(0.01)
>>> prop_2, _, _ = family.propagator_for(0.02)   # cheap: no recomposition
>>> u1 = prop_1.apply(u, x_grid, kx)

propagator_for(dt)[source]

Return the concrete one-step propagator for a numeric dt.

Cheap: a single sympy.subs call – no compose_asymptotic recursion. Returns the same (prop, is_matrix, size) triple as build_propagator, so it is a drop-in replacement at the call site.

psiop.psiop_solver.block_matrix_second_order(s_expr)[source]

Convert a second-order-in-time operator symbol S into a first-order block companion system suitable for solve_first_order.

The second-order equation

∂²u/∂t² = Op[S](u)

is rewritten as the first-order system

∂/∂t [u] [ 0 I ] [u]
[v] = [ S 0 ] [v]

where v = ∂u/∂t. For a k×k matrix symbol S, the companion matrix has dimension 2k × 2k:

M = [ 0_k I_k ]
[ S 0_k ]

Parameters:: s_expr (sympy.Expr or sympy.Matrix) – The operator symbol S (scalar or k×k matrix).
Returns:: The block companion matrix M.
Return type:: sympy.Matrix, shape (2k, 2k)
Raises:: ValueError – If s_expr is a non-square matrix.

Examples

>>> M = block_matrix_second_order(-xi**2)
>>> M.shape
(2, 2)

psiop.psiop_solver.build_propagator(s_expr, vars_x, dt, order=3, quantization='kohn-nirenberg', mode_composition='kn', apply_backend='peetre', do_simplify=True)[source]

Build the one-step numerical propagator exp(dt · Op[s]) for a pseudo-differential operator via truncated asymptotic exponentiation.

Given a symbol s(x, ξ) (scalar or matrix-valued), this function:

Wraps it into the appropriate operator class.
Computes the symbol of exp(dt · P) via exponential_symbol(t=dt, order=order):

exp(dt · P) ≈ I + dt·P + (dt²/2!)·P∘P + ⋯ + (dtⁿ/n!)·P^{∘n}

where each power P^{∘n} is obtained through asymptotic composition.
Wraps the resulting symbol into a new operator ready for apply().

Parameters:

s_expr (sympy.Expr or sympy.Matrix) – Symbol expression. A scalar expression produces a PseudoDifferentialOperator; a matrix (or nested list) produces a MatrixPseudoDifferentialOperator.
vars_x (list of sympy.Symbol) – Spatial variables ([x] for 1D, [x, y] for 2D).
dt (float) – Time-step size used as the evolution parameter t in exp(t·P).
order (int, default 3) – Truncation order for both the Taylor series in t and the asymptotic composition at each power.
quantization ({'kohn-nirenberg', 'weyl'}, default 'kohn-nirenberg') – Quantization convention for the resulting propagator.
mode_composition ({'kn', 'weyl'}, default 'kn') – Composition rule used inside exponential_symbol.
apply_backend ({'peetre', 'direct'}, default 'peetre') – Numerical backend attached to the propagator operator.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Returns:

prop (PseudoDifferentialOperator or MatrixPseudoDifferentialOperator) – The propagator operator such that prop.apply(u, …) advances u by one time step dt.
is_matrix (bool) – True if the propagator is matrix-valued.
size (int or None) – Matrix dimension N if is_matrix is True, else None.

Notes

The propagator is constructed once and reused across all time steps. For constant-coefficient symbols the composition is exact (all derivative corrections vanish), so the only error is the Taylor truncation in dt.

This function is backed by a process-wide cache of PropagatorFamily objects (see above), keyed on everything except dt. So calling it repeatedly with the same s_expr/vars_x/order/quantization/ mode_composition/apply_backend but a different dt (e.g. a step-size sweep, or an adaptive stepper) reuses the expensive exponential_symbol/compose_asymptotic recursion instead of redoing it from scratch – only a cheap sympy.subs is performed per new dt. The very first call for a given symbol still pays the full cost.

Examples

>>> prop, is_mat, sz = build_propagator(xi**2, [x], dt=0.01, order=4)
>>> u_next = prop.apply(u, x_grid, kx)

psiop.psiop_solver.make_grid_1d(L=10.0, N=256)[source]

Construct a uniform periodic spatial grid and its associated FFT-ordered angular frequency grid in one dimension.

The spatial domain is [−L, L) discretised into N equally spaced points, and the frequency grid covers the discrete angular wavenumbers compatible with the FFT ordering:

x_j = −L + j·Δx, j = 0, …, N−1, Δx = 2L/N k_m = 2π · fftfreq(N, Δx)

Parameters:

L (float, default 10.0) – Half-length of the spatial domain. The full period is 2L.
N (int, default 256) – Number of grid points.

Returns:

x (ndarray, shape (N,)) – Spatial coordinates in [−L, L).
kx (ndarray, shape (N,)) – Angular frequency grid in FFT order (radians per unit length).

Examples

>>> x, kx = make_grid_1d(L=5.0, N=128)
>>> x[0], x[-1]
(-5.0, 4.921875)

psiop.psiop_solver.make_grid_2d(L=10.0, N=128)[source]

Construct uniform periodic spatial grids and their associated FFT-ordered angular frequency grids in two dimensions.

Both axes share the same half-length L and resolution N, yielding a square domain [−L, L)² with N × N grid points.

x_i = −L + i·Δx, y_j = −L + j·Δy, Δx = Δy = 2L/N kx_m = 2π · fftfreq(N, Δx), ky_n = 2π · fftfreq(N, Δy)

Parameters:

L (float, default 10.0) – Half-length of the spatial domain along each axis.
N (int, default 128) – Number of grid points per axis.

Returns:

x (ndarray, shape (N,)) – Spatial coordinates along x.
y (ndarray, shape (N,)) – Spatial coordinates along y.
kx (ndarray, shape (N,)) – Angular frequency grid along x in FFT order.
ky (ndarray, shape (N,)) – Angular frequency grid along y in FFT order.

Examples

>>> x, y, kx, ky = make_grid_2d(L=5.0, N=64)

psiop.psiop_solver.make_grids(vars_x, L, N)[source]

Build spatial + frequency grids and the meshgrid-ed spatial coordinates used to evaluate initial conditions – factors out the grid-setup boilerplate that used to be copy-pasted verbatim in every solve_* function below.

Returns:

X, Y (ndarray, ndarray or None) – Meshgrid-ed spatial coordinates (‘ij’ indexing). Y is None in 1D.
x_grid, y_grid (ndarray, ndarray or None) – 1D spatial axes (y_grid is None in 1D).
kx, ky (ndarray, ndarray or None) – Frequency axes (ky is None in 1D).
grids (tuple) – (x, kx) in 1D or (x, y, kx, ky) in 2D – what callers return.

psiop.psiop_solver.run_time_loop(step_fn, U0, dt, n_steps, save_every, check_finite=True)[source]

Repeatedly apply step_fn(U) -> U_next, saving a snapshot every save_every steps (plus the final step and t=0) – factors out the save-cadence bookkeeping that used to be copy-pasted verbatim in every solve_* function below.

Raises:: FloatingPointError – If check_finite is True and a non-finite value shows up – avoids silently returning a diverged/garbage trajectory.

psiop.psiop_solver.solve_first_order(s_expr, vars_x, f, dt, n_steps, order=3, L=10.0, N=256, apply_kwargs=None, save_every=1, quantization='kohn-nirenberg', apply_backend='peetre', check_finite=True, do_simplify=True)[source]

Solve the first-order evolution equation

∂u/∂t = Op[s](u), u(x, 0) = f(x)

by repeated application of the asymptotic propagator exp(dt · Op[s]).

At each time step the field is advanced via

u^{n+1} = exp(dt · Op[s]) u^n ≈ (I + dt·P + (dt²/2!)P∘P + ⋯) u^n

Parameters:

s_expr (sympy.Expr or sympy.Matrix) – Symbol of the spatial operator P. Scalar for a single-field equation; matrix-valued for a coupled system.
vars_x (list of sympy.Symbol) – Spatial variables.
f (callable) – Initial condition. Must accept (X,) in 1D or (X, Y) in 2D and return an ndarray (scalar case) or a list/tuple of ndarrays (matrix case with N components).
dt (float) – Time-step size.
n_steps (int) – Total number of time steps to evolve.
order (int, default 3) – Asymptotic expansion order for the propagator construction.
L (float, default 10.0) – Half-length of the periodic spatial domain.
N (int, default 256) – Number of grid points per spatial axis.
apply_kwargs (dict, optional) – Extra keyword arguments forwarded to prop.apply() at every step (e.g. boundary_condition, freq_window, clamp).
save_every (int, default 1) – Store the solution snapshot every save_every steps.
quantization (str, default 'kohn-nirenberg') – Quantization convention.
apply_backend (str, default 'peetre') – Numerical application backend.
check_finite (bool, default True) – Raise FloatingPointError as soon as a NaN/Inf appears, instead of silently returning a diverged trajectory.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Returns:

t (ndarray, shape (n_saved,)) – Time values at which snapshots were recorded.
U (ndarray) – Solution snapshots. Shape (n_saved, N) for scalar 1D, (n_saved, N, N) for scalar 2D, or (n_saved, size, N…) for matrix-valued systems.
grids (tuple) – The spatial and frequency grids used: (x, kx) in 1D or (x, y, kx, ky) in 2D.

Raises:

NotImplementedError – If vars_x has length other than 1 or 2.
ValueError – If f returns the wrong number of components for a matrix system.

Examples

>>> t, U, (x, kx) = solve_first_order(xi**2, [x], lambda X: np.exp(-X**2),
...                                     dt=0.01, n_steps=100, N=256)

psiop.psiop_solver.solve_matrix_field(s_expr, vars_x, F, dt, n_steps, order=3, L=10.0, N=256, apply_kwargs=None, save_every=1, quantization='kohn-nirenberg', apply_backend='peetre', check_finite=True, do_simplify=True)[source]

Time-step the matrix-field evolution equation ∂ₜU = P U, where P is the pseudo-differential operator with N×N matrix symbol s_expr and U(x) is itself an N×N matrix at every spatial point (e.g. a density matrix or matrix Green’s function), with P acting on U only from the left: (P U)_ik = Σⱼ Op[P_ij](U_jk). This repeatedly applies the exponential propagator Op(exp(dt·s)) built by build_propagator, via MatrixPseudoDifferentialOperator.apply_matrix_field, exactly as solve does for vector fields via apply.

Parameters:

s_expr (sympy.MatrixBase or nested list of sympy.Expr) – N×N matrix symbol S(x, ξ) of the generator P; must be matrix-valued (matrix left-multiplication only makes sense at N > 1 – use solve for a scalar generator).
vars_x (list of sympy symbols) – Spatial variables (length 1 or 2).
F (callable) – Initial matrix field U(·, 0), called as F(X) in 1D or F(X, Y) in 2D on the meshgrid-ed spatial coordinates, and expected to return an N×N array/nested list of grid-shaped components (F(…)[j][k], or an ndarray of shape (N, N, *grid_shape)).
dt (float) – Time step.
n_steps (int) – Number of propagator applications (time steps) to take.
order (int, optional) – Truncation order of the exponential symbol expansion. Default 3.
L (float, optional) – Half-width of the spatial domain. Default 10.0.
N (int, optional) – Number of grid points per axis. Default 256.
apply_kwargs (dict, optional) – Extra keyword arguments forwarded to apply_matrix_field.
save_every (int, optional) – Save the solution every save_every steps (plus the final step and t=0). Default 1 (save every step).
quantization (str, optional) – Quantization convention. Default ‘kohn-nirenberg’.
apply_backend (str, optional) – Numerical application backend. Default ‘peetre’.
check_finite (bool, default True) – Raise FloatingPointError as soon as a NaN/Inf appears, instead of silently returning a diverged trajectory.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Returns:

t_list (ndarray) – Saved time points, starting at 0.
U_list (ndarray, shape (n_saved, N, N, *grid_shape)) – Saved matrix-field snapshots U(t).
grids (tuple of ndarray) – (x, kx) in 1D or (x, y, kx, ky) in 2D, as returned by make_grid_1d/make_grid_2d.

Raises:

NotImplementedError – If vars_x has a length other than 1 or 2.
ValueError – If s_expr is not matrix-valued, or F does not return an N×N field.

psiop.psiop_solver.solve_ricci_flow_conformal_2d(phi0, dt, n_steps, order=3, L=8.0, N=64, save_every=1, quantization='kohn-nirenberg', apply_backend='peetre', check_finite=True)[source]

Integrate 2D Ricci flow in conformal gauge on a flat, doubly periodic background.

Writing the metric as g = e^{2φ}(dx² + dy²), the Gauss curvature is K = −e^{−2φ}Δφ and, since R_ij = K·g_ij in two dimensions, the tensorial flow ∂ₜg_ij = −2R_ij collapses to the scalar quasi-linear heat equation

∂ₜφ = e^{−2φ} Δφ ,

with Δ = ∂ₓ² + ∂ᵧ² the flat Laplacian. This is NOT handled by psiop’s ordinary linear/matrix machinery: the coefficient e^{−2φ} depends on the evolving solution itself, so no fixed sympy symbol p(x, ξ) describes the operator ahead of time.

Instead, each step uses an IMEX/Lie splitting that still reuses psiop’s exact exponential propagator for the stiff part:

Explicit correction — using the current coefficient field c(x) = e^{−2φ(x)}, compute Δφ once via a plain (x-independent) Laplacian PseudoDifferentialOperator, and take one explicit Euler sub-step with the deviation of c from its spatial average c₀: residual = (c − c₀) Δφ.
Stiff step — propagate the spatially averaged, x-independent (Fourier-multiplier) generator c₀·Δ applied to the explicitly-corrected field. Since c0 * lap_symbol is a pure Fourier multiplier, exp(dt · c0 · Δ) has the closed form exp(-dt · c0 · |k|^2) on this grid – applied directly via FFT instead of rebuilding an asymptotic exponential-symbol propagator through build_propagator on every step. This is both cheaper (no symbolic recomputation per step) and more accurate (exact exponential instead of an order-n Taylor truncation in dt); order is kept only for backward compatibility and no longer affects the stiff step.

This freezes the quasi-linear coefficient once per step (a Rothe-type linearization), so accuracy in dt is limited by that freezing, not by the stiff step itself, which is now exact for the frozen (constant-coefficient) part at every step.

Parameters:

phi0 (callable) – Initial conformal factor, called as phi0(X, Y) on the meshgrid-ed spatial coordinates; must return a real-valued array.
dt (float) – Time step.
n_steps (int) – Number of steps to take.
order (int, optional) – Kept for backward compatibility; no longer affects the (now closed-form) stiff step. Default 3.
L (float, optional) – Half-width of the (periodic) spatial domain along each axis. Default 8.0.
N (int, optional) – Number of grid points per axis. Default 64.
save_every (int, optional) – Save every save_every steps (plus the final step and t=0). Default 1.
quantization (str, optional) – Quantization convention used for the explicit Laplacian evaluation. Default ‘kohn-nirenberg’.
apply_backend (str, optional) – Numerical application backend for the explicit Laplacian evaluation. Default ‘peetre’.
check_finite (bool, default True) – Raise FloatingPointError as soon as a NaN/Inf appears, instead of silently returning a diverged trajectory.

Returns:

t_list (ndarray) – Saved time points, starting at 0.
phi_list (ndarray, shape (n_saved, N, N)) – Saved conformal-factor snapshots φ(x, y, t); the metric at each saved time is g(t) = e^{2·phi_list[k]} (dx² + dy²).
grids (tuple of ndarray) – (x, y) spatial grids, as returned by make_grid_2d.

psiop.psiop_solver.solve_second_order(s_expr, vars_x, f, g, dt, n_steps, order=3, L=10.0, N=256, apply_kwargs=None, save_every=1, quantization='kohn-nirenberg', apply_backend='peetre', do_simplify=True, scheme='propagator', check_finite=True)[source]

Solve the second-order-in-time evolution equation

∂²u/∂t² = Op[S](u), u(x, 0) = f(x), ∂u/∂t(x, 0) = g(x)

The system is split into ∂u/∂t = v, ∂v/∂t = Op[S](u), and time-stepped according to scheme:

'propagator' (default): reduce to a first-order 2k-dimensional block companion system and time-step it via solve_first_order, i.e. repeated application of the amortized asymptotic/exact exponential propagator. For constant-coefficient S this is exact at any dt (unconditionally stable, zero phase error). For x-dependent S it is accurate to the asymptotic truncation order, but the truncated exponential is not symplectic in general, so very long integrations can accumulate a slow energy drift.
'leapfrog': explicit kick-drift-kick (Störmer-Verlet) scheme, applying Op[S] directly – no exponential/compose_asymptotic machinery is built, so order and do_simplify are ignored. It is only conditionally stable: writing the local frequency symbol as ω(x, ξ) = sqrt(-S(x, ξ)), dt must roughly satisfy dt · max ω < 2 (a violation typically shows up immediately as the FloatingPointError below). In exchange, being symplectic and time-reversible, its energy error stays bounded over exponentially long integration times instead of drifting secularly – the main reason to prefer it over ‘propagator’ for long-time runs with genuinely x-dependent S. Plain forward Euler is deliberately not offered as a scheme here: applied to this companion system (purely imaginary eigenvalues for an oscillatory S) its amplification factor exceeds 1 for every dt > 0, so it is unconditionally unstable; leapfrog is the cheapest explicit scheme that avoids that failure mode.

Per-step KDK update, from (u^n, v^n) to (u^{n+1}, v^{n+1}):

v^{n+1/2} = v^n + (dt/2)·Op[S](u^n) # kick u^{n+1} = u^n + dt ·v^{n+1/2} # drift v^{n+1} = v^{n+1/2} + (dt/2)·Op[S](u^{n+1}) # kick

equivalent to the textbook u^{n+1} = 2u^n − u^{n−1} + dt²·Op[S](u^n) form but self-starting from (u^0, v^0) with no separate initialization step, and it yields a velocity that is naturally synchronized with u at every saved time (the average of the two half-kicks straddling it).

Parameters:

s_expr (sympy.Expr or sympy.Matrix) – Symbol of the spatial operator S.
vars_x (list of sympy.Symbol) – Spatial variables.
f (callable) – Initial displacement u(x, 0). Signature: f(X) in 1D, f(X, Y) in 2D. Must return one component per row of S.
g (callable) – Initial velocity ∂u/∂t(x, 0). Same signature and component structure as f.
dt (float) – Time-step size.
n_steps (int) – Number of time steps.
order (int, default 3) – Asymptotic order for the propagator. Ignored if scheme=’leapfrog’.
L (float, default 10.0) – Spatial domain half-length.
N (int, default 256) – Grid points per axis.
apply_kwargs (dict, optional) – Forwarded to apply() at each step.
save_every (int, default 1) – Snapshot cadence.
quantization (str, default 'kohn-nirenberg') – Quantization convention.
apply_backend (str, default 'peetre') – Numerical backend.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree. Ignored if scheme=’leapfrog’.
scheme ({'propagator', 'leapfrog'}, default 'propagator') – Time-integration scheme; see above.
check_finite (bool, default True) – Raise FloatingPointError as soon as a NaN/Inf appears, instead of silently returning a diverged trajectory. Under scheme=’leapfrog’ this is usually how a CFL violation first manifests.

Returns:

t (ndarray, shape (n_saved,)) – Time values of recorded snapshots.
U (ndarray) – Displacement field snapshots u(x, t).
V (ndarray) – Velocity field snapshots ∂u/∂t(x, t).
grids (tuple) – Spatial and frequency grids: (x, kx) or (x, y, kx, ky).

Raises:

ValueError – If f or g produce the wrong number of components, or if scheme is not one of ‘propagator’/’leapfrog’.

Examples

>>> t, U, V, (x, kx) = solve_second_order(
...     -xi**2, [x], lambda X: np.exp(-X**2), lambda X: 0*X,
...     dt=0.01, n_steps=200, scheme='leapfrog')

psiop.psiop_solver.solve_sylvester_field(P_expr, Q_expr, vars_x, F, dt, n_steps, order=3, splitting='strang', L=10.0, N=256, apply_kwargs=None, save_every=1, quantization='kohn-nirenberg', apply_backend='peetre', check_finite=True, do_simplify=True)[source]

Time-step the Sylvester-type matrix-field evolution equation ∂ₜU = P U − U Q, where P and Q are pseudo-differential operators with N×N matrix symbols and U(x) is an N×N matrix at every spatial point.

Left-multiplication by P and right-multiplication by Q always commute as operations ((P U) Q == P (U Q)), so when P and Q are x-independent (Fourier multipliers), the exact solution over a step dt is the closed-form

U(t) = exp(t P) U(0) exp(-t Q) ,

obtained by applying the left-propagator Op(exp(dt·P)) (apply_matrix_field) and the right-propagator Op(exp(-dt·Q)) (apply_matrix_field_right), in either order. When P and/or Q depend on x, Op[P_ij] and Op[Q_jk] need not commute with each other, so the two sub-steps no longer combine exactly; this function then falls back to a standard Lie-Trotter (first-order, O(dt) splitting error) or Strang (second-order, O(dt^2)) operator splitting between the left and right exponential propagators.

Parameters:

P_expr (sympy.MatrixBase or nested list of sympy.Expr) – N×N matrix symbol P(x, ξ) acting on U from the left.
Q_expr (sympy.MatrixBase or nested list of sympy.Expr) – N×N matrix symbol Q(x, ξ) acting on U from the right (with a minus sign, as in ∂ₜU = P U − U Q); must be the same size as P_expr.
vars_x (list of sympy symbols) – Spatial variables (length 1 or 2).
F (callable) – Initial matrix field U(·, 0), called as F(X) in 1D or F(X, Y) in 2D, returning an N×N array/nested list of grid-shaped components (as for solve_matrix_field).
dt (float) – Time step.
n_steps (int) – Number of splitting steps (each advancing U by dt).
order (int, optional) – Truncation order of each exponential-symbol expansion. Default 3.
splitting (str, {'lie', 'strang'}, optional) –
Operator-splitting scheme between the P (left) and Q (right) sub-steps:
- ’lie’one full left step exp(dt·P), then one full right
  step exp(-dt·Q) – first order accurate, O(dt).
- ’strang’half left step exp(dt/2·P), full right step
  exp(-dt·Q), half left step exp(dt/2·P) – second order accurate, O(dt^2). Default.
L (float, optional) – Half-width of the spatial domain. Default 10.0.
N (int, optional) – Number of grid points per axis. Default 256.
apply_kwargs (dict, optional) – Extra keyword arguments forwarded to apply_matrix_field and apply_matrix_field_right.
save_every (int, optional) – Save the solution every save_every steps (plus the final step and t=0). Default 1.
quantization (str, optional) – Quantization convention. Default ‘kohn-nirenberg’.
apply_backend (str, optional) – Numerical application backend. Default ‘peetre’.
check_finite (bool, default True) – Raise FloatingPointError as soon as a NaN/Inf appears, instead of silently returning a diverged trajectory.
do_simplify (bool, default True) – Whether to call sympy’s simplify() while assembling the propagator symbol (once when building it, and once inside every compose_asymptotic() call in the asymptotic expansion loop). This does not change the operator being applied – lambdify evaluates the same function on an unsimplified expression – it only affects how much symbolic cleanup happens before that. simplify() is the dominant cost of build_propagator() for symbols mixing trigonometric and polynomial terms, and its cost grows with order; set to False to skip it and speed up propagator construction, at the risk of a larger (but numerically equivalent) unsimplified expression tree.

Returns:

t_list (ndarray) – Saved time points, starting at 0.
U_list (ndarray, shape (n_saved, N, N, *grid_shape)) – Saved matrix-field snapshots U(t).
grids (tuple of ndarray) – (x, kx) in 1D or (x, y, kx, ky) in 2D.

Raises:

NotImplementedError – If vars_x has a length other than 1 or 2.
ValueError – If P_expr/Q_expr are not matrix-valued of matching size, if F does not return an N×N field, or if splitting is not ‘lie’ or ‘strang’.