asymptotic

asymptotic.py — Large-parameter asymptotics for oscillatory and Laplace-type integrals

Overview

This module provides symbolic-numerical tools for computing the asymptotic behaviour of parameter-dependent integrals of the form

I(λ) = ∫ a(x) exp(iλφ(x)) dx, λ → +∞

as the large parameter λ grows without bound. The nature of the phase function φ determines which asymptotic method applies:

φ	Integral form	Method
real	oscillatory	Stationary phase
purely imaginary	exponentially damped	Laplace
genuinely complex	oscillatory + damped	Saddle-point (method of steepest descent)

The correct method is selected automatically from the symbolic expression of φ when the analyzer is initialised (method=IntegralMethod.AUTO, the default). It can also be set explicitly.

Mathematical background

All three methods share the same underlying idea: the dominant contribution to I(λ) as λ → ∞ comes from a small neighbourhood of a critical point x_c where ∇φ(x_c) = 0. Away from x_c, rapid oscillations (or exponential decay) make the integrand self-cancelling.

Stationary phase (φ real)

The leading term at a non-degenerate critical point (det ∇²φ ≠ 0, Morse point) is

I(λ) ≈ (2π/λ)^(n/2) · a(x_c) · exp(iλφ(x_c)) · exp(iπμ/4) / √|det ∇²φ(x_c)|

where n is the dimension and μ = n − 2σ is the Maslov index (σ = number of negative eigenvalues of ∇²φ). Degenerate critical points require special treatment: corank-1 singularities with a non-zero cubic term yield Airy integrals (decay O(λ^(−1/3)) in 1D, O(λ^(−5/6)) in 2D); those with a vanishing cubic but non-zero quartic term yield Pearcey integrals (decay O(λ^(−3/4))).

Laplace’s method (φ = iψ, ψ real)

The integrand concentrates exponentially around the minimum x_c of ψ. The leading term is identical to the stationary-phase formula with the oscillatory factor replaced by a real Gaussian:

I(λ) ≈ (2π/λ)^(n/2) · a(x_c) · exp(−λψ(x_c)) / √|det ∇²ψ(x_c)|

A second-order correction O(λ^(−n/2−1)) involving amplitude derivatives and phase anharmonicity (cubic/quartic tensors) is also computed.

Saddle-point / steepest descent (φ complex)

The integration contour is deformed into ℂⁿ to pass through a saddle point z_c ∈ ℂⁿ satisfying ∇φ(z_c) = 0. On the steepest-descent contour through z_c the phase Im(λφ) is stationary and Re(λφ) grows as fast as possible, making the integrand a complex Gaussian. The asymptotic formula is formally identical to the Morse case:

I(λ) ≈ (2π/λ)^(n/2) · a(z_c) · exp(iλφ(z_c)) / √det ∇²φ(z_c)

where the square root is taken on the principal branch. Limitation: this implementation uses a naive continuation strategy (minimising |∇φ(z)|² over ℝ^(2n)) and does NOT verify that the original contour can be deformed through the found saddle (Picard- Lefschetz theory). A RuntimeWarning is always emitted.

References

class asymptotic.Analyzer(phase_expr, amplitude_expr, variables, domain=None, tolerance=1e-06, cubic_threshold=None, method: IntegralMethod = IntegralMethod.AUTO)[source]

Bases: object

Handles symbolic analysis of phase and amplitude functions.

Supports three integration paradigms selected via method:

IntegralMethod.STATIONARY_PHASE:
Oscillatory integrals I(λ) = ∫ a(x) exp(iλφ(x)) dx, φ real. All singularity types (Morse, Airy, Pearcey) are supported.
IntegralMethod.LAPLACE:
Exponentially damped integrals I(λ) = ∫ a(x) exp(-λφ(x)) dx, φ real. Only non-degenerate minima of φ are relevant.
IntegralMethod.SADDLE_POINT:
Mixed integrals I(λ) = ∫ a(x) exp(iλφ(x)) dx, φ genuinely complex. Saddle points are searched in ℂⁿ by continuation from real guesses.
IntegralMethod.AUTO (default):
The analyzer inspects φ symbolically and chooses automatically among the three concrete methods. The resolved method is written back to self.method after __init__ completes.

