Elliptic Integrals Example

Demonstrates complete and incomplete elliptic integrals using specialfunctions.ellipk(), specialfunctions.ellipe(), specialfunctions.ellipkinc(), and specialfunctions.ellipeinc(). The application section computes the arc length of an ellipse with semi-axes \(a=2\), \(b=1\) via \(4a\,E(e^2)\).

Example Code

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

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import ellipk, ellipe, ellipkinc, ellipeinc

# ---------------------------------------------------------------------------
# Known values
# ---------------------------------------------------------------------------
print(f"K(0)    = {ellipk(0.0):.10f}  (expected: {np.pi / 2:.10f}  = π/2)")
print(f"K(0.5)  = {ellipk(0.5):.10f}")
print(f"E(0)    = {ellipe(0.0):.10f}  (expected: {np.pi / 2:.10f}  = π/2)")
print(f"E(0.5)  = {ellipe(0.5):.10f}")

# Incomplete integrals at φ = π/2 must recover complete values
phi_half = np.pi / 2
k_test = 0.5
print(f"F(π/2, 0.5) = {ellipkinc(phi_half, k_test):.10f}  (should equal K(0.5))")
print(f"E(π/2, 0.5) = {ellipeinc(phi_half, k_test):.10f}  (should equal E(0.5))")

# ---------------------------------------------------------------------------
# Arc length of an ellipse  a=2, b=1  via E(k)
# Perimeter ≈ 4a · E(e²)  where eccentricity e = √(1 - (b/a)²)
# ---------------------------------------------------------------------------
a, b = 2.0, 1.0
e2 = 1.0 - (b / a) ** 2          # k parameter (= e²)
arc_length = 4.0 * a * ellipe(e2)
print(f"\nEllipse a=2, b=1:")
print(f"  e² (k)      = {e2:.6f}")
print(f"  Arc length  = {arc_length:.6f}  (Ramanujan approx ≈ 9.6884)")

# ---------------------------------------------------------------------------
# Plots
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle("Elliptic Integrals", fontsize=14)

# --- Panel 1: K(k) and E(k) vs k in [0, 0.99] ---
k_vals = np.linspace(0.0, 0.99, 400)
K_vals = ellipk(k_vals)
E_vals = ellipe(k_vals)

axes[0].plot(k_vals, K_vals, "b-", lw=1.8, label=r"$K(k)$")
axes[0].plot(k_vals, E_vals, "r-", lw=1.8, label=r"$E(k)$")
axes[0].axhline(np.pi / 2, color="gray", ls="--", lw=0.8, label=r"$\pi/2$")
axes[0].set_xlim(0, 1)
axes[0].set_ylim(0, 6)
axes[0].set_title("Complete Elliptic Integrals")
axes[0].set_xlabel("k")
axes[0].set_ylabel("Value")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# --- Panel 2: Arc length of ellipse vs eccentricity ---
ecc_vals = np.linspace(0.0, 0.999, 300)          # eccentricity e
k_ecc = ecc_vals ** 2                             # k = e²
# Avoid k=1 (complete integral diverges in K but E(1)=1)
arc_vals = 4.0 * a * ellipe(np.clip(k_ecc, 0, 0.999))

axes[1].plot(ecc_vals, arc_vals, "g-", lw=1.8)
axes[1].axhline(2 * np.pi * a, color="gray", ls="--", lw=0.8, label=f"Circle (a={a})")
axes[1].set_title(r"Arc Length of Ellipse ($a=2$) vs Eccentricity")
axes[1].set_xlabel("Eccentricity e")
axes[1].set_ylabel("Perimeter")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

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

Plot Output

Left: \(K(k)\) and \(E(k)\) as functions of the modulus k. Both equal \(\pi/2\) at \(k=0\); \(K(k)\) diverges as \(k\to 1\) while \(E(k)\to 1\). Right: arc length of the ellipse \(a=2\) as eccentricity increases from 0 (circle, perimeter \(2\pi a\)) to 1 (degenerate line).

Console Output

K(0)    = 1.5707963268  (expected: 1.5707963268  = π/2)
K(0.5)  = 1.8540746773
E(0)    = 1.5707963268  (expected: 1.5707963268  = π/2)
E(0.5)  = 1.3506438810
F(π/2, 0.5) = 1.8540746773  (should equal K(0.5))
E(π/2, 0.5) = 1.3506438810  (should equal E(0.5))

Ellipse a=2, b=1:
  e┬▓ (k)      = 0.750000
  Arc length  = 9.688448  (Ramanujan approx Γëê 9.6884)

Saved: test_output/example_imsl_elliptic_integrals.svg