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

[1]

Hörmander, L. The Analysis of Linear Partial Differential Operators I, Springer, 1983. Chapter 7: Oscillatory Integrals.

[2]

Olver, F. W. J. Asymptotics and Special Functions, Academic Press, 1974 (reprinted A K Peters, 1997).

[3]

Wong, R. Asymptotic Approximations of Integrals, Academic Press, 1989.

[4]

Bleistein, N. & Handelsman, R. Asymptotic Expansions of Integrals, Holt, Rinehart & Winston, 1975.

[5]

Berry, M. V. & Howls, C. J. “High orders of the Weyl expansion for quantum billiards”, Physical Review E 50(5), 3577–3595, 1994.

[6]

Delabaere, E. & Howls, C. J. “Global asymptotics for multiple integrals with boundaries”, Duke Mathematical Journal 112(2), 199–264, 2002.

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

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:

1. Initialize with symbolic SymPy expressions (method=AUTO by default).

2. Find critical / saddle points using find_critical_points().

3. Analyze each point using analyze_point().

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

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]

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)

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)

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

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)

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

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:

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

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:

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

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

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

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)

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)

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

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

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:

• leading_term: I₀(λ), real for real a and φ.

• correction_term: I₁(λ), the O(λ^(-n/2-1)) correction.

• total_value: I₀ + I₁.

• order_leading: n/2.

• method: IntegralMethod.LAPLACE.

Return type:

AsymptoticContribution with

Raises:

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

class asymptotic.SaddlePointEvaluator(tolerance: float = 1e-08)

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)

1. Start from real initial guesses x₀ ∈ ℝⁿ (supplied by the caller).

2. Analytically continue into ℂⁿ by minimising |∇φ(z)|² over the 2n real degrees of freedom (Re z, Im z), using scipy.optimize.minimize.

3. Accept a point z_c if |∇φ(z_c)|² < tolerance.

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

[1]

Bleistein & Handelsman, “Asymptotic Expansions of Integrals” (1975)

[2]

Delabaere & Howls, “Global asymptotics for multiple integrals with boundaries” (2002)

evaluate(cp: CriticalPoint, lam: float) → AsymptoticContribution

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]

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)

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)

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

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