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.

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.

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

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

References

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

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, method='fast')[source]

Bohr–Sommerfeld quantisation for 1D bound states.

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

Parameters:
  • H (sympy expression) – Hamiltonian H(x,p).

  • n_max (int) – Maximum quantum number.

  • x_range (tuple) – Spatial range for turning points.

  • hbar (float) – Planck constant.

  • method (str) – Ignored, kept for compatibility.

Returns:

Quantised energies and actions.

Return type:

dict

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

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:

Contains: - ‘implicit’: the symbol expression - ‘equation’: sympy Eq(symbol,0) - ‘explicit’: list of explicit solutions ξ(x) (1D only) or None - ‘function’: a callable that evaluates the symbol

Return type:

dict

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

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, p=None)[source]

Compute the Maslov index for a single trajectory. traj : dict returned by bicharacteristic_flow (must contain J11..J22 or J) p : unused (kept for compatibility)

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

Find caustics (envelope of bicharacteristics) in 1D.

Simplified condition: where d²p/dξ² ≈ 0 (turning points in frequency).

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

Plot bicharacteristic curves.

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

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

Plot the characteristic variety.

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

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)[source]

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 = sp.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 = sp.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)[source]

Propagate singular support along bicharacteristics.

microlocal.test_bicharacteristic_flow()[source]
microlocal.test_bohr_sommerfeld()[source]
microlocal.test_characteristic_variety()[source]
microlocal.test_wkb()[source]