The main workflow is:

Initialize with symbolic SymPy expressions (method=AUTO by default).
Find critical / saddle points using find_critical_points().
Analyze each point using analyze_point().
Use AsymptoticEvaluator (unified façade) to compute contributions.

phase_expr: SymPy expression for the phase function φ(x).

amplitude_expr: SymPy expression for the amplitude function a(x).

variables: List of SymPy symbols representing the integration variables.

dim

Dimension of the integration domain.

Type:: int

domain

Optional bounds [(min, max), …] for each variable.

Type:: Optional[List[Tuple]]

tolerance

Numerical tolerance for detecting zeros and critical points.

Type:: float

cubic_threshold

Absolute threshold for distinguishing Airy vs Pearcey singularities. Default: max(1e-5, 10*tolerance). Unused for LAPLACE.

Type:: float

method

Resolved integration method (never AUTO after __init__).

Type:: IntegralMethod

analyze_point(xc) → CriticalPoint[source]

Perform complete analysis of a critical point (real or complex).

For STATIONARY_PHASE and LAPLACE the coordinates xc are real (numpy float array). For SADDLE_POINT, xc may be complex (numpy complex array produced by SaddlePointEvaluator.find_saddle_points).

Computes all geometric and analytical properties needed for asymptotic evaluation:

Phase value φ(x_c) and amplitude value a(x_c)
Hessian matrix ∇²φ and its properties (determinant, eigenvalues, signature)
Higher-order derivatives: D3 and D4 tensors of φ, gradients and Hessians of a
Classification of singularity type (Morse, Airy, Pearcey, etc.)
Canonical coefficients for degenerate cases
Hessian inverse (for Morse points only)

Parameters:: xc – Coordinates of the critical point. Real numpy array for STATIONARY_PHASE / LAPLACE; complex numpy array for SADDLE_POINT.
Returns:: CriticalPoint object containing all computed properties.

find_critical_points(initial_guesses=None) → List[ndarray][source]

Locate critical points where ∇φ(x) = 0.

Delegates to caustics.find_critical_points_numerical, the shared numerical kernel (L-BFGS-B minimisation of |∇φ|² + DBSCAN dedup).

Parameters:: initial_guesses – List of starting coordinate arrays for optimization. If None, uses [0, …] and domain center (if domain is specified). Provide multiple guesses to search for multiple critical points.
Returns:: List of unique critical point coordinates (as numpy arrays) found within the specified tolerance. Empty list if no critical points are found.

class asymptotic.AsymptoticContribution(leading_term: complex, correction_term: complex, total_value: complex, point: CriticalPoint, order_leading: float, method: IntegralMethod = IntegralMethod.STATIONARY_PHASE)[source]

Bases: object

Represents the calculated asymptotic contribution from a specific critical point.

The total asymptotic expansion typically has the form:: I(λ) ≈ leading_term + correction_term + O(λ^(-order_leading - 2))

For Morse points, the correction term is of order λ^(-n/2-1) relative to λ^(-n/2). For degenerate singularities (Airy, Pearcey), typically only the leading term is computed, as correction terms require more sophisticated analysis.

leading_term

The dominant term, O(λ^(-order_leading)): for Morse in 2D: O(λ^(-1)), for Airy 1D: O(λ^(-1/3)), for Airy 2D: O(λ^(-5/6)), for Pearcey: O(λ^(-3/4))

Type:: complex

correction_term

The next-order correction term: for Morse: O(λ^(-n/2-1)), for degenerate cases: typically 0j (not computed)

Type:: complex

total_value

Sum of leading_term + correction_term.

Type:: complex

point

The source critical point for this contribution.

Type:: CriticalPoint

order_leading

The exponent p in the scaling λ^(-p) of the leading term.

Type:: float

method

The asymptotic method used to compute this contribution (STATIONARY_PHASE, LAPLACE or SADDLE_POINT).

Type:: IntegralMethod

correction_term: complex

leading_term: complex

method: IntegralMethod = 'stationary_phase'

order_leading: float

point: CriticalPoint

total_value: complex

class asymptotic.AsymptoticEvaluator(tolerance: float = 1e-08)[source]

Bases: object

