microlocal

microlocal.py — Unified microlocal analysis toolkit for 1D and 2D problems

Overview

The microlocal module provides a high‑level interface for studying the propagation of singularities and constructing semiclassical approximations for linear partial differential equations. It builds upon the dedicated wkb (WKB approximations) and caustics (catastrophe classification and ray caustic detection) modules, adding dimension‑agnostic functions for the core concepts of microlocal analysis:

• Characteristic variety Char(P) = {(x,ξ) : p(x,ξ)=0} – the set of phase‑space points where the principal symbol vanishes, indicating where singularities can propagate.

• Bicharacteristic flow – Hamilton’s equations for the principal symbol, whose integral curves (bicharacteristics) govern the propagation of wave‑front sets.

• Wavefront set WF(u) – a refinement of the singular support that also records the directions (frequencies) in which the singularity occurs. The module visualises how an initial wavefront set evolves under the flow.

• WKB approximation (re‑exported from wkb) – asymptotic solutions of the form u ≈ A e^{iS/ε}.

• Caustics and Maslov index – detection of caustics (where rays focus) and computation of the associated Maslov phase shifts, crucial for correct semiclassical quantisation.

• Bohr–Sommerfeld quantisation (1D) – semiclassical energy levels for bound states.

• Operator visualisation and interactive analysis – a comprehensive suite for plotting symbol amplitude, phase, micro‑support, characteristic gradients, and launching interactive ipywidgets dashboards.

• Matrix-valued symbols and pseudospectra – extraction of characteristic branches via eigenvalue computation for matrix operators, and visualisation of ε‑pseudospectra.

• PDE solution rendering – rendering and animation tools for scalar, matrix, and wave equation fields in 1D and 2D space-time.

All functions automatically detect the spatial dimension (1 or 2) from the input data, making the module usable for both one‑dimensional and two‑dimensional problems without changing the calling syntax.

Mathematical background

In microlocal analysis, a linear partial differential operator P is studied via its principal symbol p(x,ξ), a function on the cotangent bundle T*ℝⁿ. The characteristic variety is the zero set of p. Singularities of a distribution u satisfying P u ≈ 0 are confined to the characteristic variety and propagate along bicharacteristics – integral curves of the Hamiltonian vector field

X_p = ( ∂p/∂ξ , –∂p/∂x ).

The wavefront set WF(u) is a closed conic subset of T*ℝⁿ {0} that records both the location x and the direction ξ of the singularity. If (x₀,ξ₀) ∉ WF(u), then u is smooth in a neighbourhood of x₀ in the direction ξ₀. The fundamental theorem of microlocal analysis states that WF(Pu) ⊆ WF(u) and that WF(u) WF(Pu) is contained in the characteristic variety and is invariant under the bicharacteristic flow.

The WKB method seeks solutions of the form u(x) = e^{iS(x)/ε} (a₀(x) + ε a₁(x) + …). The phase S satisfies the eikonal equation p(x,∇S)=0, and the amplitudes a_k satisfy transport equations along bicharacteristics. This construction breaks down at caustics, where rays focus; the Maslov index μ (a signed count of caustic crossings) provides a phase correction e^{iμπ/2} that restores uniformity.

The module integrates these concepts into a coherent toolkit, allowing the user to:

• Symbolically compute the characteristic variety.

• Numerically integrate bicharacteristics with symplectic integrators.

• Visualise the evolution of wavefront sets, including 2D and 3D singularity animations.

• Obtain semiclassical spectra via Bohr–Sommerfeld quantisation (1D).

• Detect caustics and compute the Maslov index (using the caustics module).

• Interactively explore symbol properties such as amplitude, phase, and micro-support.

• Analyse matrix-valued symbols by computing their characteristic branches.

• Plot ε-pseudospectra and overlay complex eigenvalues.

• Render scalar, matrix, and wave PDE solutions in space-time.

References

[1]

Hörmander, L. The Analysis of Linear Partial Differential Operators I, Springer, 1983. Chapter 8: Wave Front Sets.

[2]

Duistermaat, J. J. Fourier Integral Operators, Courant Institute Lecture Notes, 1996.

[3]

Maslov, V. P. & Fedoriuk, M. V. Semi‑Classical Approximation in Quantum Mechanics, Reidel, 1981.

[4]

Zworski, M. Semiclassical Analysis, American Mathematical Society, 2012. Chapter 3: Propagation of Singularities.

[5]

Taylor, M. E. Partial Differential Equations II, Springer, 2011. Chapter 8: Microlocal Analysis.

