Asymptotic Methods

Many problems in physics cannot be solved exactly but yield to asymptotic analysis — extracting the dominant behavior in some limit (large quantum numbers, small coupling, high frequency, deep tunneling). Asymptotic methods transform impossible exact calculations into powerful approximate ones. This chapter develops the key techniques: asymptotic series, the WKB approximation, the method of steepest descent, and the method of stationary phase.

1. Asymptotic Series

An asymptotic expansion is not a convergent series in the usual sense. Rather, it provides increasingly accurate approximations up to a point, after which the terms grow. Formally:

Definition: Asymptotic Expansion

The series $\sum_{n=0}^{\infty} a_n / x^n$ is an asymptotic expansion of $f(x)$ as$x \to \infty$ if for each fixed $N$:

$$f(x) - \sum_{n=0}^{N} \frac{a_n}{x^n} = O\left(\frac{1}{x^{N+1}}\right) \quad \text{as } x \to \infty$$

We write $f(x) \sim \sum_{n=0}^{\infty} a_n / x^n$.

The classic example is the exponential integral:

$$E_1(x) = \int_x^{\infty} \frac{e^{-t}}{t} dt \sim \frac{e^{-x}}{x}\sum_{n=0}^{\infty} \frac{(-1)^n n!}{x^n}$$

The coefficients $(-1)^n n!$ grow factorially, so the series diverges for every $x$. Yet truncating at the optimal number of terms ($N \approx x$) gives exponentially accurate results:

$$\text{error} \sim e^{-x} \sqrt{\frac{2\pi}{x}} \quad \text{(superasymptotic accuracy)}$$