Unified façade that dispatches asymptotic evaluation to the appropriate evaluator based on the integration method stored in a CriticalPoint.

This is the recommended entry point for end users. Routing table:

cp.method == STATIONARY_PHASE → StationaryPhaseEvaluator
cp.method == LAPLACE → LaplaceEvaluator
cp.method == SADDLE_POINT → SaddlePointEvaluator

The method is determined automatically when the analyzer is constructed with method=AUTO (the default).

Usage

>>> analyzer = Analyzer(phi, amp, [x, y])  # AUTO by default
>>> print(analyzer.method)   # e.g. IntegralMethod.SADDLE_POINT
>>> pts = analyzer.find_critical_points([np.array([0., 0.])])
>>> cp  = analyzer.analyze_point(pts[0])
>>> result = AsymptoticEvaluator().evaluate(cp, lam=100)
>>> print(result.method, result.total_value)

For SADDLE_POINT, use SaddlePointEvaluator.find_saddle_points() to obtain complex coordinates before calling analyze_point():

>>> sp_eval = SaddlePointEvaluator()
>>> saddles = sp_eval.find_saddle_points(analyzer, real_guesses)
>>> cp = analyzer.analyze_point(saddles[0])
>>> result = AsymptoticEvaluator().evaluate(cp, lam=100)

sp_evaluator

Type:: StationaryPhaseEvaluator

laplace_evaluator

Type:: LaplaceEvaluator

saddle_evaluator

Type:: SaddlePointEvaluator

evaluate(cp: CriticalPoint, lam: float) → AsymptoticContribution[source]

Evaluate the asymptotic contribution at parameter λ.

The evaluation method is selected from cp.method:

STATIONARY_PHASE Full singularity-type dispatch (Morse, Airy, Pearcey …)
LAPLACE Laplace formula with O(1/λ) corrections
SADDLE_POINT Complex Morse formula (naive; see SaddlePointEvaluator)

Parameters:

cp (CriticalPoint) – Critical/saddle point from Analyzer.analyze_point().
lam (float) – Large asymptotic parameter λ > 0.

Return type:

AsymptoticContribution

Raises:

ValueError – If cp.method is None, AUTO, or unrecognised.

class asymptotic.AsymptoticVisualizer(analyzer: Analyzer)[source]

Bases: object

Visualization toolkit for asymptotic analysis — supports all three integration methods (STATIONARY_PHASE, LAPLACE, SADDLE_POINT).

Provides three diagnostic plots, each adapted to the nature of the phase φ:

plot_phase_landscape 2D contour map of φ with critical/saddle-point overlay.
- φ real (STATIONARY_PHASE): single panel, Re(φ).
- φ imag (LAPLACE): single panel, Im(φ) = ψ (the damping potential).
- φ complex (SADDLE_POINT): two panels side-by-side, Re(φ) and Im(φ).
The display parameter overrides the automatic choice ('real', 'imag', 'both', 'abs', 'arg').
plot_integrand 2D map of the integrand f(x,y) = a(x,y)·exp(iλφ(x,y)) at a given λ.
- STATIONARY_PHASE: single panel, Re(f) — pure oscillation, |f| = const.
- LAPLACE: single panel, Re(f) = a·exp(-λψ) — exponential envelope.
- SADDLE_POINT: two panels, Re(f) and |f| = |a|·exp(-λ Im φ), revealing both the oscillation pattern and the exponential damping.
plot_asymptotic_convergence Log-log plot of |I₀(λ)| and |I₁(λ)| vs λ for any dimension and any method. Overlays the theoretical decay slope λ^(-p) for verification.

Notes

plot_phase_landscape and plot_integrand require dim = 2.
plot_asymptotic_convergence works for any dimension.

plot_asymptotic_convergence(cp: CriticalPoint, lambda_start: float = 10, lambda_end: float = 1000, num_points: int = 50) → None[source]

Log-log convergence diagnostic for any method and any dimension.

Plots |I₀(λ)| and |I₁(λ)| vs λ and compares the empirical slope with the theoretical decay exponent −p:

Method / type	Theoretical slope −p
Morse (any dim n)	−n/2
Airy 1D	−1/3
Airy 2D	−5/6
Pearcey	−3/4
Laplace (any n)	−n/2
Saddle-point	−n/2 (complex Morse)