microlocal.animate_operator_singularity(op, xi0=5.0, eta0=0.0, x0=0.0, y0=0.0, tmax=4.0, n_frames=100, projection=None)

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).

microlocal.animate_scalar_1d(t, U, x, quantity='real', interval=40, save_path=None)

Animate a scalar 1D solution as a line plot evolving in time.

Parameters:

• t (ndarray, shape (n_times,)) – Time samples.

• U (ndarray, shape (n_times, Nx)) – Solution values u(x, t) sampled on the grid.

• x (ndarray, shape (Nx,)) – Spatial grid.

• quantity ({'real', 'imag', 'abs'}, default='real') – Which part of U to plot; also sets the fixed y-axis limits from the min/max of that quantity over the whole trajectory.

• interval (int, default=40) – Delay between animation frames, in milliseconds.

• save_path (str, optional) – If given, the animation is saved to this path via anim.save.

Returns:

The animation object (figure is not shown automatically).

Return type:

matplotlib.animation.FuncAnimation

microlocal.animate_singularity(s_expr, vars_x, x0=0.0, xi0=5.0, tmax=4.0, n_frames=100, projection=None, branches='all', labels=None, interval=50, contours=True, solution=None, quantity='abs', save_path=None)

Animate the propagation of singularities along bicharacteristic trajectories, projected onto a 2D phase-space plane.

Internally calls integrate_singularity to obtain the trajectories, then draws a growing dashed trail plus a moving point per branch via the shared _trail_animation helper.

Parameters:

• s_expr (sympy.Expr) – Principal symbol used to build the Hamiltonian(s).

• vars_x (list of sympy.Symbol) – Spatial variables.

• x0 (float or array_like, default=0.0) – Initial spatial position(s).

• xi0 (float or array_like, default=5.0) – Initial frequency (momentum) component(s).

• tmax (float, default=4.0) – Final integration time.

• n_frames (int, default=100) – Number of time samples used for the trajectory and the animation.

• projection ({'phase', 'position', 'frequency'} or None, default=None) – Which plane to draw. In 1D, defaults to ‘phase’ (x vs xi); in 2D, defaults to ‘position’ (x vs y) and (‘frequency’/’phase’ are not selectable in 2D – the projection is always (x, y)).

• branches ('all', int, or list of int, default='all') – Which characteristic branches to animate.

• labels (list of str, optional) – Currently unused (reserved for per-branch legend labels).

• interval (int, default=50) – Delay between animation frames, in milliseconds.

• contours (bool, default=True) – Currently unused (reserved for background contour overlays).

• solution (optional) – Currently unused (reserved for overlaying a PDE solution field).

• quantity (str, default='abs') – Currently unused (reserved alongside solution).

• save_path (str, optional) – If given, the animation is saved to this path via anim.save.

Returns:

The animation object.

Return type:

matplotlib.animation.FuncAnimation

Raises:

ValueError – If projection is not one of the supported values for a 1D symbol.

microlocal.animate_singularity_3d(s_expr, vars_x, x0=0.0, xi0=5.0, tmax=4.0, n_frames=100, projection=None, branches='all', labels=None, interval=50, save_path=None)

Animate bicharacteristic trajectories in a 3D matplotlib plot.

Internally calls integrate_singularity to obtain the trajectories, then draws a growing dashed trail plus a moving point per branch via the shared _trail_animation helper. In 1D, the third axis is time t; in 2D (or higher), the first three phase-space coordinates (x, y, …) are used directly.

Parameters:

• s_expr (sympy.Expr) – Principal symbol used to build the Hamiltonian(s).

• vars_x (list of sympy.Symbol) – Spatial variables.

• x0 (float or array_like, default=0.0) – Initial spatial position(s).

• xi0 (float or array_like, default=5.0) – Initial frequency (momentum) component(s).

• tmax (float, default=4.0) – Final integration time.

• n_frames (int, default=100) – Number of time samples used for the trajectory and the animation.

• projection (optional) – Currently unused for the 3D case (reserved for API parity with animate_singularity); the axes are always chosen as described above.

• branches ('all', int, or list of int, default='all') – Which characteristic branches to animate.

• labels (list of str, optional) – Currently unused (reserved for per-branch legend labels).

• interval (int, default=50) – Delay between animation frames, in milliseconds.

• save_path (str, optional) – If given, the animation is saved to this path via anim.save.

Returns:

The animation object.

Return type:

matplotlib.animation.FuncAnimation