Why divergent series are useful: In perturbative QFT, the perturbation series in the coupling constant is typically asymptotic (Dyson's argument). Nevertheless, truncating at the optimal order gives the extraordinary precision of QED (anomalous magnetic moment of the electron computed to 10 significant figures). The divergence at high orders signals non-perturbative effects (instantons, tunneling).

Stokes phenomenon. An asymptotic expansion can change form discontinuously across certain rays in the complex plane (Stokes lines). This is deeply connected to the analytic continuation of solutions and the appearance of exponentially small corrections.

2. The WKB Approximation

The WKB (Wentzel-Kramers-Brillouin) approximation is a semiclassical method for solving differential equations when a parameter (typically $\hbar$) is small. We derive it from the time-independent Schrodinger equation:

$$-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi = E\psi$$

Step 1: The WKB ansatz. Write the wave function as:

$$\psi(x) = A(x) \exp\left(\frac{i}{\hbar} S(x)\right)$$

where $S(x)$ is a rapidly varying phase and $A(x)$ is a slowly varying amplitude. Substituting into the Schrodinger equation and separating real and imaginary parts:

Step 2: Expand in powers of $\hbar$. Write $S = S_0 + \hbar S_1 + \hbar^2 S_2 + \cdots$. At leading order ($\hbar^0$):

$$\left(\frac{dS_0}{dx}\right)^2 = 2m(E - V(x)) \equiv p(x)^2$$

This is the Hamilton-Jacobi equation! The solution is $S_0 = \pm\int p(x)\,dx$, where$p(x) = \sqrt{2m(E-V(x))}$ is the classical momentum.

Step 3: Next order ($\hbar^1$). This gives:

$$2S_0' S_1' + S_0'' = 0 \quad \Rightarrow \quad S_1 = -\frac{1}{2}\ln|S_0'| = -\frac{1}{2}\ln|p(x)|$$

Combining, the WKB wave function in the classically allowed region ($E > V$) is:

$$\psi_{\mathrm{WKB}}(x) = \frac{1}{\sqrt{p(x)}}\left[C_+ \exp\left(\frac{i}{\hbar}\int^x p(x')\,dx'\right) + C_- \exp\left(-\frac{i}{\hbar}\int^x p(x')\,dx'\right)\right]$$

In the classically forbidden region ($E < V$), $p(x)$ becomes imaginary,$\kappa(x) = \sqrt{2m(V-E)}$, and the wave function decays exponentially:

$$\psi_{\mathrm{WKB}}(x) = \frac{1}{\sqrt{\kappa(x)}}\left[D_+ \exp\left(\frac{1}{\hbar}\int^x \kappa(x')\,dx'\right) + D_- \exp\left(-\frac{1}{\hbar}\int^x \kappa(x')\,dx'\right)\right]$$

Connection formulas at turning points. The WKB approximation breaks down at classical turning points where $p(x) = 0$. Near a turning point $x_0$, we linearize $V(x) \approx V(x_0) + V'(x_0)(x-x_0)$ and solve exactly in terms of Airy functions. Matching gives the famous connection formulas:

$$\frac{2}{\sqrt{\kappa}}\cos\left(\frac{1}{\hbar}\int_{x_0}^{x} p\,dx' - \frac{\pi}{4}\right) \quad \longleftrightarrow \quad \frac{1}{\sqrt{\kappa}} \exp\left(-\frac{1}{\hbar}\int_{x}^{x_0} \kappa\,dx'\right)$$

Bohr-Sommerfeld quantization. For a bound state between turning points$a$ and $b$, requiring single-valuedness gives:

$$\oint p(x)\,dx = 2\int_a^b p(x)\,dx = 2\pi\hbar\left(n + \frac{1}{2}\right), \quad n = 0, 1, 2, \ldots$$

The $1/2$ is the Maslov correction from the turning points. For the harmonic oscillator$V = \frac{1}{2}m\omega^2 x^2$, this gives the exact energy levels$E_n = \hbar\omega(n + 1/2)$.

3. WKB and Quantum Tunneling

One of the most important applications of the WKB method is computing tunneling rates. For a particle incident on a barrier where $V(x) > E$ between turning points$a$ and $b$, the WKB transmission coefficient is:

$$T \approx \exp\left(-\frac{2}{\hbar}\int_a^b \sqrt{2m(V(x) - E)}\,dx\right)$$

This formula explains the extreme sensitivity of tunneling to barrier width and height. Applications include:

Alpha decay: The Gamow factor explains why nuclear lifetimes span$10^{-7}$ to $10^{17}$ seconds — small changes in the Coulomb barrier produce enormous changes in $T$.
Scanning tunneling microscope: The exponential sensitivity of tunneling current to tip-sample distance gives atomic resolution.
Instantons in QFT: Tunneling between degenerate vacua is described by instanton solutions, with action given by the WKB integral in imaginary time.

Double well and energy splitting. For a symmetric double well, tunneling lifts the degeneracy of the two ground states. The energy splitting is:

$$\Delta E \approx \frac{\hbar\omega}{\pi}\exp\left(-\frac{1}{\hbar}\int_a^b \sqrt{2m(V(x) - E)}\,dx\right)$$

This splitting is the origin of the ammonia maser frequency and the mechanism behind band structure in solids (tight-binding model).

4. Laplace's Method

Laplace's method evaluates integrals of the form:

$$I(\lambda) = \int_a^b f(x) e^{\lambda \phi(x)} dx, \quad \lambda \to \infty$$

The key idea is that as $\lambda \to \infty$, the integrand is dominated by the neighborhood of the maximum of $\phi(x)$. If $\phi$ has a unique maximum at $x_0 \in (a,b)$with $\phi''(x_0) < 0$:

Derivation. Expand $\phi(x)$ about $x_0$:

$$\phi(x) \approx \phi(x_0) + \frac{1}{2}\phi''(x_0)(x - x_0)^2$$

Substituting and extending the limits to $\pm\infty$ (justified since the Gaussian contribution away from $x_0$ is exponentially suppressed):

$$I(\lambda) \approx f(x_0) e^{\lambda\phi(x_0)} \int_{-\infty}^{\infty} e^{\frac{\lambda}{2}\phi''(x_0)(x-x_0)^2} dx = f(x_0) e^{\lambda\phi(x_0)} \sqrt{\frac{2\pi}{\lambda|\phi''(x_0)|}}$$

Classic application: Stirling's formula. Consider$n! = \Gamma(n+1) = \int_0^\infty t^n e^{-t}\,dt$. Writing $t = nx$:

$$n! = n^{n+1}\int_0^\infty e^{n(\ln x - x)}\,dx$$

The exponent $\phi(x) = \ln x - x$ has maximum at $x_0 = 1$ with$\phi(1) = -1$ and $\phi''(1) = -1$. Laplace's method gives:

$$n! \approx n^{n+1} e^{-n} \sqrt{\frac{2\pi}{n}} = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$

This is Stirling's approximation, accurate to within a fraction of a percent even for small $n$.

5. Method of Steepest Descent (Saddle Point Method)

The method of steepest descent generalizes Laplace's method to complex integrals:

$$I(\lambda) = \int_C g(z) e^{\lambda f(z)} dz, \quad \lambda \to \infty$$

where $C$ is a contour in the complex plane. The idea is to deform $C$ to pass through a saddle point $z_0$ (where $f'(z_0) = 0$) along the path of steepest descent of $\mathrm{Re}(f)$.

Why steepest descent? At a saddle point, write $f(z) = f(z_0) + \frac{1}{2}f''(z_0)(z-z_0)^2 + \cdots$. Let $f''(z_0) = |f''(z_0)|e^{i\alpha}$. Along the direction $z - z_0 = re^{i\theta}$:

$$\mathrm{Re}\left[\frac{1}{2}f''(z_0)(z-z_0)^2\right] = \frac{r^2}{2}|f''(z_0)|\cos(\alpha + 2\theta)$$

The steepest descent direction satisfies $\cos(\alpha + 2\theta) = -1$, i.e.,$\theta = (\pi - \alpha)/2$. Along this path, $\mathrm{Im}(f)$ is constant (no oscillations), and $\mathrm{Re}(f)$ decreases most rapidly. The result is:

$$I(\lambda) \approx g(z_0) e^{\lambda f(z_0)} \sqrt{\frac{2\pi}{\lambda |f''(z_0)|}} e^{i\theta_0}$$

where $\theta_0 = (\pi - \alpha)/2$ is the steepest descent angle.

Application to path integrals. In quantum mechanics, the propagator is:

$$K(x_f, x_i; T) = \int \mathcal{D}[x(t)] \, e^{iS[x]/\hbar}$$

The steepest descent approximation corresponds to expanding around the classical path ($\delta S = 0$), giving the semiclassical (Van Vleck) propagator. The Gaussian integral over fluctuations produces the one-loop correction.

6. Method of Stationary Phase

The method of stationary phase applies to oscillatory integrals:

$$I(\lambda) = \int_a^b g(x) e^{i\lambda \phi(x)} dx, \quad \lambda \to \infty$$

For large $\lambda$, rapid oscillations cause cancellation everywhere except near stationary points where $\phi'(x_0) = 0$. This is the mathematical realization of the principle of stationary phase in optics and wave mechanics.

Expanding $\phi(x) \approx \phi(x_0) + \frac{1}{2}\phi''(x_0)(x-x_0)^2$ and evaluating the resulting Fresnel integral:

$$I(\lambda) \approx g(x_0) e^{i\lambda\phi(x_0)} \sqrt{\frac{2\pi}{\lambda|\phi''(x_0)|}} \, e^{i\pi\,\mathrm{sgn}(\phi'')/4}$$

The phase factor $e^{i\pi\,\mathrm{sgn}(\phi'')/4}$ is characteristic of stationary phase and appears in Fresnel diffraction.

Connection to classical limit. The Feynman path integral has the form$\int \mathcal{D}[x] \, e^{iS[x]/\hbar}$. In the classical limit $\hbar \to 0$, stationary phase selects paths where $\delta S = 0$ — precisely the classical trajectories. This is the deepest justification of classical mechanics as the$\hbar \to 0$ limit of quantum mechanics.

Multiple stationary points. If there are several stationary points, each contributes independently to leading order:

$$I(\lambda) \approx \sum_{k} g(x_k) e^{i\lambda\phi(x_k)} \sqrt{\frac{2\pi}{\lambda|\phi''(x_k)|}} e^{i\pi\,\mathrm{sgn}(\phi''(x_k))/4}$$

Interference between contributions from different stationary points gives rise to semiclassical interference patterns (e.g., the Airy pattern near a caustic in optics).

7. Beyond Leading Order: Systematic Corrections

The leading-order results above can be systematically improved. For Laplace's method, expanding$\phi(x)$ and $f(x)$ to higher orders around $x_0$:

$$I(\lambda) = e^{\lambda\phi(x_0)}\sqrt{\frac{2\pi}{\lambda|\phi''(x_0)|}}\left[f(x_0) + \frac{c_1}{\lambda} + \frac{c_2}{\lambda^2} + \cdots\right]$$

where the coefficients $c_k$ involve derivatives of $\phi$ and $f$ at $x_0$. The first correction is:

$$c_1 = \frac{f''(x_0)}{2|\phi''(x_0)|} - \frac{f(x_0)\phi^{(4)}(x_0)}{8[\phi''(x_0)]^2} + \frac{5f(x_0)[\phi'''(x_0)]^2}{24[\phi''(x_0)]^3}$$

For the WKB approximation, higher-order corrections in $\hbar$ are obtained from$S_2, S_3, \ldots$ in the expansion $S = \sum_{n=0}^\infty \hbar^n S_n$. The next correction modifies the quantization condition:

$$\oint p(x)\,dx = 2\pi\hbar\left(n + \frac{1}{2}\right) + \hbar^2 \oint \frac{V''(x)}{48\,p(x)^3}\,dx + \cdots$$

Borel summation provides a way to extract meaningful results from divergent asymptotic series. The Borel transform replaces $n!$ by 1:

$$f(x) = \int_0^\infty e^{-t} B(t/x) \, dt, \quad B(z) = \sum_{n=0}^\infty \frac{a_n}{n!} z^n$$

If $B(z)$ converges and can be analytically continued, the integral recovers $f(x)$. Singularities of $B(z)$ on the positive real axis (renormalons, instantons) signal non-perturbative effects and are a major topic in modern quantum field theory.

8. WKB in Higher Dimensions and the Eikonal Approximation

In higher dimensions, the WKB method generalizes to the eikonal approximation. For the Helmholtz equation $\nabla^2\psi + k^2 n^2(\mathbf{r})\psi = 0$ with slowly varying refractive index $n(\mathbf{r})$, write $\psi = A(\mathbf{r})e^{ik S(\mathbf{r})}$:

$$|\nabla S|^2 = n^2(\mathbf{r}) \quad \text{(eikonal equation)}$$

$$2\nabla S \cdot \nabla A + A\nabla^2 S = 0 \quad \text{(transport equation)}$$

The eikonal equation is the Hamilton-Jacobi equation of ray optics. Its characteristics are the classical rays (geodesics in the medium). The transport equation determines the amplitude variation along rays, giving conservation of energy flux in ray tubes.

Caustics occur where neighboring rays converge and the WKB amplitude diverges. Near a caustic, the eikonal approximation breaks down and must be replaced by uniform asymptotic expansions involving Airy functions (for fold caustics) or Pearcey functions (for cusp caustics). The Maslov index tracks the phase accumulated at each caustic crossing.

9. Comparison and Connections

The asymptotic methods form a unified framework:

Method	Integral Type	Key Idea	Physics
Laplace	$\int e^{\lambda\phi}$ (real)	Maximum of exponent	Stat. mech.
Steepest descent	$\int e^{\lambda f(z)}$ (complex)	Saddle point	Path integrals
Stationary phase	$\int e^{i\lambda\phi}$ (oscillatory)	Phase cancellation	Classical limit
WKB	ODE with small parameter	Slowly varying envelope	Semiclassics

Steepest descent and stationary phase are related by a Wick rotation ($\lambda \to i\lambda$). In physics, this corresponds to the passage between Euclidean and Minkowski signature. Laplace's method is the special case of steepest descent for purely real integrals.

10. Numerical Exploration: Asymptotic Methods

This simulation demonstrates asymptotic series truncation, WKB energy levels, the Laplace method for Stirling's approximation, and quantum tunneling through a barrier.

Python

script.py242 lines

#!/usr/bin/env python3
"""
asymptotic_methods_simulation.py
Numerical demonstrations of asymptotic methods using numpy.
"""
import numpy as np
np.set_printoptions(precision=10, suppress=True, linewidth=100)

# ===== PART 1: Asymptotic Series - Exponential Integral =====
print("="*70)
print("PART 1: Asymptotic Series for E_1(x) = int_x^inf e^{-t}/t dt")
print("="*70)
print()

def E1_numerical(x, N_quad=10000):
    """Compute E_1(x) by numerical integration."""
    # Use substitution t = x + u, integrate u from 0 to large cutoff
    u_max = 50.0 + 10*x  # Ensure convergence
    u = np.linspace(0, u_max, N_quad)
    du = u[1] - u[0]
    integrand = np.exp(-(x + u)) / (x + u)
    return np.sum(integrand) * du

def E1_asymptotic(x, N_terms):
    """Asymptotic series: E_1(x) ~ (e^{-x}/x) * sum_{n=0}^{N} (-1)^n n! / x^n."""
    total = 0.0
    term = 1.0
    for n in range(N_terms):
        total += term
        term *= -(n + 1) / x
    return np.exp(-x) / x * total

x_test = 10.0
E1_exact = E1_numerical(x_test, 100000)
print(f"x = {x_test}")
print(f"E_1(x) numerical = {E1_exact:.12e}")
print()
print(f"{'N terms':>8}  {'Asymptotic':>18}  {'Abs Error':>14}  {'Rel Error':>14}")
print("-" * 65)

best_err = 1e10
best_N = 0
for N in range(1, 25):
    approx = E1_asymptotic(x_test, N)
    err = abs(approx - E1_exact)
    rel_err = err / abs(E1_exact)
    marker = ""
    if err < best_err:
        best_err = err
        best_N = N
    else:
        if N == best_N + 1:
            marker = " <-- divergence begins"
    print(f"{N:>8}  {approx:>18.12e}  {err:>14.6e}  {rel_err:>14.6e}{marker}")

print()
print(f"Optimal truncation at N = {best_N} terms (approx = x = {x_test:.0f})")
print(f"Best error = {best_err:.2e} (superasymptotic accuracy)")
print()

# ===== PART 2: Stirling's Approximation via Laplace's Method =====
print("="*70)
print("PART 2: Stirling's Approximation (Laplace's Method)")
print("="*70)
print()

print(f"{'n':>5}  {'n!':>20}  {'Stirling':>20}  {'Rel Error':>14}")
print("-" * 65)

factorial_exact = 1.0
for n in range(1, 16):
    factorial_exact *= n
    stirling = np.sqrt(2 * np.pi * n) * (n / np.e)**n
    rel_err = abs(stirling - factorial_exact) / factorial_exact
    print(f"{n:>5}  {factorial_exact:>20.1f}  {stirling:>20.4f}  {rel_err:>14.6e}")
print()
print("Stirling is remarkably accurate even for small n!")
print()

# ===== PART 3: WKB Energy Levels =====
print("="*70)
print("PART 3: WKB (Bohr-Sommerfeld) Energy Levels")
print("="*70)
print()

# Harmonic oscillator: V(x) = 0.5 * x^2
# Bohr-Sommerfeld: int_a^b sqrt(2(E - V)) dx = pi*(n + 1/2)  [hbar=m=1]
# Turning points: x = +/- sqrt(2E)
# Integral = pi*E
# So E_n = n + 1/2  (EXACT for harmonic oscillator!)

print("Harmonic Oscillator V(x) = x^2/2:")
print(f"{'n':>5}  {'E_WKB':>12}  {'E_exact':>12}  {'Error':>14}")
print("-" * 50)
for n in range(10):
    E_exact = n + 0.5
    # WKB: integral = pi*E = pi*(n + 1/2), so E = n + 1/2
    E_wkb = n + 0.5
    print(f"{n:>5}  {E_wkb:>12.6f}  {E_exact:>12.6f}  {abs(E_wkb - E_exact):>14.2e}")
print("(WKB is exact for the harmonic oscillator!)")
print()

# Quartic oscillator: V(x) = x^4
# Bohr-Sommerfeld: 2 * int_0^{E^{1/4}} sqrt(2(E - x^4)) dx = pi*(n + 1/2)
print("Quartic Oscillator V(x) = x^4:")
print(f"{'n':>5}  {'E_WKB':>12}  {'E_numerical':>12}  {'Rel Error':>14}")
print("-" * 55)

def quartic_bs_integral(E, N_pts=10000):
    """Compute Bohr-Sommerfeld integral for V=x^4."""
    x_tp = E**0.25  # turning point
    x = np.linspace(0, x_tp, N_pts)
    dx = x[1] - x[0]
    integrand = np.sqrt(np.maximum(2*(E - x**4), 0))
    return 2 * np.sum(integrand) * dx

def solve_quartic_wkb(n, E_range=(0.1, 100)):
    """Find WKB energy for quartic oscillator."""
    target = np.pi * (n + 0.5)
    # Bisection
    E_lo, E_hi = E_range
    for _ in range(100):
        E_mid = (E_lo + E_hi) / 2
        val = quartic_bs_integral(E_mid)
        if val < target:
            E_lo = E_mid
        else:
            E_hi = E_mid
    return (E_lo + E_hi) / 2

# "Exact" energies from matrix diagonalization
def quartic_exact(n_max, basis_size=200):
    """Compute quartic oscillator energies using harmonic oscillator basis."""
    # H = p^2/2 + x^4 in HO basis
    dim = basis_size
    H = np.zeros((dim, dim))
    for i in range(dim):
        # <i|p^2/2|i> = (2i+1)/4 * omega (with omega=1)
        # But we use the full matrix elements of x^4
        pass
    # Direct approach: use DVR (discrete variable representation)
    x_max = 10.0
    N = 500
    dx = 2 * x_max / N
    x_grid = np.linspace(-x_max, x_max, N+1)[1:-1]  # exclude endpoints
    # Kinetic energy matrix (finite difference)
    T = np.zeros((len(x_grid), len(x_grid)))
    for i in range(len(x_grid)):
        T[i, i] = np.pi**2 / (3 * dx**2)
        for j in range(i+1, len(x_grid)):
            T[i, j] = 2 * (-1)**(i-j) / ((i-j)**2 * dx**2)
            T[j, i] = T[i, j]
    V_diag = np.diag(x_grid**4)
    H = 0.5 * T + V_diag
    evals = np.sort(np.linalg.eigvalsh(H))
    return evals[:n_max]

exact_energies = quartic_exact(8)

for n in range(8):
    E_wkb = solve_quartic_wkb(n)
    E_ex = exact_energies[n]
    rel_err = abs(E_wkb - E_ex) / E_ex
    print(f"{n:>5}  {E_wkb:>12.6f}  {E_ex:>12.6f}  {rel_err:>14.6e}")
print()

# ===== PART 4: WKB Tunneling Probability =====
print("="*70)
print("PART 4: Quantum Tunneling (WKB vs Exact)")
print("="*70)
print()

# Square barrier: V(x) = V0 for 0 < x < a, 0 otherwise
# Exact transmission: T = 1 / (1 + V0^2 sinh^2(kappa*a) / (4E(V0-E)))
# WKB: T ~ exp(-2*kappa*a)

V0 = 5.0
a = 2.0
print(f"Square barrier: V0 = {V0}, width a = {a}")
print(f"{'E/V0':>8}  {'T_exact':>14}  {'T_WKB':>14}  {'Ratio':>10}")
print("-" * 55)

for E_frac in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]:
    E = E_frac * V0
    kappa = np.sqrt(2 * (V0 - E))

# Exact
    sinh_val = np.sinh(kappa * a)
    T_exact = 1.0 / (1.0 + V0**2 * sinh_val**2 / (4 * E * (V0 - E)))

# WKB
    T_wkb = np.exp(-2 * kappa * a)

ratio = T_wkb / T_exact if T_exact > 0 else float('inf')
    print(f"{E_frac:>8.1f}  {T_exact:>14.6e}  {T_wkb:>14.6e}  {ratio:>10.4f}")

print()
print("WKB captures the exponential dependence correctly.")
print("The ratio approaches 1 for thick barriers (deep tunneling regime).")
print()

# ===== PART 5: Stationary Phase =====
print("="*70)
print("PART 5: Stationary Phase Approximation")
print("="*70)
print()

# Compute I(lambda) = int_{-inf}^{inf} exp(i*lambda*(x^2/2 - x^4/12)) dx
# Stationary points: phi'(x) = x - x^3/3 = 0 => x = 0, +/-sqrt(3)
# At x=0: phi(0) = 0, phi''(0) = 1

def oscillatory_integral_numerical(lam, N_pts=50000):
    """Compute int exp(i*lam*(x^2/2 - x^4/12)) dx numerically."""
    x_max = 6.0
    x = np.linspace(-x_max, x_max, N_pts)
    dx = x[1] - x[0]
    phase = lam * (x**2 / 2 - x**4 / 12)
    # Add convergence factor for x^4 term
    damping = np.exp(-0.01 * x**4)
    integrand = np.exp(1j * phase) * damping
    return np.sum(integrand) * dx

def stationary_phase_approx(lam):
    """Stationary phase at x=0: phi(0)=0, phi''(0)=1."""
    return np.sqrt(2 * np.pi / lam) * np.exp(1j * np.pi / 4)

print("Oscillatory integral: I(lam) = int exp(i*lam*(x^2/2 - x^4/12)) dx")
print("Stationary phase at x=0: I ~ sqrt(2*pi/lam) * exp(i*pi/4)")
print()
print(f"{'lambda':>8}  {'|I_numerical|':>14}  {'|I_approx|':>14}  {'Rel Error':>14}")
print("-" * 58)

for lam in [5, 10, 20, 50, 100, 200]:
    I_num = oscillatory_integral_numerical(lam, 100000)
    I_sp = stationary_phase_approx(lam)
    rel_err = abs(abs(I_num) - abs(I_sp)) / abs(I_sp)
    print(f"{lam:>8}  {abs(I_num):>14.8f}  {abs(I_sp):>14.8f}  {rel_err:>14.6e}")

print()
print("Stationary phase becomes more accurate as lambda increases.")
print("This is the mathematical basis of the classical limit in quantum mechanics.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Special Functions Variational Calculus →