Parameters:

cp (CriticalPoint) – Critical / saddle point (any method, any dimension).
lambda_start (float) – Minimum λ for the convergence sweep (default 10).
lambda_end (float) – Maximum λ (default 1000).
num_points (int) – Number of log-spaced λ samples (default 50).

Notes

A straight line on the log-log plot confirms the asymptotic regime.
Deviations at small λ indicate pre-asymptotic behaviour.
The correction term |I₁| should be parallel to |I₀| but shifted down by slope −1 for Morse / Laplace points.

plot_integrand(lam_value: float, bounds: Tuple[Tuple[float, float], Tuple[float, float]] = ((-3, 3), (-3, 3)), points_per_axis: int = 200) → None[source]

Visualize the structure of the integrand f = a(x,y)·exp(iλφ(x,y)) at a given parameter λ.

The panels shown depend on the integration method:

Method	Panels
`STATIONARY_PHASE`	Re(f) — oscillation with unit-modulus envelope
`LAPLACE`	f = a·exp(-λψ) — real exponential concentration
`SADDLE_POINT`	Re(f) \| \|f\| = \|a\|·exp(-λ Im φ) side-by-side: left shows oscillations, right shows the exponential damping envelope

Parameters:

lam_value (float) – Parameter λ. Larger values produce finer oscillations / sharper concentration.
bounds (pair of (min, max) pairs) – Spatial domain ((x_min, x_max), (y_min, y_max)).
points_per_axis (int) – Grid resolution (default 200). Increase for large λ to resolve fine oscillations.

Notes

For STATIONARY_PHASE, as λ increases the oscillations become finer everywhere except near stationary points (∇φ = 0), where the phase is locally flat — this is the geometric core of the method.
For LAPLACE, the integrand concentrates sharply around the minimum of Im(φ) = ψ, illustrating why only a small neighbourhood contributes.
For SADDLE_POINT, |f| reveals the exponential ridge structure while Re(f) shows the additional rapid oscillations along the ridge.

plot_phase_landscape(critical_points: List[CriticalPoint], bounds: Tuple[Tuple[float, float], Tuple[float, float]] = ((-3, 3), (-3, 3)), points_per_axis: int = 120, display: str | None = None) → None[source]

Visualize the phase function φ(x,y) with critical/saddle-point overlay.

The panels shown depend on the integration method (or on display):

Method / display	Panels
`STATIONARY_PHASE`	Re(φ)
`display='real'`
`LAPLACE`	Im(φ) = ψ (the damping potential)
`display='imag'`
`SADDLE_POINT`	Re(φ) \| Im(φ) side-by-side
`display='both'`
`display='abs'`	\|φ\|
`display='arg'`	arg(φ)

Parameters:

critical_points (list of CriticalPoint) – Points to overlay (real or complex coordinates accepted).
bounds (pair of (min, max) pairs) – Spatial domain ((x_min, x_max), (y_min, y_max)).
points_per_axis (int) – Grid resolution (default 120).
display (str or None) – Override automatic panel selection. One of: 'real', 'imag', 'both', 'abs', 'arg'.

Notes

Marker conventions (shared with plot_integrand):
- ○ red — Morse (non-degenerate)
- ★ orange — Airy (corank 1, cubic)
- ◆ magenta — Pearcey (corank 1, quartic)
- □ gray — Higher-order / unclassified

class asymptotic.CriticalPoint(position: ndarray, phase_value: complex, amplitude_value: complex, singularity_type: SingularityType, hessian_matrix: ndarray, hessian_inv: ndarray | None = None, hessian_det: float = 0.0, signature: int = 0, eigenvalues: ndarray = <factory>, eigenvectors: ndarray = <factory>, grad_amp: ndarray | None = None, hess_amp: ndarray | None = None, phase_d3: ndarray | None = None, phase_d4: ndarray | None = None, canonical_coefficients: Dict | None = None, method: IntegralMethod = None)[source]

Bases: object

Stores all geometric and analytical properties of a critical point.