microlocal.bicharacteristic_flow(symbol, z0, tspan, dim=None, method='symplectic', n_steps=1000)

Integrate the bicharacteristic flow on the cotangent bundle.

Hamilton’s equations: ẋ = ∂p/∂ξ, ξ̇ = -∂p/∂x (1D) or ẋ = ∂p/∂ξ, ẏ = ∂p/∂η, ξ̇ = -∂p/∂x, η̇ = -∂p/∂y (2D)

Parameters:

• symbol (sympy expression) – Principal symbol.

• z0 (tuple) – Initial condition on T*M. For 1D: (x₀, ξ₀); for 2D: (x₀, y₀, ξ₀, η₀).

• tspan (tuple) – (t_start, t_end).

• dim (int, optional) – Dimension. If None, inferred from length of z0.

• method (str) – Integration method: ‘rk45’, ‘symplectic’, or ‘verlet’ (2D only).

• n_steps (int) – Number of time steps.

Returns:

Trajectory data with keys ‘t’, ‘x’, ‘xi’ (1D) and also ‘y’,’eta’ (2D), plus ‘symbol_value’.

Return type:

dict

microlocal.bohr_sommerfeld_quantization(H, n_max=10, x_range=(-10, 10), hbar=1.0, E_range=(1e-06, 50.0))

Compute Bohr-Sommerfeld quantization for 1D bound states.

Solves (1/(2π)) ∮ p dx = ℏ(n + α) with α = 1/2 (Maslov index).

Parameters:

• H (sympy.Expr) – Hamiltonian H(x,p) as a symbolic expression.

• n_max (int, optional) – Maximum quantum number to compute. Default is 10.

• x_range (tuple of float, optional) – Spatial range for finding turning points. Default is (-10, 10).

• hbar (float, optional) – Reduced Planck constant. Default is 1.0.

• E_range (tuple of float, optional) – Energy range to scan for quantized levels. Default is (1e-6, 50.0).

Returns:

A dictionary containing: - ‘n’ : numpy.ndarray (Quantum numbers) - ‘E_n’ : numpy.ndarray (Quantized energy levels) - ‘actions’ : numpy.ndarray (Corresponding classical actions) - ‘hbar’ : float (The reduced Planck constant used) - ‘alpha’ : float (The Maslov index correction used, 0.5)

Return type:

dict

microlocal.characteristic_hamiltonians(s_expr, vars_x, vars_xi=None)

Extract the characteristic Hamiltonian functions H(x, ξ) from a (possibly matrix-valued) operator symbol.

For a scalar symbol p(x, ξ), the single Hamiltonian is

H(x, ξ) = Re(p(x, ξ))

For a matrix symbol P(x, ξ), the eigenvalues λ_k(x, ξ) are computed symbolically and each branch yields

H_k(x, ξ) = Re(λ_k(x, ξ))

These Hamiltonians generate the bicharacteristic (ray) flow via Hamilton’s equations:

ẋ = ∂H/∂ξ, ξ̇ = −∂H/∂x

Parameters:

• s_expr (sympy.Expr or sympy.Matrix) – Operator symbol (scalar or matrix-valued).

• vars_x (list of sympy.Symbol) – Spatial variables.

• vars_xi (list of sympy.Symbol, optional) – Frequency variables. If None, they are inferred from the free symbols of s_expr that are not in vars_x, and ordered canonically via _order_freq_vars.

Returns:

• H_list (list of sympy.Expr) – One Hamiltonian per characteristic branch (eigenvalue).

• xs (list of sympy.Symbol) – Canonical spatial symbols (real=True) used in H_list.

• xis (list of sympy.Symbol) – Canonical frequency symbols (real=True) used in H_list.

Notes

The classical Hamiltonian governing bicharacteristic flow is the (real part of the) principal symbol itself. A previous version of this function erroneously computed Re(i * p), which vanishes identically for any real-valued symbol and therefore produced trivial (stationary) trajectories in the flow visualization.

The substitution to fresh canonical symbols ensures consistent differentiation even if the input expression uses symbols with different assumptions.

microlocal.characteristic_variety(symbol, dim=None, tol=1e-08)

Compute the characteristic variety of a pseudo-differential operator.

Char(P) = { (x,ξ) in T*ℝ : p(x,ξ)=0 } (1D) or Char(P) = { (x,y,ξ,η) in T*ℝ² : p(x,y,ξ,η)=0 } (2D)

Parameters:

• symbol (sympy expression) – Principal symbol p.

• dim (int, optional) – Dimension (1 or 2). If None, inferred from symbol.

