More Orthogonal Polynomials & Inverse Error Functions Example

Demonstrates chebyshev_u(), laguerre_poly(), jacobi_poly(), erfinv(), and erfcinv().

Example Code

"""Example: IMSL Special Functions — Chebyshev U, Laguerre, Jacobi, erfinv/erfcinv."""
from __future__ import annotations

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import chebyshev_u, laguerre_poly, jacobi_poly, erfinv, erfcinv, erf, erfc

# ---------------------------------------------------------------------------
# Known-value checks
# ---------------------------------------------------------------------------
roundtrip = erfinv(erf(0.5))
print(f"erfinv(erf(0.5))  = {roundtrip:.10f}  (expected: 0.5000000000  — round-trip)")
print(f"erfinv(0)         = {erfinv(0.0):.10f}  (expected: 0.0)")
print(f"erfcinv(1)        = {erfcinv(1.0):.10f}  (expected: 0.0  — erfc(0)=1)")
print(f"erfcinv(erfc(0.7)) = {erfcinv(erfc(0.7)):.10f}  (expected: 0.7000000000)")
print(f"chebyshev_u(0, 0.5) = {chebyshev_u(0, 0.5):.10f}  (expected: 1.0)")
print(f"chebyshev_u(1, 0.5) = {chebyshev_u(1, 0.5):.10f}  (expected: 1.0  = 2*0.5)")
print(f"laguerre_poly(0, 1) = {laguerre_poly(0, 1.0):.10f}  (expected: 1.0)")
print(f"laguerre_poly(1, 1) = {laguerre_poly(1, 1.0):.10f}  (expected: 0.0  = 1-x at x=1)")
print(f"jacobi_poly(0, 0.5, 0.5, 0.5) = {jacobi_poly(0, 0.5, 0.5, 0.5):.10f}  (expected: 1.0)")

# ---------------------------------------------------------------------------
# Plotting
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle("Chebyshev U, Laguerre, Jacobi Polynomials and Inverse Error Functions", fontsize=13)

# --- Panel 1: Chebyshev U polynomials on [-1, 1] ---------------------------
ax = axes[0, 0]
x_cheb = np.linspace(-1, 1, 400)
colours = ["tab:blue", "tab:orange", "tab:green", "tab:red", "tab:purple"]
for n, colour in enumerate(colours):
    ax.plot(x_cheb, chebyshev_u(n, x_cheb), color=colour, lw=1.8, label=f"U{n}(x)")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_title("Chebyshev U polynomials Uₙ(x)")
ax.set_xlabel("x")
ax.set_ylabel("Uₙ(x)")
ax.set_ylim(-5.5, 5.5)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# --- Panel 2: Laguerre polynomials on [0, 10] ------------------------------
ax = axes[0, 1]
x_lag = np.linspace(0, 10, 400)
for n, colour in enumerate(colours):
    y = laguerre_poly(n, x_lag)
    ax.plot(x_lag, y, color=colour, lw=1.8, label=f"L{n}(x)")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_title("Laguerre polynomials Lₙ(x)")
ax.set_xlabel("x")
ax.set_ylabel("Lₙ(x)")
ax.set_ylim(-20, 20)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# --- Panel 3: Jacobi polynomials P(n, alpha=0.5, beta=0.5) on [-1, 1] -----
ax = axes[1, 0]
x_jac = np.linspace(-1, 1, 400)
alpha_j, beta_j = 0.5, 0.5
for n, colour in enumerate(colours[:4]):
    ax.plot(x_jac, jacobi_poly(n, alpha_j, beta_j, x_jac),
            color=colour, lw=1.8, label=f"P{n}^(0.5,0.5)(x)")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_title("Jacobi polynomials Pₙ^(0.5,0.5)(x)")
ax.set_xlabel("x")
ax.set_ylabel("Pₙ^(α,β)(x)")
ax.set_ylim(-3, 3)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# --- Panel 4: erfinv and erfcinv (inverse error functions) -----------------
ax = axes[1, 1]
# erfinv domain: (-1, 1) — plot (-0.99, 0.99)
u_erf = np.linspace(-0.99, 0.99, 400)
ax.plot(u_erf, erfinv(u_erf), "tab:blue", lw=1.8, label="erfinv(u)")

# erfcinv domain: (0, 2) — plot (0.01, 1.99)
u_erfc = np.linspace(0.01, 1.99, 400)
ax.plot(u_erfc - 1, erfcinv(u_erfc), "tab:orange", lw=1.8, ls="--",
        label="erfcinv(1+u)  [shifted to same x-axis]")

ax.axhline(0, color="k", lw=0.5, ls=":")
ax.axvline(0, color="k", lw=0.5, ls=":")
ax.set_title("Inverse Error Functions\nerfinv and erfcinv")
ax.set_xlabel("Argument (erfinv: u ∈ (−1,1); erfcinv: v−1, v ∈ (0,2))")
ax.set_ylabel("Inverse value")
ax.set_ylim(-3, 3)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("test_output/example_imsl_more_polynomials.svg", bbox_inches="tight")
print("Saved: test_output/example_imsl_more_polynomials.svg")

Plot Output

Chebyshev U polynomials Uₙ(x), Laguerre polynomials Lₙ(x), Jacobi polynomials Pₙ^(0.5,0.5)(x), and the inverse error functions erfinv and erfcinv.

Console Output

erfinv(erf(0.5))  = 0.5000000000  (expected: 0.5000000000  ΓÇö round-trip)
erfinv(0)         = 0.0000000000  (expected: 0.0)
erfcinv(1)        = -0.0000000000  (expected: 0.0  ΓÇö erfc(0)=1)
erfcinv(erfc(0.7)) = 0.7000000000  (expected: 0.7000000000)
chebyshev_u(0, 0.5) = 1.0000000000  (expected: 1.0)
chebyshev_u(1, 0.5) = 1.0000000000  (expected: 1.0  = 2*0.5)
laguerre_poly(0, 1) = 1.0000000000  (expected: 1.0)
laguerre_poly(1, 1) = 0.0000000000  (expected: 0.0  = 1-x at x=1)
jacobi_poly(0, 0.5, 0.5, 0.5) = 1.0000000000  (expected: 1.0)
Saved: test_output/example_imsl_more_polynomials.svg