A critical point x_c satisfies ∇φ(x_c) = 0. This class contains all the data needed to compute its asymptotic contribution to the stationary phase integral ∫ a(x) exp(iλφ(x)) dx as λ → ∞.

position

Coordinates of the critical point x_c.

Type:: np.ndarray

phase_value

Value of phase function φ(x_c).

Type:: complex

amplitude_value

Value of amplitude function a(x_c).

Type:: complex

singularity_type

Classification determining which formula to use.

Type:: SingularityType

hessian_matrix

The Hessian matrix ∇²φ at x_c.

Type:: np.ndarray

hessian_inv

Inverse of the Hessian (for Morse points only).

Type:: Optional[np.ndarray]

hessian_det

Determinant of the Hessian.

Type:: float

signature

Number of negative eigenvalues of the Hessian (Morse index).

Type:: int

eigenvalues

Eigenvalues of the Hessian matrix.

Type:: np.ndarray

eigenvectors

Eigenvectors of the Hessian matrix.

Type:: np.ndarray

grad_amp

Gradient of amplitude ∇a at x_c.

Type:: Optional[np.ndarray]

hess_amp

Hessian of amplitude ∇²a at x_c.

Type:: Optional[np.ndarray]

phase_d3

Rank-3 tensor of 3rd derivatives of φ.

Type:: Optional[np.ndarray]

phase_d4

Rank-4 tensor of 4th derivatives of φ.

Type:: Optional[np.ndarray]

canonical_coefficients

Coefficients for normal forms (Airy/Pearcey canonical representations). Contains keys like ‘cubic’, ‘quartic’, ‘quadratic_transverse’ depending on singularity type.

Type:: Optional[Dict]

amplitude_value: complex

canonical_coefficients: Dict | None = None

eigenvalues: ndarray

eigenvectors: ndarray

grad_amp: ndarray | None = None

hess_amp: ndarray | None = None

hessian_det: float = 0.0

hessian_inv: ndarray | None = None

hessian_matrix: ndarray

method: IntegralMethod = None

phase_d3: ndarray | None = None

phase_d4: ndarray | None = None

phase_value: complex

position: ndarray

signature: int = 0

singularity_type: SingularityType

class asymptotic.IntegralMethod(*values)[source]

Bases: Enum

Selects which asymptotic method to apply to the integral I(λ).

The three concrete methods correspond to the three possible natures of the phase function φ(x) appearing in the exponential exp(iλφ(x)):

STATIONARY_PHASE:
φ is purely real. I(λ) = ∫ a(x) exp(iλφ(x)) dx, φ ∈ ℝ, λ → +∞ The integrand oscillates with unit modulus; contributions arise from stationary points ∇φ = 0. Decay rate: O(λ^(-n/2)). Typical applications: wave optics, quantum mechanics, Fourier integrals.
LAPLACE:
φ is purely imaginary, i.e. φ(x) = i·ψ(x) with ψ ∈ ℝ. I(λ) = ∫ a(x) exp(iλ·iψ(x)) dx = ∫ a(x) exp(-λψ(x)) dx, λ → +∞ The integrand is real and exponentially damped; contributions arise from minima of ψ where ∇ψ = 0 and ∇²ψ > 0. Typical applications: large deviations, Bayesian inference, statistical mechanics partition functions.
SADDLE_POINT:
φ is genuinely complex, φ = φ_R + i·φ_I with both φ_R ≠ 0 and φ_I ≠ 0. I(λ) = ∫ a(x) exp(iλφ_R(x)) exp(-λφ_I(x)) dx, λ → +∞ The integrand both oscillates and is exponentially modulated. Contributions come from saddle points in ℂⁿ found by analytically continuing ∇φ(z) = 0 into the complex plane. The integration contour must be deformed to pass through these saddle points along the steepest-descent direction. Note: this implementation uses a naive continuation strategy; see SaddlePointEvaluator for limitations.
AUTO:
automatic detection (default). The analyzer inspects φ symbolically (via sympy.im / sympy.re) and falls back to a numerical test if the symbolic check is inconclusive. The detected method is stored back in Analyzer.method after __init__ so the user can always query which method was chosen.

Hierarchy

SADDLE_POINT is the general case; the other two are special cases:
SADDLE_POINT with φ_I ≡ 0 → STATIONARY_PHASE
SADDLE_POINT with φ_R ≡ 0 → LAPLACE