• tol (float) – Tolerance for zero detection (unused, kept for compatibility).

Returns:

A dictionary containing the following keys:

• ’implicit’sympy expression

  The symbol expression.

• ’equation’sympy Eq

  The equation symbol = 0.

• ’explicit’list or None

  List of explicit solutions ξ(x) (1D only), or None if no explicit solution exists.

• ’function’callable

  A callable that evaluates the symbol.

Return type:

dict

microlocal.compute_caustics_2d(p, initial_curve, tmax, n_rays=None, **kwargs)

Compute caustics for a 2D Hamiltonian given an initial curve.

Parameters:

• p (sympy expression) – Hamiltonian symbol p(x, y, xi, eta).

• initial_curve (dict) – Must contain keys ‘x’, ‘y’, ‘xi’, ‘eta’ with array values of equal length.

• tmax (float) – Maximum integration time.

• n_rays (int, optional) – Number of rays to use. If None, use all points in the initial curve.

• **kwargs (additional arguments passed to bicharacteristic_flow) – (e.g., method=’symplectic’, n_steps=200).

Returns:

Each event contains position, time, momentum, Arnold type, etc.

Return type:

list of CausticEvent

microlocal.compute_maslov_index(traj)

Compute the Maslov index for a single trajectory.

Parameters:

traj (dict) – Returned by bicharacteristic_flow. Must contain J11..J22 components.

Returns:

Maslov index (number of sign changes of det(J) or caustic crossings).

Return type:

int

microlocal.find_caustics_1d(symbol, x_range, xi_range, resolution=100)

Find caustics (envelope of bicharacteristics) in 1D.

Uses a simplified condition: identifies regions where d²p/dξ² ≈ 0, which correspond to turning points in frequency.

Parameters:

• symbol (sympy.Expr) – The 1D principal symbol p(x, ξ).

• x_range (tuple of float) – The (min, max) range for the spatial variable x.

• xi_range (tuple of float) – The (min, max) range for the frequency variable ξ.

• resolution (int, optional) – Number of grid points per axis for the evaluation mesh. Default is 100.

Returns:

A dictionary containing: - ‘x_grid’ : numpy.ndarray (2D meshgrid of x coordinates) - ‘xi_grid’ : numpy.ndarray (2D meshgrid of ξ coordinates) - ‘caustic_indicator’ : numpy.ndarray (Absolute value of d²p/dξ²) - ‘threshold’ : float (10th percentile of the indicator, useful for masking)

Return type:

dict

microlocal.group_velocity_field(op, xlim=(-2, 2), klim=(-10, 10), density=30)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

Raises:

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

Notes

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

• Used for analyzing wave propagation properties and dispersion relations.

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

microlocal.integrate_singularity(s_expr, vars_x, x0=0.0, xi0=5.0, tmax=4.0, n_frames=100, vars_xi=None, branches='all', method='RK45', **ivp_kwargs)

Numerically integrate bicharacteristic (Hamiltonian ray) trajectories from an initial phase-space point (x₀, ξ₀).

For each characteristic branch H_k, Hamilton’s equations

ẋ = ∂H_k/∂ξ, ξ̇ = −∂H_k/∂x

are integrated using scipy.integrate.solve_ivp over [0, tmax].

Parameters:

• s_expr (sympy.Expr or sympy.Matrix) – Operator symbol from which Hamiltonians are extracted.

• vars_x (list of sympy.Symbol) – Spatial variables.

• x0 (float or array_like, default 0.0) – Initial spatial position(s). Scalar for 1D, sequence for 2D.

• xi0 (float or array_like, default 5.0) – Initial frequency (momentum) component(s).

• tmax (float, default 4.0) – Final integration time.

• n_frames (int, default 100) – Number of output time samples in [0, tmax].

• vars_xi (list of sympy.Symbol, optional) – Explicit frequency variables (inferred if None).

• branches ('all', int, or list of int, default 'all') – Which characteristic branches to integrate. ‘all’ integrates every branch; an int or list selects specific ones.

• method (str, default 'RK45') – ODE solver method passed to solve_ivp.

• **ivp_kwargs – Additional keyword arguments forwarded to solve_ivp (e.g. rtol, atol, max_step).

Returns:

• H_list (list of sympy.Expr) – Hamiltonian expressions for the integrated branches.

• xs (list of sympy.Symbol) – Canonical spatial symbols.

• xis (list of sympy.Symbol) – Canonical frequency symbols.

