# Copyright 2026 Philippe Billet assisted by LLMs in free mode: chatGPT, Qwen, Deepseek, Gemini, Claude, le chat Mistral.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
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.
"""
from imports import *
[docs]
class IntegralMethod(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
"""
STATIONARY_PHASE = "stationary_phase"
LAPLACE = "laplace"
SADDLE_POINT = "saddle_point"
AUTO = "auto"
[docs]
class SingularityType(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
"""
MORSE = "morse" # Non-degenerate (det H != 0)
AIRY_1D = "airy_1d" # Corank 1, cubic term != 0 (1D)
AIRY_2D = "airy_2d" # Corank 1, cubic term != 0 (2D)
PEARCEY = "pearcey" # Corank 1, cubic = 0, quartic != 0
HIGHER_ORDER = "higher_order"
[docs]
@dataclass
class CriticalPoint:
"""
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 λ → ∞.
Attributes:
position (np.ndarray): Coordinates of the critical point x_c.
phase_value (complex): Value of phase function φ(x_c).
amplitude_value (complex): Value of amplitude function a(x_c).
singularity_type (SingularityType): Classification determining which formula to use.
hessian_matrix (np.ndarray): The Hessian matrix ∇²φ at x_c.
hessian_inv (Optional[np.ndarray]): Inverse of the Hessian (for Morse points only).
hessian_det (float): Determinant of the Hessian.
signature (int): Number of negative eigenvalues of the Hessian (Morse index).
eigenvalues (np.ndarray): Eigenvalues of the Hessian matrix.
eigenvectors (np.ndarray): Eigenvectors of the Hessian matrix.
grad_amp (Optional[np.ndarray]): Gradient of amplitude ∇a at x_c.
hess_amp (Optional[np.ndarray]): Hessian of amplitude ∇²a at x_c.
phase_d3 (Optional[np.ndarray]): Rank-3 tensor of 3rd derivatives of φ.
phase_d4 (Optional[np.ndarray]): Rank-4 tensor of 4th derivatives of φ.
canonical_coefficients (Optional[Dict]): Coefficients for normal forms
(Airy/Pearcey canonical representations). Contains keys like 'cubic',
'quartic', 'quadratic_transverse' depending on singularity type.
"""
position: np.ndarray
phase_value: complex
amplitude_value: complex
singularity_type: SingularityType
hessian_matrix: np.ndarray
hessian_inv: Optional[np.ndarray] = None
hessian_det: float = 0.0
signature: int = 0
eigenvalues: np.ndarray = field(default_factory=lambda: np.array([]))
eigenvectors: np.ndarray = field(default_factory=lambda: np.array([]))
# Higher order derivatives (stored as numpy tensors)
grad_amp: Optional[np.ndarray] = None
hess_amp: Optional[np.ndarray] = None
phase_d3: Optional[np.ndarray] = None
phase_d4: Optional[np.ndarray] = None
canonical_coefficients: Optional[Dict] = None
# Integration method that produced this critical point
method: 'IntegralMethod' = None # set by the analyzer; forward ref resolved at runtime
[docs]
@dataclass
class AsymptoticContribution:
"""
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.
Attributes:
leading_term (complex): 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))
correction_term (complex): The next-order correction term:
for Morse: O(λ^(-n/2-1)),
for degenerate cases: typically 0j (not computed)
total_value (complex): Sum of leading_term + correction_term.
point (CriticalPoint): The source critical point for this contribution.
order_leading (float): The exponent p in the scaling λ^(-p) of the leading term.
method (IntegralMethod): The asymptotic method used to compute this contribution
(STATIONARY_PHASE, LAPLACE or SADDLE_POINT).
"""
leading_term: complex
correction_term: complex
total_value: complex
point: CriticalPoint
order_leading: float
method: IntegralMethod = IntegralMethod.STATIONARY_PHASE # default for backward compat
# --- Analyzer (Symbolic -> Numerical) ---
[docs]
class Analyzer:
"""
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.
Attributes:
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 (int): Dimension of the integration domain.
domain (Optional[List[Tuple]]): Optional bounds [(min, max), ...] for each variable.
tolerance (float): Numerical tolerance for detecting zeros and critical points.
cubic_threshold (float): Absolute threshold for distinguishing Airy vs Pearcey
singularities. Default: max(1e-5, 10*tolerance). Unused for LAPLACE.
method (IntegralMethod): Resolved integration method (never AUTO after __init__).
"""
def __init__(self, phase_expr, amplitude_expr, variables, domain=None,
tolerance=1e-6, cubic_threshold=None,
method: IntegralMethod = IntegralMethod.AUTO):
"""
Initialize the analyzer.
Args:
phase_expr: SymPy expression for phase φ(x).
The nature of φ (real / imaginary / complex) determines which
asymptotic method is appropriate; use method=AUTO to detect it
automatically.
amplitude_expr: SymPy expression for amplitude a(x). Can be complex.
variables: List of SymPy symbols [x, y, ...] or a single symbol for 1D.
domain: Optional list of tuples [(min, max), ...] specifying search bounds
for each variable when finding critical points.
tolerance: Numerical tolerance for zero-detection and optimization
(default: 1e-6).
cubic_threshold: Absolute threshold for classifying cubic terms
(STATIONARY_PHASE only). If None, defaults to max(1e-5, 10*tolerance).
method: One of IntegralMethod.{AUTO, STATIONARY_PHASE, LAPLACE,
SADDLE_POINT}. AUTO (default) inspects φ symbolically and
resolves to one of the three concrete values.
"""
self.phase_expr = phase_expr
self.amplitude_expr = amplitude_expr
self.variables = list(variables) if isinstance(variables, (list, tuple)) else [variables]
self.dim = len(self.variables)
self.domain = domain
self.tolerance = tolerance
self.cubic_threshold = cubic_threshold if cubic_threshold is not None else max(1e-5, 10 * tolerance)
# Resolve AUTO before preparing derivatives so that _detect_method
# can inspect the raw expression before any lambdification.
if method == IntegralMethod.AUTO:
self.method = self._detect_method(phase_expr)
else:
self.method = method
self._prepare_derivatives()
self._create_numerical_functions()
# ------------------------------------------------------------------
# Method auto-detection
# ------------------------------------------------------------------
def _detect_method(self, phase_expr) -> IntegralMethod:
"""
Inspect the phase expression symbolically to select the integration method.
Strategy
--------
1. Symbolic test (fast, exact when SymPy can simplify):
- Compute re_part = re(φ) and im_part = im(φ) after
assuming all variables are real (refine with Q.real).
- If re_part simplifies to zero → LAPLACE
- If im_part simplifies to zero → STATIONARY_PHASE
- Otherwise → SADDLE_POINT
2. Numerical fallback (if symbolic test is inconclusive, i.e. SymPy
cannot decide):
- Sample n_samples random real points in [-2, 2]^dim.
- Evaluate φ numerically at each sample.
- If max|Re φ| < tol * max|Im φ| → LAPLACE
- If max|Im φ| < tol * max|Re φ| → STATIONARY_PHASE
- Otherwise → SADDLE_POINT
The numerical threshold uses tol = 1e-6 relative to the dominant part.
Returns
-------
IntegralMethod
One of STATIONARY_PHASE, LAPLACE, or SADDLE_POINT (never AUTO).
"""
# --- Step 1: symbolic test ---
real_assumptions = {v: True for v in self.variables}
# Replace variables with real-stamped symbols for re / im
real_vars = [Symbol(str(v), real=True) for v in self.variables]
expr_real = phase_expr.subs(dict(zip(self.variables, real_vars)))
try:
re_part = simplify(re(expr_real))
im_part = simplify(im(expr_real))
re_is_zero = (re_part == S.Zero)
im_is_zero = (im_part == S.Zero)
if re_is_zero and not im_is_zero:
return IntegralMethod.LAPLACE
if im_is_zero and not re_is_zero:
return IntegralMethod.STATIONARY_PHASE
if re_is_zero and im_is_zero:
# Constant zero phase — stationary phase is the safest default
return IntegralMethod.STATIONARY_PHASE
# Both parts are symbolically non-zero → fall through to numeric
except Exception:
pass # SymPy failed to simplify; proceed to numerical fallback
# --- Step 2: numerical fallback ---
n_samples = 12
rng = np.random.default_rng(seed=0) # deterministic seed for reproducibility
samples = rng.uniform(-2.0, 2.0, size=(n_samples, self.dim))
func_phase_num = lambdify(tuple(self.variables), phase_expr, 'numpy')
re_magnitudes = []
im_magnitudes = []
for pt in samples:
try:
val = complex(func_phase_num(*pt))
re_magnitudes.append(abs(val.real))
im_magnitudes.append(abs(val.imag))
except Exception:
pass
if not re_magnitudes:
# Could not evaluate at all — default to stationary phase
warnings.warn(
"AUTO method detection: could not evaluate φ numerically. "
"Defaulting to STATIONARY_PHASE.",
RuntimeWarning,
)
return IntegralMethod.STATIONARY_PHASE
max_re = max(re_magnitudes)
max_im = max(im_magnitudes)
scale = max(max_re, max_im, 1e-30)
rel_tol = 1e-6
if max_re < rel_tol * scale:
detected = IntegralMethod.LAPLACE
elif max_im < rel_tol * scale:
detected = IntegralMethod.STATIONARY_PHASE
else:
detected = IntegralMethod.SADDLE_POINT
return detected
def _prepare_derivatives(self):
"""
Symbolically compute all necessary derivatives of phase and amplitude functions.
Computes and stores as SymPy expressions:
1. Gradient ∇φ and Hessian ∇²φ of phase
2. Gradient ∇a and Hessian ∇²a of amplitude
3. Third-order tensor D3[i,j,k] = ∂³φ/∂xi∂xj∂xk (for Airy classification
and second-order Morse correction)
4. Fourth-order tensor D4[i,j,k,l] = ∂⁴φ/∂xi∂xj∂xk∂xl (for Pearcey
classification and second-order Morse correction)
These symbolic derivatives are later converted to numerical functions via lambdify.
"""
# 1. Phase Gradient & Hessian
self.grad_sym = [diff(self.phase_expr, v) for v in self.variables]
self.hess_sym = [[diff(self.phase_expr, v1, v2) for v2 in self.variables]
for v1 in self.variables]
# 2. Amplitude Gradient & Hessian
self.grad_amp_sym = [diff(self.amplitude_expr, v) for v in self.variables]
self.hess_amp_sym = [[diff(self.amplitude_expr, v1, v2) for v2 in self.variables]
for v1 in self.variables]
# 3. Higher order tensors for Phase
# D3 Tensor (Rank 3)
self.d3_indices = []
self.d3_sym = []
import itertools
for idx in itertools.product(range(self.dim), repeat=3):
self.d3_indices.append(idx)
var_seq = [self.variables[i] for i in idx]
self.d3_sym.append(diff(self.phase_expr, *var_seq))
# D4 Tensor (Rank 4) - Required for 2nd order Morse correction
self.d4_indices = []
self.d4_sym = []
for idx in itertools.product(range(self.dim), repeat=4):
self.d4_indices.append(idx)
var_seq = [self.variables[i] for i in idx]
self.d4_sym.append(diff(self.phase_expr, *var_seq))
def _create_numerical_functions(self):
"""
Convert symbolic expressions to fast numerical functions using SymPy's lambdify.
Creates NumPy-compatible functions for efficient numerical evaluation:
- func_phase, func_amp: Evaluate φ(x) and a(x)
- func_grad, func_hess: Evaluate ∇φ and ∇²φ
- func_grad_amp, func_hess_amp: Evaluate ∇a and ∇²a
- func_d3, func_d4: Evaluate third and fourth order derivatives of φ
These lambdified functions are much faster than evaluating SymPy expressions directly.
"""
vars_tuple = tuple(self.variables)
self.func_phase = lambdify(vars_tuple, self.phase_expr, 'numpy')
self.func_amp = lambdify(vars_tuple, self.amplitude_expr, 'numpy')
self.func_grad = lambdify(vars_tuple, self.grad_sym, 'numpy')
self.func_hess = lambdify(vars_tuple, self.hess_sym, 'numpy')
self.func_grad_amp = lambdify(vars_tuple, self.grad_amp_sym, 'numpy')
self.func_hess_amp = lambdify(vars_tuple, self.hess_amp_sym, 'numpy')
self.func_d3 = lambdify(vars_tuple, self.d3_sym, 'numpy')
self.func_d4 = lambdify(vars_tuple, self.d4_sym, 'numpy')
[docs]
def find_critical_points(self, initial_guesses=None) -> List[np.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).
Args:
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.
"""
from caustics import find_critical_points_numerical
if initial_guesses is None:
initial_guesses = [np.zeros(self.dim)]
if self.domain:
centers = [0.5 * (d[0] + d[1]) for d in self.domain]
initial_guesses.append(np.array(centers))
return find_critical_points_numerical(
grad_func=self.func_grad,
initial_guesses=initial_guesses,
tolerance=self.tolerance,
domain=self.domain,
)
[docs]
def analyze_point(self, 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)
Args:
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.
"""
args = tuple(xc)
# Evaluate Hessian — complex-safe: cast to complex array so that
# complex saddle-point coordinates do not silently drop imaginary parts.
H = np.array(self.func_hess(*args), dtype=complex)
# Eigendecomposition — use np.linalg.eig for complex matrices
# (eigh requires Hermitian; a complex Hessian at a saddle point is
# symmetric but not necessarily Hermitian).
if np.iscomplexobj(H) and np.any(np.imag(H) != 0):
vals, vecs = np.linalg.eig(H)
else:
H = np.real(H)
vals, vecs = np.linalg.eigh(H)
# Reconstruct higher order tensors from flattened symbolic output
d3_flat = self.func_d3(*args)
D3 = np.zeros((self.dim,)*3, dtype=complex)
for k, idx in enumerate(self.d3_indices):
D3[idx] = d3_flat[k]
d4_flat = self.func_d4(*args)
D4 = np.zeros((self.dim,)*4, dtype=complex)
for k, idx in enumerate(self.d4_indices):
D4[idx] = d4_flat[k]
grad_a = np.array(self.func_grad_amp(*args), dtype=complex)
hess_a = np.array(self.func_hess_amp(*args), dtype=complex)
det = np.prod(vals)
# Rank: number of eigenvalues with non-negligible magnitude
rank = int(np.sum(np.abs(vals) > self.tolerance))
# Signature: number of eigenvalues with strictly negative real part
signature = int(np.sum(np.real(vals) < -self.tolerance))
cp = CriticalPoint(
position=np.asarray(xc),
phase_value=complex(self.func_phase(*args)),
amplitude_value=complex(self.func_amp(*args)),
singularity_type=SingularityType.MORSE, # Default, may be overridden below
hessian_matrix=H,
hessian_det=complex(det),
signature=signature,
eigenvalues=vals,
eigenvectors=vecs,
grad_amp=grad_a,
hess_amp=hess_a,
phase_d3=D3,
phase_d4=D4,
method=self.method,
)
if rank == self.dim:
cp.singularity_type = SingularityType.MORSE
cp.hessian_inv = np.linalg.inv(H)
# For Laplace's method the Hessian must be positive definite (minimum of φ).
if self.method == IntegralMethod.LAPLACE and np.any(np.real(vals) <= self.tolerance):
warnings.warn(
f"Laplace method: critical point at {xc} has a non-positive Hessian "
"eigenvalue (it may be a saddle point or maximum). "
"The Laplace approximation requires a strict minimum of φ.",
RuntimeWarning,
)
elif self.dim == 1 and rank == 0: # 1D degenerate
coeffs = self._project_degenerate_coeffs(cp)
cp.canonical_coefficients = coeffs
if abs(coeffs['cubic']) > self.cubic_threshold:
cp.singularity_type = SingularityType.AIRY_1D
elif abs(coeffs['quartic']) > self.tolerance:
cp.singularity_type = SingularityType.PEARCEY
else:
cp.singularity_type = SingularityType.HIGHER_ORDER
elif self.dim == 2 and rank == 1: # 2D corank 1
coeffs = self._project_degenerate_coeffs(cp)
cp.canonical_coefficients = coeffs
if abs(coeffs['cubic']) > self.cubic_threshold:
cp.singularity_type = SingularityType.AIRY_2D
elif abs(coeffs['quartic']) > self.tolerance:
cp.singularity_type = SingularityType.PEARCEY
else:
cp.singularity_type = SingularityType.HIGHER_ORDER
else:
cp.singularity_type = SingularityType.HIGHER_ORDER
return cp
def _project_degenerate_coeffs(self, cp: CriticalPoint) -> Dict[str, float]:
"""
Projects derivatives onto the eigenvectors to find canonical form coefficients.
This handles cases where the singularity is not aligned with the axes
(e.g., phi = (x+y)^3). It identifies the 'null' direction and computes
directional derivatives along it.
Returns:
Dictionary with 'cubic', 'quartic' and 'quadratic_transverse' coefficients.
"""
null_idx = np.argmin(np.abs(cp.eigenvalues))
v_null = cp.eigenvectors[:, null_idx]
# AIRY CORRECTION: alpha = D^3/2 for phi ~ alpha * u^3/3 (since d^3(alpha * u^3/3)/du^3 = 2 * alpha)
alpha = np.einsum('ijk,i,j,k->', cp.phase_d3, v_null, v_null, v_null) / 2.0
# PEARCEY CORRECTION: gamma = D^4/24 for phi ~ gamma * u^4/4
# (since d^4(gamma * u^4/4)/du^4 = 6 * gamma -> gamma = D^4/6,
# but the canonical normal form uses u^4/4 -> coefficient = D^4/24)
# gamma_coeff = np.einsum('ijkl,i,j,k,l->', cp.phase_d4, v_null, v_null, v_null, v_null) / 24.0
gamma_coeff = np.einsum('ijkl,i,j,k,l->', cp.phase_d4, v_null, v_null, v_null, v_null) / 6.0
# Transverse quadratic term (2D only)
quadratic_transverse = None
if self.dim > 1:
non_null_idxs = np.where(np.abs(cp.eigenvalues) > self.tolerance)[0]
if len(non_null_idxs) > 0:
# The eigenvalue IS the coefficient beta in phi ~ beta*v^2/2
# (since d²phi/dv² = beta for this form)
quadratic_transverse = cp.eigenvalues[non_null_idxs[0]]
return {
'cubic': alpha, # alpha for phi = alpha * u^3/3
'quartic': gamma_coeff, # gamma for phi = gamma * u^4/4
'quadratic_transverse': quadratic_transverse
}
# --- Evaluator (Asymptotic Calculation) ---
[docs]
class StationaryPhaseEvaluator:
"""
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)
Attributes:
tolerance (float): Numerical tolerance for detecting near-zero coefficients.
"""
def __init__(self, tolerance=1e-8):
self.tolerance = tolerance # ← Addition required
[docs]
def evaluate(self, cp: CriticalPoint, lam: float) -> AsymptoticContribution:
"""
Dispatch evaluation to the appropriate method based on singularity type.
Args:
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.
"""
if cp.singularity_type == SingularityType.MORSE:
result = self._eval_morse_order2(cp, lam)
elif cp.singularity_type == SingularityType.AIRY_1D:
result = self._eval_airy_1d(cp, lam)
elif cp.singularity_type == SingularityType.AIRY_2D:
result = self._eval_airy_2d(cp, lam)
elif cp.singularity_type == SingularityType.PEARCEY:
result = self._eval_pearcey(cp, lam)
else: # HIGHER_ORDER or unknown type
warnings.warn(
f"Unhandled singularity type {cp.singularity_type.value} at {cp.position}. "
f"Returning zero contribution (no asymptotic formula available).",
RuntimeWarning
)
result = AsymptoticContribution(
leading_term=0j,
correction_term=0j,
total_value=0j,
point=cp,
order_leading=float('inf'), # Indicates negligible contribution
method=IntegralMethod.STATIONARY_PHASE,
)
result.method = IntegralMethod.STATIONARY_PHASE
return result
def _eval_morse_order2(self, cp: CriticalPoint, lam: float) -> AsymptoticContribution:
"""
Evaluate the asymptotic contribution for a non-degenerate (Morse) critical point.
This method implements the standard stationary phase formula with second-order
correction terms for oscillatory integrals of the form:
I(λ) = ∫ a(x) exp(iλφ(x)) dx
as λ → ∞. The asymptotic expansion is:
I(λ) ≈ I₀(λ) + I₁(λ) + O(λ^(-n/2-2))
where I₀ is the leading term (order λ^(-n/2)) and I₁ is the first correction
(order λ^(-n/2-1)).
Leading Term (Order 0)
----------------------
The dominant contribution from the critical point x_c where ∇φ(x_c) = 0:
I₀(λ) = (2π)^(n/2) / (λ^(n/2) √|det H|) × exp(iλφ(x_c)) × a(x_c) × exp(iπμ/4)
Components:
- (2π/λ)^(n/2): Gaussian prefactor from the quadratic approximation
- √|det H|: Determinant of the Hessian matrix H = ∇²φ(x_c)
- exp(iλφ(x_c)): Rapid oscillation at the critical point
- a(x_c): Amplitude function evaluated at the critical point
- exp(iπμ/4): Maslov phase correction, where μ = n - 2σ is the Morse index
(n = dimension, σ = signature = number of negative eigenvalues of H)
Correction Term (Order 1)
-------------------------
The next-order contribution accounts for:
1. Non-constant amplitude (amplitude derivatives)
2. Cubic phase anharmonicity (third derivatives of φ)
3. Quartic phase anharmonicity (fourth derivatives of φ)
I₁(λ) = I₀(λ) / (iλ) × C
where the correction factor C is:
C = (1/2) Tr(H⁻¹ ∇²a) - (1/2) ⟨H⁻¹∇a, V⟩ + (a(x_c)/24) (5S₃ - 3S₄)
Term breakdown:
1. Amplitude Laplacian term: (1/2) Tr(H⁻¹ ∇²a)
- Captures the effect of amplitude curvature at the critical point
- H⁻¹ "twists" the Laplacian by the phase geometry
2. Mixed amplitude-phase term: -(1/2) ⟨H⁻¹∇a, V⟩
- Couples amplitude gradient with cubic phase nonlinearity
- V_k = Σᵢⱼ (H⁻¹)ᵢⱼ ∂³φ/∂xᵢ∂xⱼ∂xₖ
3. Pure phase anharmonicity: (a(x_c)/24) (5S₃ - 3S₄)
- S₄: Quartic term = Σᵢⱼₖₗ (H⁻¹)ᵢⱼ (H⁻¹)ₖₗ ∂⁴φ/∂xᵢ∂xⱼ∂xₖ∂xₗ
- S₃: Cubic term = Σᵢⱼₖₗₘₙ (H⁻¹)ᵢⱼ (H⁻¹)ₖₗ (H⁻¹)ₘₙ D³φᵢₖₘ D³φⱼₗₙ
- The coefficients 5 and 3 come from Feynman diagram combinatorics
Mathematical Background
-----------------------
The correction terms arise from expanding the integrand to higher orders in
(x - x_c) around the critical point and performing Gaussian integrals. The
coefficients are determined by the topology of Feynman diagrams:
- S₃ corresponds to "theta graph" diagrams (three-loop)
- S₄ corresponds to "sunset" diagrams (two-loop with quartic vertex)
The factors of 1/2 in the amplitude terms come from the expansion of the
Gaussian measure, while the 5 and 3 in the phase term arise from diagram
symmetry factors.
Parameters
----------
cp : CriticalPoint
Critical point with non-zero Hessian determinant (det H ≠ 0).
Must contain: position, phase_value, amplitude_value, hessian_matrix,
hessian_inv, signature, grad_amp, hess_amp, phase_d3, phase_d4.
lam : float
Large frequency parameter λ. The asymptotic approximation improves
as λ → ∞. Typically valid for λ ≳ 10.
Returns
-------
AsymptoticContribution
Object containing:
- leading_term: I₀(λ), the dominant O(λ^(-n/2)) contribution
- correction_term: I₁(λ), the O(λ^(-n/2-1)) correction
- total_value: I₀(λ) + I₁(λ)
- point: Reference to the input critical point
- order_leading: n/2 (the decay exponent)
Notes
-----
The correction term becomes negligible for large λ. The ratio
|I₁/I₀| ~ O(λ⁻¹) should decrease linearly on a log-log plot, which
can be verified using the convergence diagnostic tools.
For dimension n=2, the leading term scales as O(λ⁻¹) and the correction
as O(λ⁻²), providing rapid asymptotic convergence.
References
----------
.. [1] Hörmander, L. "The Analysis of Linear Partial Differential Operators I"
(1983), Chapter 7: Oscillatory Integrals
.. [2] Berry, M.V. & Howls, C.J. "High orders of the Weyl expansion for quantum
billiards" Physical Review E 50.5 (1994): 3577-3595
.. [3] Wong, R. "Asymptotic Approximations of Integrals" (1989), Chapter 2
Examples
--------
>>> # For a Gaussian phase φ = x²/2 + y²/2 with constant amplitude a = 1
>>> # at the critical point (0, 0), the leading term is:
>>> # I₀(λ) = 2π/λ (exact for Gaussian)
>>> evaluator = StationaryPhaseEvaluator()
>>> contribution = evaluator._eval_morse_order2(cp, lam=100)
>>> print(f"Leading: {contribution.leading_term:.4e}")
>>> print(f"Correction: {contribution.correction_term:.4e}")
>>> print(f"Ratio: {abs(contribution.correction_term/contribution.leading_term):.2%}")
"""
# dim = cp.position.shape[0]
dim = len(cp.position)
# ============================================================================
# LEADING TERM (Order λ^(-n/2))
# ============================================================================
# Compute the dominant Gaussian contribution from the quadratic approximation
# of the phase near the critical point.
# Gaussian prefactor: (2π/λ)^(n/2)
# This comes from the n-dimensional Gaussian integral formula
prefactor = (2 * np.pi / lam) ** (dim / 2.0)
# Maslov phase: exp(iπμ/4) where μ = n - 2σ (Morse index)
# Accounts for the topology of the phase function at the critical point
# σ = signature = number of negative eigenvalues of the Hessian
maslov = np.exp(1j * np.pi / 4 * (dim - 2 * cp.signature))
# Rapid oscillatory factor: exp(iλφ(x_c))
# This is the phase evaluated at the critical point
phase_osc = np.exp(1j * lam * cp.phase_value)
# Geometric factor: 1/√|det H|
# The Hessian determinant measures the "curvature volume" at the critical point
denom = np.sqrt(np.abs(cp.hessian_det))
# Amplitude at critical point: a(x_c)
leading_amp = cp.amplitude_value
# Combine all factors for the leading term
term_0 = (prefactor / denom) * phase_osc * maslov * leading_amp
# ============================================================================
# CORRECTION TERM (Order λ^(-n/2-1))
# ============================================================================
# Compute next-order corrections from amplitude derivatives and phase
# anharmonicity (cubic and quartic terms in the Taylor expansion).
# Inverse Hessian matrix: H⁻¹ = (∇²φ)⁻¹
# Used to "propagate" corrections through the phase geometry
H_inv = cp.hessian_inv
# ------------------------------------------------------------------------
# Term 1: Amplitude Laplacian Contribution
# ------------------------------------------------------------------------
# Measures how the amplitude curvature affects the integral
# Formula: (1/2) Tr(H⁻¹ ∇²a)
#
# Physical interpretation: If the amplitude has negative curvature along
# directions where the phase is flat (small eigenvalues of H), this term
# can become significant.
term_amp = 0.5 * np.einsum('ij,ij->', H_inv, cp.hess_amp)
# ------------------------------------------------------------------------
# Term 2: Mixed Amplitude-Phase Contribution
# ------------------------------------------------------------------------
# Couples the amplitude gradient with cubic phase terms
#
# Step 1: Contract H⁻¹ with D³φ to get effective vector V_k
# V_k = Σᵢⱼ (H⁻¹)ᵢⱼ ∂³φ/∂xᵢ∂xⱼ∂xₖ
# This vector represents the "cubic force" felt by the amplitude gradient
V = np.einsum('ij,ijk->k', H_inv, cp.phase_d3)
# Step 2: Inner product of (H⁻¹∇a) with V, scaled by -1/2
# Formula: -(1/2) ⟨H⁻¹∇a, V⟩
#
# Physical interpretation: If amplitude increases along directions where
# phase has strong cubic nonlinearity, this coupling enhances the contribution
term_mix = -0.5 * np.dot(np.dot(H_inv, cp.grad_amp), V)
# ------------------------------------------------------------------------
# Term 3: Pure Phase Anharmonicity
# ------------------------------------------------------------------------
# Captures corrections from non-quadratic phase terms (cubic and quartic)
# These arise from expanding exp(iλφ(x)) beyond the Gaussian approximation
# S₄: Quartic contraction
# Contract the fourth derivative tensor D⁴φ with two copies of H⁻¹
# Formula: Σᵢⱼₖₗ (H⁻¹)ᵢⱼ (H⁻¹)ₖₗ ∂⁴φ/∂xᵢ∂xⱼ∂xₖ∂xₗ
#
# Corresponds to Feynman diagrams with a single quartic vertex
S4 = np.einsum('ij,kl,ijkl->', H_inv, H_inv, cp.phase_d4)
# S₃: Cubic contraction (Theta graph)
# Contract two copies of D³φ with three copies of H⁻¹
# Formula: Σᵢⱼₖₗₘₙ (H⁻¹)ᵢⱼ (H⁻¹)ₖₗ (H⁻¹)ₘₙ D³φᵢₖₘ D³φⱼₗₙ
#
# Corresponds to Feynman diagrams with two cubic vertices connected in a loop
# (the "theta graph" topology)
S3 = np.einsum('ij,kl,mn,ikm,jln->', H_inv, H_inv, H_inv, cp.phase_d3, cp.phase_d3)
# Combine cubic and quartic contributions with diagram symmetry factors
# Formula: (a₀/24) (5S₃ - 3S₄)
#
# The coefficients 5 and 3 arise from:
# - Combinatorial factors in the Taylor expansion
# - Symmetry factors of Feynman diagrams (vertex permutations)
# - Wick's theorem for Gaussian integrals
term_phase = (cp.amplitude_value / 24.0) * (5.0 * S3 - 3.0 * S4)
# Sum all three correction contributions
correction_factor = term_amp + term_mix + term_phase
# Scale correction by 1/(iλ) relative to leading term
# The factor of i comes from ∫ x² exp(iλφ) dx ∝ -i/λ ∫ ∂²/∂λ² exp(iλφ) dx
val_correction = (prefactor / denom) * phase_osc * maslov * (correction_factor / (1j * lam))
# ============================================================================
# RETURN ASYMPTOTIC CONTRIBUTION
# ============================================================================
return AsymptoticContribution(
leading_term=term_0,
correction_term=val_correction,
total_value=term_0 + val_correction,
point=cp,
order_leading=dim/2.0
)
def _eval_airy_1d(self, cp: CriticalPoint, lam: float) -> AsymptoticContribution:
"""
Evaluate the 1D Airy catastrophe integral contribution.
For a canonical Airy integral of the form:
∫ exp(iλ α x³/3) dx
The exact asymptotic formula is:
I(λ) = 2π Ai(0) × (3λ|α|)^(-1/3) × exp(iπ/6 × sign(α))
where:
- Ai(0) ≈ 0.355028... is the Airy function at zero
- α is the cubic coefficient in the canonical form φ ~ α x³/3
- The scaling is O(λ^(-1/3)), which is slower decay than Morse O(λ^(-1/2))
The Maslov phase exp(iπ/6 × sign(α)) accounts for the orientation of the
integration contour in the complex plane:
∫ exp(i t³/3) dt = 2π Ai(0)
∫ exp(-i t³/3) dt = 2π Ai(0) × exp(-iπ/3)
Reference:
- Olver, "Asymptotics and Special Functions" (1997), §7.3
- Hörmander, "The Analysis of Linear Partial Differential Operators I" (1983)
Args:
cp: CriticalPoint classified as AIRY_1D with canonical coefficients.
lam: Large parameter λ.
Returns:
AsymptoticContribution with exact Airy scaling O(λ^(-1/3)).
No correction term is computed (set to 0j).
"""
coeffs = cp.canonical_coefficients
alpha = coeffs['cubic'] # Coefficient of x³/3 in the normal form
if abs(alpha) < self.tolerance:
warnings.warn("Cubic coefficient near zero in Airy evaluation")
return AsymptoticContribution(0j, 0j, 0j, cp, 1/3)
# Exact value of Ai(0)
Ai0 = airy(0)[0] # ≈ 0.3550280538878172
# Scale factor: (λ|α|)^(-1/3)
# Derivation: sub t = (λ|α|)^(1/3)·x → (λ|α|)^(-1/3) ∫ exp(it³/3) dt
# = (λ|α|)^(-1/3) · 2πAi(0) [Olver, "Asymptotics & Special Functions", §7.3]
# NOTE: (3λ|α|)^(-1/3) is WRONG — there is no factor of 3 here.
scale = (lam * abs(alpha)) ** (-1/3)
# No Maslov phase: ∫ exp(iλαx³/3) dx over ℝ is purely real for any sign of α,
# because cos(λαx³/3) is even (contributes) and sin(λαx³/3) is odd (cancels).
# Total contribution
val = 2 * np.pi * Ai0 * scale * cp.amplitude_value
return AsymptoticContribution(
leading_term=val,
correction_term=0j,
total_value=val,
point=cp,
order_leading=1/3 # λ^(-1/3)
)
def _eval_airy_2d(self, cp: CriticalPoint, lam: float) -> AsymptoticContribution:
"""
Evaluate the 2D Airy catastrophe integral contribution (corank 1).
For a canonical 2D Airy integral of the form:
∫∫ exp(iλ [α u³/3 + β v²/2]) du dv
The asymptotic formula is obtained by combining a 1D Airy integral
with a transverse 1D Gaussian integral:
I(λ) = [2π Ai(0) (3λ|α|)^(-1/3) e^{iπ/6 sign(α)}] * [sqrt(2π/(λ|β|)) e^{iπ/4 sign(β)}]
where:
- α is the cubic coefficient of the degenerate direction.
- β is the quadratic coefficient of the transverse (non-degenerate) direction.
- The scaling is O(λ^(-5/6)), combining O(λ^(-1/3)) and O(λ^(-1/2)).
Args:
cp: CriticalPoint classified as AIRY_2D with canonical coefficients.
lam: Large parameter λ.
Returns:
AsymptoticContribution with scaling O(λ^(-5/6)).
No correction term is computed (set to 0j).
"""
coeffs = cp.canonical_coefficients
alpha = coeffs['cubic']
beta = coeffs['quadratic_transverse']
# Transverse Gaussian Integral
# For ∫ exp(iλβv²/2) dv = √(2π/(λβ)) exp(iπ/4 sign(β))
scale_v = np.sqrt(2 * np.pi / (lam * np.abs(beta)))
phase_v = np.exp(1j * np.pi/4 * np.sign(beta))
# Degenerate Airy Integral (REAL, no Maslov phase):
# ∫ exp(iλα u³/3) du = 2π Ai(0) · (λ|α|)^(-1/3)
# Sub t = (λ|α|)^(1/3)·u → (λ|α|)^(-1/3) ∫ exp(it³/3) dt = (λ|α|)^(-1/3)·2πAi(0)
# The integral over ℝ is REAL: cos(λu³/3) is even → contributes; sin is odd → cancels.
# NOTE: (3λ|α|)^(-1/3) would be wrong — there is NO factor of 3.
Ai0 = airy(0)[0]
scale_u = 2 * np.pi * Ai0 * (lam * np.abs(alpha))**(-1.0/3.0)
# No phase_u: integral over ℝ of exp(iλαu³/3) is purely real for any sign of α.
val = (cp.amplitude_value * np.exp(1j * lam * cp.phase_value) * scale_u * scale_v * phase_v)
return AsymptoticContribution(
leading_term=val,
correction_term=0j,
total_value=val,
point=cp,
order_leading=5.0/6.0 # 1/3 (Airy) + 1/2 (Transverse Gaussian)
)
def _eval_pearcey(self, cp: CriticalPoint, lam: float) -> AsymptoticContribution:
"""
Evaluate the Pearcey catastrophe integral contribution (corank 1, quartic).
For a canonical Pearcey integral of the form:
∫∫ exp(iλ [γ u⁴/4 + β v²/2]) du dv
The asymptotic formula uses the specific value of the Pearcey integral at the origin:
I(λ) = [0.5 Γ(1/4) (λ|γ|)^(-1/4) e^{iπ/8 sign(γ)}] * [sqrt(2π/(λ|β|)) e^{iπ/4 sign(β)}]
where:
- γ is the quartic coefficient in the canonical form φ ~ γ u⁴/4.
- β is the quadratic coefficient of the transverse direction.
- The scaling is O(λ^(-3/4)), combining O(λ^(-1/4)) and O(λ^(-1/2)).
Note:
This evaluates the "cusp" catastrophe (A3) at its singular point.
Args:
cp: CriticalPoint classified as PEARCEY with canonical coefficients.
lam: Large parameter λ.
Returns:
AsymptoticContribution with scaling O(λ^(-3/4)).
No correction term is computed (set to 0j).
"""
coeffs = cp.canonical_coefficients
gamma_coeff = coeffs['quartic']
beta_coeff = coeffs['quadratic_transverse']
if abs(gamma_coeff) < self.tolerance or abs(beta_coeff) < self.tolerance:
warnings.warn("Near-zero coefficients in Pearcey evaluation")
return AsymptoticContribution(0j, 0j, 0j, cp, 0.75)
# Exact asymptotic constant for ∫ exp(iλγu⁴/4) du:
# Sub t = (λ|γ|)^{1/4}·u → (λ|γ|)^{-1/4} ∫ exp(it⁴/4) dt
# ∫_{-∞}^∞ exp(it⁴/4) dt = 4^{1/4} · ∫ exp(iv⁴) dv [sub v = t/4^{1/4}]
# = 4^{1/4} · (1/2)·Γ(1/4)·exp(iπ/8)
# Therefore: ∫ exp(iλγu⁴/4) du = (4/(λ|γ|))^{1/4} · (1/2)·Γ(1/4)·exp(iπ sign(γ)/8)
# NOTE: (1/(λ|γ|))^{1/4} is WRONG — the correct factor is (4/(λ|γ|))^{1/4} = √2/(λ|γ|)^{1/4}
pearcey_factor = (4.0 / (lam * abs(gamma_coeff)))**0.25 * 0.5 * gamma(0.25)
# Transverse Gaussian factor: ∫ exp(iλβv²/2) dv = √(2π/(λ|β|)) exp(iπ sign(β)/4)
gaussian_factor = np.sqrt(2.0 * np.pi / (lam * abs(beta_coeff)))
# Maslov phases
maslov_degen = np.exp(1j * np.pi * np.sign(gamma_coeff) / 8.0)
maslov_trans = np.exp(1j * np.pi * np.sign(beta_coeff) / 4.0)
leading = (cp.amplitude_value * np.exp(1j * lam * cp.phase_value) *
pearcey_factor * gaussian_factor * maslov_degen * maslov_trans)
return AsymptoticContribution(
leading_term=leading,
correction_term=0j,
total_value=leading,
point=cp,
order_leading=0.75
)
[docs]
class LaplaceEvaluator:
"""
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.
"""
[docs]
def evaluate(self, 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³φ)
Args:
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:
AsymptoticContribution: 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.
Raises:
ValueError: If the Hessian is singular (det H ≈ 0).
"""
n = len(cp.position)
inv_lam = 1.0 / lam
# --- Leading term (Order 0) ---
# I₀ = a(x_c) · exp(-λ φ(x_c)) · (2π/λ)^(n/2) / √|det H|
exponent = np.exp(-lam * np.real(cp.phase_value))
det_h_abs = np.abs(cp.hessian_det)
if det_h_abs < 1e-15:
raise ValueError(
"Hessian is singular at the critical point: Laplace method fails. "
"Use a higher-order evaluator for degenerate critical points."
)
prefactor = (2.0 * np.pi * inv_lam) ** (n / 2.0) / np.sqrt(det_h_abs)
term0 = cp.amplitude_value * exponent * prefactor
# --- Second-order corrections (O(1/λ)) ---
# Retrieve or recompute the inverse Hessian.
h_inv = cp.hessian_inv if cp.hessian_inv is not None else np.linalg.inv(cp.hessian_matrix)
# 1. Amplitude curvature: (1/2) Tr(H⁻¹ ∇²a)
term_lap_a = 0.5 * np.einsum('ij,ij', h_inv, cp.hess_amp)
# 2. Amplitude-gradient / cubic-phase coupling:
# −(1/2) Σᵢⱼₖ (∇a)ᵢ (H⁻¹)ⱼₖ (H⁻¹)ₖₗ D³φⱼₖₗ
# Rewritten as −(1/2) ⟨∇a, [(H⁻¹ ⊗ H⁻¹) : D³φ]⟩
t_as3 = np.einsum('i,jkl,ij,kl', cp.grad_amp, cp.phase_d3, h_inv, h_inv)
# 3. Quartic phase anharmonicity:
# −(1/8) Σᵢⱼₖₗ (H⁻¹)ᵢⱼ (H⁻¹)ₖₗ D⁴φᵢⱼₖₗ
t_s4 = np.einsum('ijkl,ij,kl', cp.phase_d4, h_inv, h_inv)
# 4. Squared cubic anharmonicity (theta-graph Feynman diagram):
# +(5/24) Σᵢⱼₖₗₘₙ (H⁻¹)ᵢₗ (H⁻¹)ⱼₘ (H⁻¹)ₖₙ D³φᵢⱼₖ D³φₗₘₙ
t_s3s3 = np.einsum('ijk,lmn,il,jm,kn', cp.phase_d3, cp.phase_d3, h_inv, h_inv, h_inv)
# Combine correction terms (all real for real-valued a and φ).
correction_factor = term_lap_a - 0.5 * t_as3 - 0.125 * t_s4 + (5.0 / 24.0) * t_s3s3
correction_val = term0 * inv_lam * correction_factor
total_value = term0 + correction_val
return AsymptoticContribution(
leading_term=complex(term0),
correction_term=complex(correction_val),
total_value=complex(total_value),
point=cp,
order_leading=n / 2.0,
method=IntegralMethod.LAPLACE,
)
# --- Saddle-Point Evaluator ---
[docs]
class SaddlePointEvaluator:
"""
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)
"""
def __init__(self, tolerance: float = 1e-8):
self.tolerance = tolerance
# ------------------------------------------------------------------
# Saddle-point search in ℂⁿ
# ------------------------------------------------------------------
[docs]
def find_saddle_points(
self,
analyzer: 'Analyzer',
initial_guesses: List[np.ndarray],
) -> List[np.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
-------
list of complex ndarray
Unique saddle points found (deduplicated within 1e-6).
Each array has shape (dim,) and dtype complex128.
"""
dim = analyzer.dim
func_grad = analyzer.func_grad
def objective(uv: np.ndarray) -> float:
"""Real objective: |∇φ(u + iv)|²."""
z = uv[:dim] + 1j * uv[dim:]
try:
g = np.array(func_grad(*z), dtype=complex)
return float(np.real(np.dot(g.conj(), g)))
except Exception:
return 1e30
saddle_points: List[np.ndarray] = []
for guess in initial_guesses:
# Initial point: real guess, zero imaginary part
uv0 = np.concatenate([np.real(guess), np.zeros(dim)])
try:
res = minimize(objective, uv0, method='L-BFGS-B',
tol=self.tolerance,
options={'maxiter': 2000, 'ftol': self.tolerance**2})
if res.fun < self.tolerance:
z_c = res.x[:dim] + 1j * res.x[dim:]
# Deduplicate: reject if too close to an existing saddle
if not any(np.linalg.norm(z_c - s) < 1e-6 for s in saddle_points):
saddle_points.append(z_c)
except Exception:
pass
return saddle_points
# ------------------------------------------------------------------
# Asymptotic formula at a single saddle point
# ------------------------------------------------------------------
[docs]
def evaluate(self, 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
-------
AsymptoticContribution:
- leading_term : saddle-point formula value.
- correction_term : 0j (not implemented for saddle points).
- order_leading : n/2.
- method : IntegralMethod.SADDLE_POINT.
"""
warnings.warn(
"SaddlePointEvaluator: contour validity (Picard-Lefschetz) is NOT "
"checked. The contribution is correct only if the integration contour "
"can be deformed through this saddle point. Verify independently.",
RuntimeWarning,
stacklevel=2,
)
dim = len(cp.position)
# Degenerate saddle: det H ≈ 0 → formula undefined
if abs(cp.hessian_det) < self.tolerance:
warnings.warn(
f"SaddlePointEvaluator: det(∇²φ) ≈ 0 at saddle {cp.position}. "
"Degenerate saddle points are not supported; returning zero contribution.",
RuntimeWarning,
stacklevel=2,
)
return AsymptoticContribution(
leading_term=0j,
correction_term=0j,
total_value=0j,
point=cp,
order_leading=dim / 2.0,
method=IntegralMethod.SADDLE_POINT,
)
# (2π/λ)^(n/2)
prefactor = (2.0 * np.pi / lam) ** (dim / 2.0)
# exp(iλφ(z_c)) — complex phase, possibly with exponential decay
phase_osc = np.exp(1j * lam * cp.phase_value)
# 1/√det(∇²φ(z_c)) — complex square root, principal branch.
# This encodes both the Gaussian curvature and the Maslov-like phase
# at the complex saddle point.
sqrt_det = np.sqrt(complex(cp.hessian_det)) # principal branch
leading = prefactor * cp.amplitude_value * phase_osc / sqrt_det
return AsymptoticContribution(
leading_term=leading,
correction_term=0j,
total_value=leading,
point=cp,
order_leading=dim / 2.0,
method=IntegralMethod.SADDLE_POINT,
)
# --- Unified Façade ---
[docs]
class AsymptoticEvaluator:
"""
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)
Attributes
----------
sp_evaluator : StationaryPhaseEvaluator
laplace_evaluator : LaplaceEvaluator
saddle_evaluator : SaddlePointEvaluator
"""
def __init__(self, tolerance: float = 1e-8):
self.sp_evaluator = StationaryPhaseEvaluator(tolerance=tolerance)
self.laplace_evaluator = LaplaceEvaluator()
self.saddle_evaluator = SaddlePointEvaluator(tolerance=tolerance)
[docs]
def evaluate(self, 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.
Returns
-------
AsymptoticContribution
Raises
------
ValueError
If cp.method is None, AUTO, or unrecognised.
"""
if cp.method is None or cp.method == IntegralMethod.AUTO:
raise ValueError(
"CriticalPoint.method is None or AUTO. "
"Make sure the point was produced by Analyzer "
"after method resolution (AUTO should have been replaced by a "
"concrete method during __init__)."
)
if cp.method == IntegralMethod.STATIONARY_PHASE:
return self.sp_evaluator.evaluate(cp, lam)
elif cp.method == IntegralMethod.LAPLACE:
return self.laplace_evaluator.evaluate(cp, lam)
elif cp.method == IntegralMethod.SADDLE_POINT:
return self.saddle_evaluator.evaluate(cp, lam)
else:
raise ValueError(f"Unknown IntegralMethod: {cp.method!r}")
[docs]
class AsymptoticVisualizer:
"""
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.
"""
# ------------------------------------------------------------------
# Marker / colour convention for critical-point overlay (shared)
# ------------------------------------------------------------------
_MARKER_STYLE: Dict[SingularityType, Tuple] = {
SingularityType.MORSE : ('o', 'red', 'Morse'),
SingularityType.AIRY_1D : ('*', 'orange', 'Airy'),
SingularityType.AIRY_2D : ('*', 'orange', 'Airy'),
SingularityType.PEARCEY : ('D', 'magenta','Pearcey'),
SingularityType.HIGHER_ORDER: ('s', 'gray', 'Higher-order'),
}
def __init__(self, analyzer: 'Analyzer'):
"""
Parameters
----------
analyzer : Analyzer
Analyzer instance (any method, any dimension).
dim = 2 is required for plot_phase_landscape and plot_integrand.
"""
self.analyzer = analyzer
if analyzer.dim != 2:
warnings.warn(
f"plot_phase_landscape and plot_integrand require dim=2 "
f"(received dim={analyzer.dim}). "
"plot_asymptotic_convergence works for any dimension.",
UserWarning,
)
# ------------------------------------------------------------------
# Internal helpers
# ------------------------------------------------------------------
def _make_grid(
self,
bounds: Tuple[Tuple[float, float], Tuple[float, float]],
n: int,
):
"""Return (X, Y, phi_val, amp_val) on a regular n×n grid."""
x_range = np.linspace(bounds[0][0], bounds[0][1], n)
y_range = np.linspace(bounds[1][0], bounds[1][1], n)
X, Y = np.meshgrid(x_range, y_range)
phi_val = np.asarray(self.analyzer.func_phase(X, Y), dtype=complex)
amp_val = np.asarray(self.analyzer.func_amp(X, Y), dtype=complex)
return X, Y, phi_val, amp_val
def _overlay_critical_points(
self,
ax,
critical_points: List[CriticalPoint],
offset: float = 0.12,
) -> None:
"""Scatter-plot critical / saddle points on an existing Axes."""
plotted_labels: set = set()
for cp in critical_points:
marker, color, label = self._MARKER_STYLE.get(
cp.singularity_type,
('s', 'gray', 'Unknown'),
)
# For complex (saddle) positions, project to real part for 2D display
px = float(np.real(cp.position[0]))
py = float(np.real(cp.position[1]))
kw = dict(c=color, s=150, marker=marker,
edgecolors='white', linewidths=1.5, zorder=10)
if label not in plotted_labels:
ax.scatter(px, py, label=label, **kw)
plotted_labels.add(label)
else:
ax.scatter(px, py, **kw)
ax.text(
px + offset, py + offset,
cp.singularity_type.value,
fontsize=8, color='white', fontweight='bold',
bbox=dict(facecolor='black', alpha=0.55, edgecolor='none', pad=1),
)
@staticmethod
def _contour_panel(
fig, ax, X, Y, Z: np.ndarray,
title: str, label: str,
cmap: str = 'viridis', n_levels: int = 40,
) -> None:
"""Draw a filled contour panel with a colour bar."""
z_min, z_max = np.nanmin(Z), np.nanmax(Z)
if z_min == z_max: # constant field — avoid degenerate levels
z_min -= 1e-10
z_max += 1e-10
levels = np.linspace(z_min, z_max, n_levels)
cf = ax.contourf(X, Y, Z, levels=levels, cmap=cmap, alpha=0.88)
ax.contour(X, Y, Z, levels=levels, colors='k', linewidths=0.35, alpha=0.25)
fig.colorbar(cf, ax=ax, label=label, shrink=0.92)
ax.set_title(title, fontsize=12)
ax.set_xlabel('$x$', fontsize=11)
ax.set_ylabel('$y$', fontsize=11)
ax.set_aspect('equal', adjustable='box')
@staticmethod
def _imshow_panel(
fig, ax, data: np.ndarray,
bounds: Tuple, title: str, label: str,
cmap: str = 'RdBu_r', symmetric: bool = True,
) -> None:
"""Draw an imshow panel with a colour bar."""
vmax = np.nanmax(np.abs(data))
vmin = -vmax if symmetric else 0.0
im = ax.imshow(
data,
extent=[bounds[0][0], bounds[0][1], bounds[1][0], bounds[1][1]],
origin='lower', cmap=cmap,
vmin=vmin, vmax=vmax,
interpolation='bilinear',
)
fig.colorbar(im, ax=ax, label=label, shrink=0.92)
ax.set_title(title, fontsize=12)
ax.set_xlabel('$x$', fontsize=11)
ax.set_ylabel('$y$', fontsize=11)
# ------------------------------------------------------------------
# Public plots
# ------------------------------------------------------------------
[docs]
def plot_phase_landscape(
self,
critical_points: List[CriticalPoint],
bounds: Tuple[Tuple[float, float], Tuple[float, float]] = ((-3, 3), (-3, 3)),
points_per_axis: int = 120,
display: Optional[str] = None,
) -> None:
"""
Visualize the phase function φ(x,y) with critical/saddle-point overlay.
The panels shown depend on the integration method (or on ``display``):
.. list-table::
:widths: 20 80
:header-rows: 1
* - 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
"""
if self.analyzer.dim != 2:
warnings.warn("plot_phase_landscape requires dim=2.", UserWarning)
return
X, Y, phi_val, _ = self._make_grid(bounds, points_per_axis)
# Determine which panels to draw
method = self.analyzer.method
if display is None:
if method == IntegralMethod.STATIONARY_PHASE:
display = 'real'
elif method == IntegralMethod.LAPLACE:
display = 'imag'
else: # SADDLE_POINT or unknown
display = 'both'
panel_map = {
'real': [(np.real(phi_val), r'$\operatorname{Re}(\phi)$',
r'$\operatorname{Re}(\phi(x,y))$', 'viridis')],
'imag': [(np.imag(phi_val), r'$\operatorname{Im}(\phi) = \psi$',
r'$\operatorname{Im}(\phi(x,y))$', 'plasma')],
'abs' : [(np.abs(phi_val), r'$|\phi|$',
r'$|\phi(x,y)|$', 'magma')],
'arg' : [(np.angle(phi_val),r'$\arg(\phi)$',
r'$\arg(\phi(x,y))$', 'hsv')],
'both': [(np.real(phi_val), r'$\operatorname{Re}(\phi)$',
r'$\operatorname{Re}(\phi(x,y))$', 'viridis'),
(np.imag(phi_val), r'$\operatorname{Im}(\phi) = \psi$',
r'$\operatorname{Im}(\phi(x,y))$', 'plasma')],
}
if display not in panel_map:
raise ValueError(f"display must be one of {list(panel_map)}; got {display!r}")
panels = panel_map[display]
ncols = len(panels)
fig, axes = plt.subplots(1, ncols, figsize=(7 * ncols, 6),
squeeze=False)
for ax, (data, cbar_label, title, cmap) in zip(axes[0], panels):
self._contour_panel(fig, ax, X, Y, data, title, cbar_label, cmap)
self._overlay_critical_points(ax, critical_points)
ax.legend(loc='upper right', fontsize=9, framealpha=0.85)
method_label = method.value.replace('_', ' ').title()
fig.suptitle(
f'Phase Landscape — {method_label}',
fontsize=14, fontweight='bold', y=0.98,
)
plt.tight_layout()
plt.show()
[docs]
def plot_integrand(
self,
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:
.. list-table::
:widths: 20 80
:header-rows: 1
* - 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.
"""
if self.analyzer.dim != 2:
warnings.warn("plot_integrand requires dim=2.", UserWarning)
return
X, Y, phi_val, amp_val = self._make_grid(bounds, points_per_axis)
method = self.analyzer.method
# Full complex integrand f = a · exp(iλφ)
f = amp_val * np.exp(1j * lam_value * phi_val)
re_f = np.real(f)
abs_f = np.abs(f)
if method == IntegralMethod.STATIONARY_PHASE:
# Pure oscillation: |f| = |a| is independent of λ; show Re(f).
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
self._imshow_panel(
fig, ax, re_f, bounds,
title=rf'$\operatorname{{Re}}[a\, e^{{i\lambda\phi}}]$'
rf' — Stationary Phase ($\lambda={lam_value}$)',
label=r'$\operatorname{Re}[f]$',
cmap='RdBu_r', symmetric=True,
)
elif method == IntegralMethod.LAPLACE:
# Real exponential: f = a·exp(-λψ); always real for real a and ψ.
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
self._imshow_panel(
fig, ax, re_f, bounds,
title=rf'$a\, e^{{-\lambda\psi}}$'
rf' — Laplace ($\lambda={lam_value}$)',
label=r'$a\,e^{-\lambda\psi}$',
cmap='hot', symmetric=False,
)
else: # SADDLE_POINT: show oscillation + damping envelope side-by-side
fig, (ax_re, ax_abs) = plt.subplots(1, 2, figsize=(15, 6))
self._imshow_panel(
fig, ax_re, re_f, bounds,
title=rf'$\operatorname{{Re}}[a\, e^{{i\lambda\phi}}]$'
rf' ($\lambda={lam_value}$)',
label=r'$\operatorname{Re}[f]$',
cmap='RdBu_r', symmetric=True,
)
self._imshow_panel(
fig, ax_abs, abs_f, bounds,
title=rf'$|a\, e^{{i\lambda\phi}}| = |a|\,e^{{-\lambda\,\operatorname{{Im}}\phi}}$'
rf' ($\lambda={lam_value}$)',
label=r'$|f|$',
cmap='hot', symmetric=False,
)
plt.tight_layout()
plt.show()
[docs]
def plot_asymptotic_convergence(
self,
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:
.. list-table::
:widths: 25 75
:header-rows: 1
* - 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.
"""
# Use the unified evaluator so all three methods are handled.
evaluator = AsymptoticEvaluator()
lams = np.logspace(np.log10(lambda_start), np.log10(lambda_end), num_points)
abs_leading = []
abs_correction = []
for lam in lams:
with warnings.catch_warnings():
# Suppress the Picard-Lefschetz warning during the sweep;
# it has already been emitted when the saddle was found.
warnings.simplefilter('ignore', RuntimeWarning)
res = evaluator.evaluate(cp, lam)
abs_leading.append(np.abs(res.leading_term))
abs_correction.append(np.abs(res.correction_term))
abs_leading = np.array(abs_leading)
abs_correction = np.array(abs_correction)
# Theoretical decay exponent (negative slope on log-log)
dim = len(cp.position)
theoretical_order: Optional[float] = {
SingularityType.MORSE : dim / 2.0,
SingularityType.AIRY_1D : 1.0 / 3.0,
SingularityType.AIRY_2D : 5.0 / 6.0,
SingularityType.PEARCEY : 3.0 / 4.0,
SingularityType.HIGHER_ORDER: None,
}.get(cp.singularity_type, dim / 2.0)
# For Laplace and Saddle-point the Morse formula applies (n/2).
if cp.method in (IntegralMethod.LAPLACE, IntegralMethod.SADDLE_POINT):
theoretical_order = dim / 2.0
fig, ax = plt.subplots(figsize=(9, 6))
ax.loglog(lams, abs_leading, 'o-',
label=r'Leading term $|I_0(\lambda)|$',
linewidth=2.5, markersize=4, alpha=0.85)
# Show correction term only if it is numerically significant
if np.any(abs_correction > 1e-15 * np.max(abs_leading)):
ax.loglog(lams, abs_correction, 's--',
label=r'Correction term $|I_1(\lambda)|$',
linewidth=2.0, markersize=3, alpha=0.80)
# Empirical slope via linear regression on log-log data
valid = abs_leading > 0
if valid.sum() >= 2:
slope_lead = np.polyfit(np.log(lams[valid]),
np.log(abs_leading[valid]), 1)[0]
else:
slope_lead = float('nan')
# Overlay theoretical reference line λ^(-p)
if theoretical_order is not None:
ref = abs_leading[0] * (lams / lams[0]) ** (-theoretical_order)
ax.loglog(lams, ref, 'k:', linewidth=1.4, alpha=0.6,
label=rf'Ref. slope $\lambda^{{-{theoretical_order:.3g}}}$')
# Annotation box with empirical vs theoretical slope
annotation = f'Empirical slope: {slope_lead:.3f}'
if theoretical_order is not None:
annotation += f'\nTheoretical: −{theoretical_order:.3g}'
ax.text(lams[4], abs_leading[4] * 2.0, annotation,
fontsize=9,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.85))
ax.grid(True, which='both', ls=':', alpha=0.6)
ax.legend(loc='best', fontsize=10)
ax.set_xlabel(r'Parameter $\lambda$ (log scale)', fontsize=12)
ax.set_ylabel(r'Magnitude $|I(\lambda)|$ (log scale)', fontsize=12)
# Build a title that works for any dimension
method_label = cp.method.value.replace('_', ' ').title()
pos_str = ', '.join(f'{np.real(v):.3g}' for v in cp.position)
ax.set_title(
f'Asymptotic Convergence — {method_label} / '
f'{cp.singularity_type.value.capitalize()}\n'
rf'$x_c = ({pos_str})$, '
rf'$\phi(x_c) = {cp.phase_value:.4g}$',
fontsize=12,
)
plt.tight_layout()
plt.show()
[docs]
def plot_contribution_decomposition(
critical_points,
analyzer,
lambda_values=None,
show_correction=True,
show_coherent_sum=True,
figsize=None,
):
"""
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 = Symbol('x')
>>> phi = x**4 - x**2 # two minima at x = ±1/√2
>>> amp = 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
... )
See Also
--------
AsymptoticVisualizer.plot_asymptotic_convergence :
Single-critical-point log-log decay plot (simpler but less informative
when multiple critical points are present).
plot_amplitude_decomposition : (in wkb.py)
Spatial WKB analogue — plots amplitude orders aₖ(x) on a grid.
"""
# ------------------------------------------------------------------
# 0. Validate inputs
# ------------------------------------------------------------------
if not critical_points:
raise ValueError("critical_points must contain at least one CriticalPoint.")
if lambda_values is None:
lambda_values = np.logspace(0, 4, 60)
lambda_values = np.asarray(lambda_values, dtype=float)
n_pts = len(critical_points)
evaluator = AsymptoticEvaluator()
# Choose a color for each critical point
colors = [cm.tab10(i % 10) for i in range(n_pts)]
# ------------------------------------------------------------------
# 1. Evaluate all contributions over the λ grid
# results[j][k] = AsymptoticContribution for point j at lambda_values[k]
# ------------------------------------------------------------------
results = []
for cp in critical_points:
row = []
for lam in lambda_values:
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
row.append(evaluator.evaluate(cp, lam))
except Exception:
row.append(None)
results.append(row)
# ------------------------------------------------------------------
# 2. Build figure layout
# ------------------------------------------------------------------
n_panels = 1 + int(show_correction) + int(show_coherent_sum)
if figsize is None:
figsize = (10, 4 * n_panels)
fig, axes_arr = plt.subplots(n_panels, 1, figsize=figsize)
if n_panels == 1:
axes_arr = [axes_arr]
axes_list = list(axes_arr)
panel_idx = 0
ax_lead = axes_list[panel_idx]; panel_idx += 1
ax_corr = axes_list[panel_idx] if show_correction else None
if show_correction:
panel_idx += 1
ax_sum = axes_list[panel_idx] if show_coherent_sum else None
# ------------------------------------------------------------------
# 3. Helper: build a human-readable label for a critical point
# ------------------------------------------------------------------
def _cp_label(j, cp):
pos_str = ", ".join(f"{float(np.real(v)):.3g}" for v in cp.position)
stype = cp.singularity_type.value.capitalize()
method = cp.method.value.replace("_", " ").title()
return f"Point {j+1} x=({pos_str}) [{stype}, {method}]"
# ------------------------------------------------------------------
# 4. Panel 1 — Leading terms
# ------------------------------------------------------------------
reference_plotted = False
for j, (cp, row) in enumerate(zip(critical_points, results)):
lead_vals = np.array([
abs(r.leading_term) if r is not None else np.nan for r in row
])
ax_lead.loglog(lambda_values, lead_vals,
color=colors[j], linewidth=2.0, marker="o",
markersize=3, label=_cp_label(j, cp))
# Overlay one theoretical reference line (from the first valid result)
if not reference_plotted:
valid_rows = [r for r in row if r is not None]
if valid_rows:
p = valid_rows[0].order_leading
ref = lead_vals[np.isfinite(lead_vals)]
if len(ref):
idx0 = np.where(np.isfinite(lead_vals))[0][0]
ref_curve = (lead_vals[idx0]
* (lambda_values / lambda_values[idx0]) ** (-p))
ax_lead.loglog(lambda_values, ref_curve,
"k:", linewidth=1.3, alpha=0.55,
label=rf"Ref. slope $\lambda^{{-{p:.3g}}}$")
reference_plotted = True
ax_lead.set_xlabel(r"$\lambda$", fontsize=12)
ax_lead.set_ylabel(r"$|I_0^{(j)}(\lambda)|$", fontsize=12)
ax_lead.set_title("Leading-term contributions per critical point", fontsize=12)
ax_lead.legend(fontsize=8, loc="best")
ax_lead.grid(True, which="both", ls=":", alpha=0.5)
# ------------------------------------------------------------------
# 5. Panel 2 — Correction terms (optional)
# ------------------------------------------------------------------
if show_correction and ax_corr is not None:
any_correction = False
for j, (cp, row) in enumerate(zip(critical_points, results)):
corr_vals = np.array([
abs(r.correction_term) if r is not None else np.nan for r in row
])
# Check whether correction is non-trivially zero
lead_vals = np.array([
abs(r.leading_term) if r is not None else np.nan for r in row
])
finite_corr = corr_vals[np.isfinite(corr_vals)]
finite_lead = lead_vals[np.isfinite(lead_vals)]
threshold = (1e-14 * np.max(finite_lead)
if len(finite_lead) else 1e-14)
if len(finite_corr) and np.max(finite_corr) > threshold:
ax_corr.loglog(lambda_values, corr_vals,
color=colors[j], linewidth=2.0,
linestyle="--", marker="s", markersize=3,
label=_cp_label(j, cp))
any_correction = True
else:
# Degenerate singularity — correction not computed
ax_corr.plot([], [],
color=colors[j], linestyle="--",
label=_cp_label(j, cp) + " [correction N/A]")
ax_corr.set_xlabel(r"$\lambda$", fontsize=12)
ax_corr.set_ylabel(r"$|I_1^{(j)}(\lambda)|$", fontsize=12)
ax_corr.set_title("Correction-term contributions per critical point",
fontsize=12)
ax_corr.legend(fontsize=8, loc="best")
ax_corr.grid(True, which="both", ls=":", alpha=0.5)
if not any_correction:
ax_corr.text(0.5, 0.5,
"No correction terms available\n(degenerate singularities)",
transform=ax_corr.transAxes,
ha="center", va="center", fontsize=11,
bbox=dict(boxstyle="round", facecolor="lightyellow",
alpha=0.8))
# ------------------------------------------------------------------
# 6. Panel 3 — Coherent sum (optional)
# ------------------------------------------------------------------
if show_coherent_sum and ax_sum is not None:
coherent = np.zeros(len(lambda_values), dtype=complex)
incoherent = np.zeros(len(lambda_values), dtype=float)
for j, row in enumerate(results):
for k, r in enumerate(row):
if r is not None:
coherent[k] += r.total_value
incoherent[k] += abs(r.total_value)
ax_sum.loglog(lambda_values, np.abs(coherent),
"b-", linewidth=2.5, label=r"Coherent sum $|\sum_j I^{(j)}|$")
ax_sum.loglog(lambda_values, incoherent,
"r--", linewidth=1.8,
label=r"Incoherent bound $\sum_j |I^{(j)}|$")
ax_sum.set_xlabel(r"$\lambda$", fontsize=12)
ax_sum.set_ylabel(r"$|I(\lambda)|$", fontsize=12)
ax_sum.set_title("Total asymptotic approximation (all critical points)",
fontsize=12)
ax_sum.legend(fontsize=10, loc="best")
ax_sum.grid(True, which="both", ls=":", alpha=0.5)
# ------------------------------------------------------------------
# 7. Overall title
# ------------------------------------------------------------------
method_label = analyzer.method.value.replace("_", " ").title()
dim_label = f"{analyzer.dim}D"
fig.suptitle(
f"Asymptotic contribution decomposition — {method_label} / {dim_label}"
f" ({n_pts} critical point{'s' if n_pts > 1 else ''})",
fontsize=13, fontweight="bold", y=1.01,
)
plt.tight_layout()
return fig, axes_list
# Backward-compatible alias so that existing code using StationaryPhaseVisualizer
# continues to work without modification.
StationaryPhaseVisualizer = AsymptoticVisualizer
# --- Execution Example ---
if __name__ == "__main__":
x, y = symbols('x y')
# =========================================================================
# Helper: shared report printer (works for all three methods)
# =========================================================================
def print_report(label: str, points, analyzer,
lam_values=(10, 100, 1000),
saddle_points=None):
"""
Print an asymptotic analysis report for the first critical/saddle point.
For SADDLE_POINT, pass the complex coordinates via `saddle_points`
(produced by SaddlePointEvaluator.find_saddle_points).
"""
evaluator = AsymptoticEvaluator()
print(f"\n{'='*60}")
print(f" {label}")
print(f" Detected method: {analyzer.method.value}")
print(f"{'='*60}")
# Select coordinate list: real points or complex saddle points
coords = saddle_points if saddle_points is not None else points
if not coords:
print(" No critical/saddle points found.")
return
cp = analyzer.analyze_point(coords[0])
print(f" Critical point : {cp.position}")
print(f" Singularity : {cp.singularity_type.value}")
print(f" Hessian det : {cp.hessian_det:.4g}")
for lam in lam_values:
try:
res = evaluator.evaluate(cp, lam)
except RuntimeWarning:
pass # saddle-point warning already emitted; continue
ratio = (np.abs(res.correction_term) / np.abs(res.leading_term)
if np.abs(res.leading_term) > 0 else float('nan'))
print(f"\n λ = {lam}")
print(f" Leading term : {res.leading_term:.4e}")
print(f" Correction term : {res.correction_term:.4e}")
print(f" Total : {res.total_value:.4e}")
print(f" |correction/lead| : {ratio:.2%}")
# =========================================================================
# 1. φ purely real → AUTO detects STATIONARY_PHASE
# I(λ) = ∫∫ (1 + x²) exp(iλ (x²/2 + y²/2 + x³/10)) dx dy
# =========================================================================
phi_real = x**2/2 + y**2/2 + x**3/10
amp_real = 1 + x**2
analyzer_sp = Analyzer(phi_real, amp_real, [x, y])
pts_sp = analyzer_find_critical_points([np.array([0., 0.])])
print_report("Case 1 — φ real (AUTO → STATIONARY_PHASE)", pts_sp, analyzer_sp)
# =========================================================================
# 2. φ purely imaginary → AUTO detects LAPLACE
# φ = i·ψ with ψ = x²/2 + y²/2 + x³/20 (ψ real, minimum at origin)
# The integral becomes ∫∫ exp(-λ ψ(x,y)) dx dy
# =========================================================================
phi_imag = I * (x**2/2 + y**2/2 + x**3/20)
amp_imag = Integer(1)
analyzer_lap = Analyzer(phi_imag, amp_imag, [x, y])
pts_lap = analyzer_lap.find_critical_points([np.array([0., 0.])])
print_report("Case 2 — φ imaginary (AUTO → LAPLACE)", pts_lap, analyzer_lap)
# =========================================================================
# 3. φ genuinely complex → AUTO detects SADDLE_POINT
# φ = (x²/2 + y²/2) + i·(x²/4 + y²/4)
# = (1/2 + i/4)(x² + y²)
# Saddle point at origin; contribution is a complex Gaussian.
# =========================================================================
phi_cplx = (x**2/2 + y**2/2) + I*(x**2/4 + y**2/4)
amp_cplx = Integer(1)
analyzer_sdl = Analyzer(phi_cplx, amp_cplx, [x, y])
# For SADDLE_POINT we use SaddlePointEvaluator to search in ℂⁿ
sdl_eval = SaddlePointEvaluator()
saddles = sdl_eval.find_saddle_points(analyzer_sdl, [np.array([0., 0.])])
print_report("Case 3 — φ complex (AUTO → SADDLE_POINT)", [],
analyzer_sdl, saddle_points=saddles)
# =========================================================================
# 4. Visualisations — one visualizer per case
# =========================================================================
print("\n--- Generating Visualizations ---")
bounds2d = ((-2, 2), (-2, 2))
# Case 1: STATIONARY_PHASE
if pts_sp:
cp_sp = analyzer_analyze_point(pts_sp[0])
viz_sp = AsymptoticVisualizer(analyzer_sp)
viz_plot_phase_landscape([cp_sp], bounds=bounds2d)
viz_plot_integrand(lam_value=50, bounds=bounds2d)
viz_plot_asymptotic_convergence(cp_sp)
# Case 2: LAPLACE
if pts_lap:
cp_lap = analyzer_lap.analyze_point(pts_lap[0])
viz_lap = AsymptoticVisualizer(analyzer_lap)
viz_lap.plot_phase_landscape([cp_lap], bounds=bounds2d)
viz_lap.plot_integrand(lam_value=10, bounds=bounds2d)
viz_lap.plot_asymptotic_convergence(cp_lap)
# Case 3: SADDLE_POINT
if saddles:
cp_sdl = analyzer_sdl.analyze_point(saddles[0])
viz_sdl = AsymptoticVisualizer(analyzer_sdl)
viz_sdl.plot_phase_landscape([cp_sdl], bounds=bounds2d)
viz_sdl.plot_integrand(lam_value=20, bounds=bounds2d)
viz_sdl.plot_asymptotic_convergence(cp_sdl)