AUTO = 'auto'

LAPLACE = 'laplace'

SADDLE_POINT = 'saddle_point'

STATIONARY_PHASE = 'stationary_phase'

class asymptotic.LaplaceEvaluator[source]

Bases: object

Evaluator for exponentially damped integrals of the form:: I(λ) = ∫ a(x) exp(-λ φ(x)) dx, λ → +∞

Uses Laplace’s method with second-order asymptotic corrections O(λ^(-n/2-1)).

The critical point must be a strict minimum of φ (positive definite Hessian). Saddle points and maxima are not supported: the Laplace method relies on the Gaussian concentration of the integrand around the minimum.

Returns an AsymptoticContribution with method=IntegralMethod.LAPLACE for consistency with StationaryPhaseEvaluator.

evaluate(cp: CriticalPoint, lam: float) → AsymptoticContribution[source]

Standard Laplace formula for n dimensions with anharmonic corrections.

Leading term (Order 0):: I₀(λ) = a(x_c) · exp(-λ φ(x_c)) · (2π/λ)^(n/2) · |det H|^(-1/2)
Correction term (Order 1, relative order O(1/λ)):: I₁(λ) = I₀(λ) · (1/λ) · C
where the real correction factor C is:: C = (1/2) Tr(H⁻¹ ∇²a) − (1/2) ⟨∇a, (H⁻¹ ⊗ H⁻¹) D³φ⟩ − (1/8) (H⁻¹ ⊗ H⁻¹) : D⁴φ + (5/24) (H⁻¹ ⊗ H⁻¹ ⊗ H⁻¹) : (D³φ ⊗ D³φ)

Parameters:

cp – CriticalPoint with a non-degenerate, positive definite Hessian. Must contain hessian_inv (or hessian_matrix), hess_amp, grad_amp, phase_d3, phase_d4.
lam – Large parameter λ. The approximation improves as λ → +∞.

Returns:

An object with the following attributes:

leading_term (float): I₀(λ), real for real a and φ.
correction_term (float): I₁(λ), the O(λ^(-n/2-1)) correction.
total_value (float): I₀ + I₁.
order_leading (float): n/2.
method (IntegralMethod): IntegralMethod.LAPLACE.

Return type:

AsymptoticContribution

Raises:

ValueError – If the Hessian is singular (det H ≈ 0).

class asymptotic.SaddlePointEvaluator(tolerance: float = 1e-08)[source]

Bases: object

Naive saddle-point evaluator for integrals with a genuinely complex phase.

Handles integrals of the form:: I(λ) = ∫ a(x) exp(iλφ(x)) dx, φ = φ_R + i·φ_I, λ → +∞

where both the real part φ_R and the imaginary part φ_I are non-trivial. The integrand simultaneously oscillates (φ_R) and is exponentially damped (φ_I). The asymptotic contribution is dominated by saddle points in ℂⁿ, i.e. solutions of ∇φ(z) = 0 with z ∈ ℂⁿ.

Strategy (naive continuation)

Start from real initial guesses x₀ ∈ ℝⁿ (supplied by the caller).
Analytically continue into ℂⁿ by minimising |∇φ(z)|² over the 2n real degrees of freedom (Re z, Im z), using scipy.optimize.minimize.
Accept a point z_c if |∇φ(z_c)|² < tolerance.
Apply the standard Morse formula with the complex Hessian:

I(λ) ≈ (2π/λ)^(n/2) · a(z_c) · exp(iλφ(z_c)) · 1/√det(∇²φ(z_c))

where the complex square root is chosen with positive real part (principal branch convention).

Limitations and warnings

Contour validity is NOT checked. The naive continuation finds a saddle point algebraically but does NOT verify that the original real integration contour can be deformed through that saddle without crossing other singularities (Picard-Lefschetz theory). A RuntimeWarning is always emitted to remind the user of this.
Branch choice. The complex square root √det H is multi-valued. This implementation uses numpy’s principal branch (argument in (-π, π]). The correct branch depends on the global topology of the steepest-descent contour.
Multiple saddles. When several saddle points are found, ALL contributions are returned; their relative signs (Stokes phenomena) are not resolved.
Degenerate saddles (det H ≈ 0) are not supported; a warning is issued and a zero contribution is returned.