• t_eval (ndarray, shape (n_frames,)) – Time samples at which trajectories are recorded.

• trajs (list of ndarray) – trajs[b] has shape (2·dim, n_frames): the first dim rows are position components, the last dim rows are momentum components.

Examples

>>> H, xs, xis, t, trajs = integrate_singularity(xi**2 + x**2, [x],
...                                              x0=0.0, xi0=3.0, tmax=6.0)

microlocal.interactive_symbol_analysis(pseudo_op, xlim=(-2, 2), ylim=(-2, 2), xi_range=(0.1, 5), eta_range=(-5, 5), density=50)

Launch an interactive dashboard for symbol exploration using ipywidgets.

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

Parameters:

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

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

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

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

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

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

Notes

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

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

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

• Supported visualization modes:

  • ‘Symbol Amplitude’ : |p(x,ξ)| or |p(x,y,ξ,η)|

  • ‘Symbol Phase’ : arg(p(x,ξ)) or similar in 2D

  • ‘Micro-Support (1/|p|)’ : Reciprocal of symbol magnitude

  • ‘Cotangent Fiber’ : Structure of symbol over frequency space at fixed x

  • ‘Characteristic Set’ : Zero set approximation {p ≈ 0}

  • ‘Characteristic Gradient’ : |∇p(x, ξ)| or |∇p(x₀, y₀, ξ, η)|

  • ‘Group Velocity Field’ : ∇_ξ p(x,ξ) or ∇_{ξ,η} p(x,y,ξ,η)

  • ‘Symplectic Vector Field’ : (∇_ξ p, -∇_x p) or similar in 2D

  • ‘Hamiltonian Flow’ : Trajectories generated by the Hamiltonian vector field

Raises:

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

• Prints –

• ------ –

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

microlocal.plot_bicharacteristics(symbol, initial_points, tspan, dim=None, projection='position', **kwargs)

Plot bicharacteristic curves (Hamiltonian flow trajectories).

For 1D: plots trajectories in the (x,ξ) phase plane. For 2D: projection can be ‘position’ (x-y), ‘frequency’ (ξ-η), or ‘mixed’ (x-ξ).

Parameters:

• symbol (sympy.Expr) – The principal symbol.

• initial_points (list of tuple) – List of initial conditions for the flow.

• tspan (tuple of float) – The (t_start, t_end) integration interval.

• dim (int, optional) – Dimension (1 or 2). Inferred from initial_points if None.

• projection (str, optional) – Projection plane for 2D plots. One of ‘position’, ‘frequency’, or ‘mixed’. Default is ‘position’.

• **kwargs (dict) – Additional arguments passed to bicharacteristic_flow.

Returns:

Displays a matplotlib plot.

Return type:

None

microlocal.plot_characteristic_set(symbol, x_range, xi_range, dim=None, resolution=200, **kwargs)

Plot the characteristic variety (zero set of the principal symbol).

For 1D: plots the contour p(x,ξ)=0 in the (x,ξ) phase plane. For 2D: plots a slice with fixed (ξ,η); you must provide xi0 and eta0 via kwargs.

Parameters:

• symbol (sympy.Expr) – The principal symbol.

• x_range (tuple of float) – Range for the first spatial variable (x).

• xi_range (tuple of float) – Range for the second spatial variable (y) in 2D, or frequency (ξ) in 1D.

• dim (int, optional) – Dimension (1 or 2). Inferred from symbol if None.

• resolution (int, optional) – Grid resolution for plotting. Default is 200.

• **kwargs (dict) – Additional arguments. For 2D, expects xi0 and eta0 (floats).

Returns:

Displays a matplotlib plot.

Return type:

None

microlocal.plot_hamiltonian_flow(op, x0=0.0, xi0=5.0, y0=0.0, eta0=0.0, tmax=1.0, n_steps=100, show_field=True)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

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

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

Notes

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

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

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

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

Raises:

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

• Displays –

• -------- –

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

microlocal.plot_matrix_1d(t, U, x, labels=None, quantity='real', save_path=None)

Plot each component of a matrix-valued 1D solution as a stacked space-time heatmap.

Parameters:

• t (ndarray, shape (n_times,)) – Time samples.

• U (ndarray, shape (n_times, size, Nx)) – Diagonal (or otherwise reduced) matrix solution components, one row of panels per index k = 0, …, size-1.

• x (ndarray, shape (Nx,)) – Spatial grid.

• labels (list of str, optional) – One label per component; defaults to [“u_1”, …, “u_size”].

