Airy Functions Example

Demonstrates Airy functions and their derivatives using specialfunctions.airy_ai(), specialfunctions.airy_bi(), specialfunctions.airy_ai_prime(), and specialfunctions.airy_bi_prime(). The first three zeros of \(\mathrm{Ai}(x)\) are marked on the plot, and the script illustrates the quantum-tunnelling connection formula role of \(\mathrm{Ai}(0)\).

Example Code

"""Example: IMSL Special Functions — Airy Functions."""
from __future__ import annotations

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import airy_ai, airy_bi, airy_ai_prime, airy_bi_prime

# ---------------------------------------------------------------------------
# First zeros of Ai(x)  (well-known tabulated values)
# ---------------------------------------------------------------------------
AI_ZEROS = [-2.338107410459767, -4.087949444130970, -5.520559828095551]
print("First 3 zeros of Ai(x):")
for i, z in enumerate(AI_ZEROS, 1):
    print(f"  a_{i} = {z:.10f}   Ai(a_{i}) = {airy_ai(z):.2e}")

# ---------------------------------------------------------------------------
# Quantum tunnelling connection formula demo
#   For a linear potential V(x) = F·x, the WKB solution switches from
#   oscillatory (x<0) to exponentially decaying (x>0) at x=0, mediated by
#   Ai.  We just print Ai(0) = 3^(-2/3)/Γ(2/3).
# ---------------------------------------------------------------------------
ai0 = airy_ai(0.0)
bi0 = airy_bi(0.0)
print(f"\nAi(0) = {ai0:.10f}")
print(f"Bi(0) = {bi0:.10f}")
print(f"Ai'(0) = {airy_ai_prime(0.0):.10f}")
print(f"Bi'(0) = {airy_bi_prime(0.0):.10f}")

# ---------------------------------------------------------------------------
# Plots
# ---------------------------------------------------------------------------
x = np.linspace(-8.0, 4.0, 800)
Ai  = airy_ai(x)
Bi  = airy_bi(x)
DAi = airy_ai_prime(x)
DBi = airy_bi_prime(x)

# Clip Bi for x > 0 to prevent the exponential blow-up from dominating
Bi_plot  = np.where((x > 0) & (np.abs(Bi)  > 3), np.nan, Bi)
DBi_plot = np.where((x > 0) & (np.abs(DBi) > 4), np.nan, DBi)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle("Airy Functions", fontsize=14)

# --- Panel 1: Ai and Bi ---
axes[0].plot(x, Ai,      "b-", lw=1.8, label=r"$\mathrm{Ai}(x)$")
axes[0].plot(x, Bi_plot, "r-", lw=1.8, label=r"$\mathrm{Bi}(x)$ (clipped $x>0$)")
axes[0].axhline(0, color="k", lw=0.5, ls="--")
for z in AI_ZEROS:
    axes[0].axvline(z, color="gray", ls=":", lw=0.8)
    axes[0].annotate(f"{z:.2f}", xy=(z, 0.05), fontsize=7,
                     ha="center", color="gray")
axes[0].set_xlim(-8, 4)
axes[0].set_ylim(-1.0, 1.5)
axes[0].set_title(r"$\mathrm{Ai}(x)$ and $\mathrm{Bi}(x)$")
axes[0].set_xlabel("x")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# --- Panel 2: Derivatives ---
axes[1].plot(x, DAi,      "b--", lw=1.8, label=r"$\mathrm{Ai}'(x)$")
axes[1].plot(x, DBi_plot, "r--", lw=1.8, label=r"$\mathrm{Bi}'(x)$ (clipped $x>0$)")
axes[1].axhline(0, color="k", lw=0.5, ls="--")
axes[1].set_xlim(-8, 4)
axes[1].set_ylim(-1.5, 1.5)
axes[1].set_title(r"Airy Derivatives $\mathrm{Ai}'(x)$, $\mathrm{Bi}'(x)$")
axes[1].set_xlabel("x")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("test_output/example_imsl_airy_functions.svg", bbox_inches="tight")
print("\nSaved: test_output/example_imsl_airy_functions.svg")

Plot Output

Left: \(\mathrm{Ai}(x)\) (blue) and \(\mathrm{Bi}(x)\) (red, clipped for \(x>0\) to prevent exponential blow-up). Dotted vertical lines mark the first three zeros of \(\mathrm{Ai}\). Right: derivatives \(\mathrm{Ai}'(x)\) and \(\mathrm{Bi}'(x)\).

Console Output

First 3 zeros of Ai(x):
  a_1 = -2.3381074105   Ai(a_1) = -4.27e-17
  a_2 = -4.0879494441   Ai(a_2) = -1.34e-16
  a_3 = -5.5205598281   Ai(a_3) = 5.14e-16

Ai(0) = 0.3550280539
Bi(0) = 0.6149266274
Ai'(0) = -0.2588194038
Bi'(0) = 0.4482883574

Saved: test_output/example_imsl_airy_functions.svg