References

evaluate(cp: CriticalPoint, lam: float) → AsymptoticContribution[source]

Evaluate the leading-order saddle-point contribution.

Uses the standard multi-dimensional Morse formula extended to a complex critical point z_c:

I(λ) ≈ (2π/λ)^(n/2) · a(z_c) · exp(iλφ(z_c)) · 1/√det(∇²φ(z_c))

The complex determinant det(∇²φ(z_c)) is evaluated with numpy’s principal square root.

Warning

This contribution is valid only if the original integration contour can be deformed through z_c along a steepest-descent path. This is NOT verified here (Picard-Lefschetz theory). Always examine the result critically.

Parameters:

cp (CriticalPoint) – Saddle point with method=SADDLE_POINT, produced by Analyzer.analyze_point() called with a complex coordinate returned by find_saddle_points().
lam (float) – Large parameter λ > 0.

Returns:

leading_term : saddle-point formula value.
correction_term : 0j (not implemented for saddle points).
order_leading : n/2.
method : IntegralMethod.SADDLE_POINT.

Return type:

AsymptoticContribution

find_saddle_points(analyzer: Analyzer, initial_guesses: List[ndarray]) → List[ndarray][source]

Search for saddle points z_c ∈ ℂⁿ satisfying ∇φ(z_c) = 0.

The search minimises the real function

F(u, v) = |∇φ(u + iv)|², u, v ∈ ℝⁿ

starting from (u₀, v₀) = (x₀, 0) for each real guess x₀.

Parameters:

analyzer (Analyzer) – Analyzer whose lambdified gradient func_grad is used. The phase must accept complex-valued arguments.
initial_guesses (list of ndarray) – Real starting points in ℝⁿ. Each is lifted to ℂⁿ by setting the imaginary part to zero.

Returns:

Unique saddle points found (deduplicated within 1e-6). Each array has shape (dim,) and dtype complex128.

Return type:

list of complex ndarray

class asymptotic.SingularityType(*values)[source]

Bases: Enum

Classification of critical points based on Hessian rank and higher-order derivatives.

The type determines which asymptotic formula applies to the stationary phase integral:

MORSE: Non-degenerate critical point (det H ≠ 0) Contribution scales as O(λ^(-n/2)) where n is the dimension
AIRY_1D: Corank-1 singularity with non-zero cubic term (1D case) Contribution scales as O(λ^(-1/3))
AIRY_2D: Corank-1 singularity with non-zero cubic term (2D case) Contribution scales as O(λ^(-5/6)) = O(λ^(-1/3-1/2))
PEARCEY: Corank-1 singularity with vanishing cubic but non-zero quartic term Contribution scales as O(λ^(-3/4)) = O(λ^(-1/4-1/2))
HIGHER_ORDER: More degenerate cases requiring special treatment Not implemented in this code

AIRY_1D = 'airy_1d'

AIRY_2D = 'airy_2d'

HIGHER_ORDER = 'higher_order'

MORSE = 'morse'

PEARCEY = 'pearcey'

class asymptotic.StationaryPhaseEvaluator(tolerance=1e-08)[source]

Bases: object

Computes asymptotic contributions from critical points for large parameter λ.

This class implements the standard stationary phase formulas for different types of critical points:

Morse points: Standard stationary phase with second-order corrections
Airy singularities (1D and 2D): Catastrophe integrals with exact formulas
Pearcey singularities: Quartic catastrophe integrals

The evaluation includes both leading-order terms and next-order corrections where applicable (primarily for Morse points). Each method returns an AsymptoticContribution object containing the computed terms.

Reference:

Wong, “Asymptotic Approximations of Integrals” (1989)
Olver, “Asymptotics and Special Functions” (1997)

tolerance

Numerical tolerance for detecting near-zero coefficients.

Type:: float

evaluate(cp: CriticalPoint, lam: float) → AsymptoticContribution[source]

Dispatch evaluation to the appropriate method based on singularity type.

Parameters:

cp – CriticalPoint object with all necessary geometric data.
lam – Large parameter λ in the oscillatory integral I(λ).