• quantity ({'real', 'imag', 'abs'}, default='real') – Which part of U to plot.

• save_path (str, optional) – If given, the figure is saved to this path (dpi=150) before closing.

Returns:

The completed figure (already closed via _finish_headless).

Return type:

matplotlib.figure.Figure

microlocal.plot_matrix_field_1d(t, U, x, quantity='abs', component='diag', labels=None, save_path=None)

Space-time heatmap(s) for a matrix-valued 1D solution.

Parameters:

• t (ndarray, shape (n_times,)) – Time samples.

• U (ndarray, shape (n_times, N, N, Nx)) – Matrix-valued solution field, as returned by solve_matrix_field / solve_sylvester_field in 1D.

• x (ndarray, shape (Nx,)) – Spatial grid.

• quantity ({'real', 'imag', 'abs'}, default='abs') – Which part of the (reduced) field to plot.

• component ('diag' | 'trace' | 'frobenius' | (i, j), default='diag') – ‘diag’ – one panel per diagonal entry U_kk(x, t). ‘trace’ – single panel, sum_k U_kk(x, t). ‘frobenius’ – single panel, ||U(x, t)||_F. (i, j) – single panel, the (i, j) entry U_ij(x, t).

• labels (list of str, optional) – Panel labels used when component=’diag’; defaults to [“U_11”, “U_22”, …].

• save_path (str, optional) – If given, the figure is saved to this path (dpi=150) before closing.

Returns:

The completed figure (already closed via _finish_headless).

Return type:

matplotlib.figure.Figure

microlocal.plot_matrix_field_2d(t, U, x, y, times=None, quantity='abs', component='trace', save_path=None)

Plot snapshot panels for a matrix-valued 2D solution field.

Parameters:

• t (numpy.ndarray) – 1D array of time samples, shape (n_times,).

• U (numpy.ndarray) – Matrix-valued solution field, shape (n_times, N, N, Nx, Ny).

• x (numpy.ndarray) – 1D arrays of spatial grid coordinates.

• y (numpy.ndarray) – 1D arrays of spatial grid coordinates.

• times (array_like of int, optional) – Indices into t selecting which snapshots to plot. Defaults to 6 evenly spaced indices.

• quantity ({'real', 'imag', 'abs'}, optional) – Which part of the reduced field to plot. Default is ‘abs’.

• component ({'trace', 'frobenius'} or tuple of int, optional) – How to reduce the matrix field to a scalar field. ‘trace’ : sum of diagonal entries. ‘frobenius’ : Frobenius norm. (i, j) : specific matrix entry. Default is ‘trace’.

• save_path (str, optional) – If provided, saves the figure to this path.

Returns:

The completed and closed figure object.

Return type:

matplotlib.figure.Figure

microlocal.plot_pseudospectrum(Lambda, resolvent_norm, sigma_min_grid, epsilon_levels, eigenvalues)

Plot pseudospectrum results.

Parameters:

• Lambda (ndarray) – Complex λ grid

• resolvent_norm (ndarray) – Resolvent norms

• sigma_min_grid (ndarray) – Smallest singular values

• epsilon_levels (list) – Contour levels

• eigenvalues (ndarray or None) – Eigenvalues to overlay

microlocal.plot_scalar_1d(t, U, x, title='u(x, t)', quantity='real', n_snapshots=6, save_path=None)

Plot a scalar 1D space-time solution as a combined heatmap and snapshot overlay.

Parameters:

• t (ndarray, shape (n_times,)) – Time samples.

• U (ndarray, shape (n_times, Nx)) – Solution values u(x, t) sampled on the grid.

• x (ndarray, shape (Nx,)) – Spatial grid.

• title (str, default="u(x, t)") – Base title used for the heatmap panel.

• quantity ({'real', 'imag', 'abs'}, default='real') – Which part of U to plot.

• n_snapshots (int, default=6) – Number of time slices drawn as line overlays in the second panel.

• save_path (str, optional) – If given, the figure is saved to this path (dpi=150) before closing.

Returns:

The completed figure (already closed via _finish_headless, so it will not display inline; use save_path or re-show it explicitly).

Return type:

matplotlib.figure.Figure

microlocal.plot_scalar_2d(t, U, x, y, times=None, quantity='real', save_path=None)

Plot a scalar 2D solution at selected time instants as a row of side-by-side pcolormesh panels.

Parameters:

• t (ndarray, shape (n_times,)) – Time samples.

• U (ndarray, shape (n_times, Nx, Ny)) – Solution values u(x, y, t) sampled on the grid.