Returns:

AsymptoticContribution containing leading term, correction (if computed), and total value. For HIGHER_ORDER or unknown types, returns zero contribution with a warning.

asymptotic.StationaryPhaseVisualizer: alias of AsymptoticVisualizer

asymptotic.plot_contribution_decomposition(critical_points, analyzer, lambda_values=None, show_correction=True, show_coherent_sum=True, figsize=None)[source]

Decompose the asymptotic expansion I(λ) into its per-critical-point contributions and plot how each term scales with the large parameter λ.

This is the asymptotic.py counterpart of wkb.plot_amplitude_decomposition. While the WKB version displays spatial amplitude orders aₖ(x) on a grid, the asymptotic expansion has no spatial structure: its “terms” are complex numbers indexed by critical point and by expansion order (leading / correction). The natural display axis is therefore λ, shown on a log scale.

The function produces a figure with three stacked panels:

Panel 1 — Leading terms (one curve per critical point): |I₀ⱼ(λ)| = |leading_term| for each critical point j, plotted on a log-log scale. A theoretical reference line λ^(-p) is overlaid using the order_leading attribute of the first evaluated contribution. The label includes the critical-point coordinates, its singularity type, and its integration method.
Panel 2 — Correction terms (one curve per critical point, optional): |I₁ⱼ(λ)| = |correction_term|. For degenerate singularities (Airy, Pearcey) where the correction is not computed the curve is omitted and a note is added to the legend. Shown only when show_correction=True.
Panel 3 — Coherent total sum (optional): |I(λ)| = |Σⱼ (leading_termⱼ + correction_termⱼ)| — the coherent sum of all contributions, including interference between critical points. A second curve shows the incoherent upper bound Σⱼ |total_valueⱼ|. Shown only when show_coherent_sum=True.

Parameters:

critical_points (list of CriticalPoint) – Critical (or saddle) points obtained from Analyzer.analyze_point(). Each element must have its method attribute set to a concrete IntegralMethod (not AUTO). Pass a single-element list for a one-critical-point problem.
analyzer (Analyzer) – The Analyzer instance that produced the critical points. Its method attribute is used to label the figure title and to dispatch evaluation through AsymptoticEvaluator.
lambda_values (array_like or None, default None) – Sequence of positive λ values at which to evaluate all contributions. If None, defaults to np.logspace(0, 4, 60) (i.e. λ ∈ [1, 10000]). Values must be strictly positive.
show_correction (bool, default True) – If True, Panel 2 (correction terms) is included in the figure. Set to False when all critical points are degenerate (Airy/Pearcey) and correction terms are identically zero.
show_coherent_sum (bool, default True) – If True, Panel 3 (coherent sum) is included. Set to False when only one critical point is present and the total is already shown in Panel 1.
figsize (tuple or None, default None) – Passed directly to matplotlib.pyplot.subplots. If None, the height is chosen automatically as 4 inches per panel.

Returns:

fig (matplotlib.figure.Figure) – The figure object, so the caller can save or further customise it.
axes (list of matplotlib.axes.Axes) – The axes objects for each panel (length 1, 2, or 3 depending on the flags).

Raises:

ValueError – If critical_points is empty.
ValueError – If any CriticalPoint.method is None or AUTO.

Examples

Basic usage with two stationary-phase critical points:

>>> import sympy as sp
>>> import numpy as np
>>> x = sp.Symbol('x')
>>> phi   = x**4 - x**2          # two minima at x = ±1/√2
>>> amp   = sp.Integer(1)
>>> analyzer = Analyzer(phi, amp, [x],
...                     method=IntegralMethod.STATIONARY_PHASE)
>>> pts = analyzer.find_critical_points(
...     [np.array([ 0.7]), np.array([-0.7])]
... )
>>> cps = [analyzer.analyze_point(p) for p in pts]
>>> fig, axes = plot_contribution_decomposition(
...     cps, analyzer,
...     lambda_values=np.logspace(1, 4, 80),
... )
>>> fig.savefig("decomposition.png", dpi=150)

Suppress the correction panel (all Airy singularities):

>>> fig, axes = plot_contribution_decomposition(
...     cps, analyzer, show_correction=False
... )