• x (ndarray) – Spatial grids along each axis.

• y (ndarray) – Spatial grids along each axis.

• times (array_like of int, optional) – Indices into t selecting which snapshots to plot. Defaults to 6 indices evenly spaced across the whole time range.

• quantity ({'real', 'imag', 'abs'}, default='real') – Which part of U to plot.

• save_path (str, optional) – If given, the figure is saved to this path (dpi=150) before closing.

Returns:

The completed figure (already closed via _finish_headless).

Return type:

matplotlib.figure.Figure

microlocal.plot_symplectic_vector_field(op, xlim=(-2, 2), klim=(-5, 5), density=30)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

Raises:

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

Notes

• Only supports one-dimensional operators.

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

• Numerical evaluation is done via lambdify with NumPy backend.

• Visualization uses matplotlib quiver plot to show vector directions.

microlocal.plot_wave_solution_1d(t, U, V, x, quantity='real', save_path=None)

Plot side-by-side space-time heatmaps of displacement and velocity.

Intended for the output of solve_second_order (scalar, 1D case).

Parameters:

• t (numpy.ndarray) – 1D array of time samples, shape (n_times,).

• U (numpy.ndarray) – Displacement field u(x, t), shape (n_times, Nx).

• V (numpy.ndarray) – Velocity field ∂_t u(x, t), shape (n_times, Nx).

• x (numpy.ndarray) – 1D array of spatial grid coordinates, shape (Nx,).

• quantity ({'real', 'imag', 'abs'}, optional) – Which part of the fields to plot. Default is ‘real’.

• save_path (str, optional) – If provided, saves the figure to this path.

Returns:

The completed and closed figure object.

Return type:

matplotlib.figure.Figure

microlocal.plot_wavefront_set(symbol, initial_sing_support, tspan, dim=None, projection='cotangent', n_steps=500, cmap='plasma', show_flow=True, show_endpoints=True, title=None, ax=None)

Plot the wavefront set WF(u) of a distribution u whose singularities propagate along bicharacteristics of the operator with symbol p.

The wavefront set is represented as a subset of the cotangent bundle T*ℝⁿ. Each initial point (x₀, ξ₀) seeds a bicharacteristic strip; the union of these strips in phase space approximates WF(u) at time tspan[1].

Parameters:

• symbol (sympy expression) – Principal symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D.

• initial_sing_support (list of tuples) – Seed points on the wavefront set at t=0. 1D: list of (x₀, ξ₀) pairs. 2D: list of (x₀, y₀, ξ₀, η₀) quadruples.

• tspan (tuple) – (t_start, t_end) for bicharacteristic integration.

• dim (int, optional) – Dimension (1 or 2). Inferred from seed points if None.

• projection (str) –

  Which subspace to visualise (dimension-dependent): - 1D:

  ’cotangent’ – (x, ξ) phase-space portrait [default] ‘position’ – x(t) projected onto ℝ (singular support)

  • 2D:

    ’cotangent’ – (x, ξ) slice (ignoring y, η) ‘position’ – (x, y) projection (singular support in ℝ²) ‘frequency’ – (ξ, η) projection (directions of non-smoothness) ‘full’ – 2×2 grid: position, frequency, (x,ξ), (y,η)

• n_steps (int) – Number of integration steps per bicharacteristic.

• cmap (str) – Matplotlib colormap used to colour individual bicharacteristics (colour encodes the index of the seed point, i.e. which part of the initial singular support it originated from).

• show_flow (bool) – If True, draw the full bicharacteristic strip (trajectory in phase space). If False, only show the endpoint scatter.

• show_endpoints (bool) – If True, mark the initial point (green •) and final point (red •) of each bicharacteristic.

• title (str, optional) – Figure title. A sensible default is generated if None.

• ax (matplotlib Axes or array of Axes, optional) – Axes to draw into. If None a new figure is created. For projection=’full’ (2D) pass an array of 4 Axes or leave None.

Returns:

• fig (matplotlib Figure)

• axes (Axes or array of Axes)

Examples

1D example – Schrödinger-type operator, horizontal line singularity:

x, xi = symbols('x xi', real=True)
p = xi**2 - (1 - x**2)           # simple potential well symbol
seeds = [(xi_val, float(xi_val)) for xi_val in np.linspace(-1, 1, 12)]
fig, ax = plot_wavefront_set(p, seeds, tspan=(0, 4), dim=1)
plt.show()

2D example – wave operator, outward circular wavefront:

x, y, xi, eta = symbols('x y xi eta', real=True)
p = xi**2 + eta**2 - 1
seeds = [(np.cos(t), np.sin(t), np.cos(t), np.sin(t))
         for t in np.linspace(0, 2*np.pi, 24, endpoint=False)]
fig, axes = plot_wavefront_set(p, seeds, tspan=(0, 2), dim=2,
                               projection='full')
plt.show()

microlocal.propagate_singularity(symbol, initial_sing_support, tspan, dim=None, n_samples=None)

Propagate singular support along bicharacteristics.

Parameters:

• symbol (sympy.Expr) – The principal symbol.

• initial_sing_support (list of tuple) – List of initial phase-space points (x₀, ξ₀) for 1D or (x₀, y₀, ξ₀, η₀) for 2D.

• tspan (tuple of float) – The (t_start, t_end) integration interval.

• dim (int, optional) – Spatial dimension (1 or 2). Inferred from symbol if None.

• n_samples (int, optional) – Currently unused. Kept for API compatibility.

Returns:

A dictionary containing: - ‘trajectories’ : list of dict (Bicharacteristic flow results) - ‘endpoints’ : list of tuple (Final phase-space coordinates) - ‘initial’ : list of tuple (The initial seed points provided)

Return type:

dict

microlocal.visualize_characteristic_gradient(op, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

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

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

Returns:

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

Return type:

None

Notes

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

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

• Numerical differentiation is performed using np.gradient.

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

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

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

Notes

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

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

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

Raises:

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

• Displays –

• ------ –

• A matplotlib contour plot showing either: –

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

microlocal.visualize_fiber(op, x_grid, xi_grid, x0=0.0, y0=0.0)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

Notes

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

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

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

Raises:

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

microlocal.visualize_micro_support(op, xlim=(-2, 2), klim=(-10, 10), threshold=0.001, density=300, xi0=0.0, eta0=0.0)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

Raises:

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

Notes

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

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

• The resulting plot helps identify characteristic sets.

microlocal.visualize_phase(op, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

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

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

• Notes

• 1D (- In)

• 2D (- In)

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

• Raises

• NotImplementedError (-)

microlocal.visualize_symbol_amplitude(op, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0)

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

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

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator whose symbol is to be analyzed. Must expose attributes: dim, p_func, vars_x, symbol, and expr.

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

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

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

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

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

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

Notes

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

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

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

microlocal.visualize_wavefront_set(op, seeds=None, tspan=(0, 3.0), projection='cotangent', n_steps=500, cmap='plasma', show_flow=True, show_endpoints=True, title=None, x0=0.0, y0=0.0, xi0=1.0, eta0=0.0, spread=2.0, n_seeds=25, radius=0.15)

Visualize the wavefront set WF(u) obtained by propagating seed singularities along the bicharacteristics of the operator’s symbol.

If seeds is not given, a default fan/point-source is built from (x0, y0) and (xi0, eta0).

Parameters:

• op (PseudoDifferentialOperator) – The pseudo-differential operator.

• seeds (list of tuple, optional) – Seed points for the wavefront set. If None, auto-generated.

• tspan (tuple of float, optional) – Integration time interval (t_start, t_end). Default is (0, 3.0).

• projection (str, optional) – Subspace to visualize. 1D: ‘cotangent’, ‘position’. 2D: ‘cotangent’, ‘position’, ‘frequency’, ‘mixed_x’, ‘mixed_y’, ‘full’.

• n_steps (int, optional) – Number of integration steps per bicharacteristic. Default is 500.

• cmap (str, optional) – Matplotlib colormap. Default is ‘plasma’.

• show_flow (bool, optional) – Whether to draw the full bicharacteristic strip. Default is True.

• show_endpoints (bool, optional) – Whether to mark initial and final points. Default is True.

• title (str, optional) – Figure title.

• x0 (float, optional) – Base spatial coordinates for auto-generated seeds.

• y0 (float, optional) – Base spatial coordinates for auto-generated seeds.

• xi0 (float, optional) – Base frequency coordinates for auto-generated seeds.

• eta0 (float, optional) – Base frequency coordinates for auto-generated seeds.

• spread (float, optional) – Frequency spread for 1D auto-generated seeds.

• n_seeds (int, optional) – Number of seed points to generate if seeds is None.

• radius (float, optional) – Radius for 2D point-source auto-generated seeds.

Returns:

• fig (matplotlib.figure.Figure) – The generated figure.

• axes (matplotlib.axes.Axes or ndarray of Axes) – The axes object(s).