Modified Bessel Functions Example

Demonstrates bessel_i0(), bessel_i1(), bessel_in(), bessel_k0(), bessel_k1(), bessel_kn(), and bessel_yn().

Example Code

"""Example: IMSL Special Functions — Modified Bessel functions (I, K) and Bessel Y."""
from __future__ import annotations

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import (
    bessel_i0, bessel_i1, bessel_in,
    bessel_k0, bessel_k1, bessel_kn,
    bessel_yn,
)

# ---------------------------------------------------------------------------
# Known-value checks
# ---------------------------------------------------------------------------
print(f"bessel_i0(0)  = {bessel_i0(0.0):.10f}  (expected: 1.0000000000)")
print(f"bessel_i1(0)  = {bessel_i1(0.0):.10f}  (expected: 0.0000000000)")
print(f"bessel_i0(1)  = {bessel_i0(1.0):.10f}")
print(f"bessel_k0(1)  = {bessel_k0(1.0):.10f}")
print(f"bessel_k0(0.01) = {bessel_k0(0.01):.6f}  (large — K0 diverges as x→0+)")
print(f"bessel_yn(0, 1) = {bessel_yn(0, 1.0):.10f}")
print(f"bessel_yn(1, 1) = {bessel_yn(1, 1.0):.10f}")

# ---------------------------------------------------------------------------
# Plotting
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.suptitle("Modified Bessel Functions and Bessel Y", fontsize=14)

# --- Panel 1: Modified Bessel I_n(x) for n = 0, 1, 2 ----------------------
ax = axes[0]
x_i = np.linspace(0, 5, 400)
ax.plot(x_i, bessel_i0(x_i), "tab:blue",  lw=1.8, label="I₀(x)")
ax.plot(x_i, bessel_i1(x_i), "tab:orange", lw=1.8, label="I₁(x)")
ax.plot(x_i, bessel_in(2, x_i), "tab:green", lw=1.8, label="I₂(x)")
ax.set_title("Modified Bessel I — first kind\nI₀(0)=1, monotone increasing")
ax.set_xlabel("x")
ax.set_ylabel("Iₙ(x)")
ax.legend()
ax.set_ylim(0, 30)
ax.grid(True, alpha=0.3)

# --- Panel 2: Modified Bessel K_n(x) for n = 0, 1, 2 ----------------------
ax = axes[1]
x_k = np.linspace(0.05, 5, 400)
ax.plot(x_k, bessel_k0(x_k), "tab:blue",  lw=1.8, label="K₀(x)")
ax.plot(x_k, bessel_k1(x_k), "tab:orange", lw=1.8, label="K₁(x)")
ax.plot(x_k, bessel_kn(2, x_k), "tab:green", lw=1.8, label="K₂(x)")
ax.set_title("Modified Bessel K — second kind\nK₀→∞ as x→0⁺, decays for large x")
ax.set_xlabel("x")
ax.set_ylabel("Kₙ(x)")
ax.legend()
ax.set_ylim(0, 10)
ax.grid(True, alpha=0.3)

# --- Panel 3: Bessel Y_n(x) for n = 0, 1, 2 --------------------------------
ax = axes[2]
x_y = np.linspace(0.1, 12, 600)
for n, colour in [(0, "tab:blue"), (1, "tab:orange"), (2, "tab:green")]:
    y = bessel_yn(n, x_y)
    ax.plot(x_y, y, color=colour, lw=1.8, label=f"Y{n}(x)")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_title("Bessel Y — second kind\nYₙ(x) diverges as x→0⁺")
ax.set_xlabel("x")
ax.set_ylabel("Yₙ(x)")
ax.legend()
ax.set_ylim(-3, 1)
ax.grid(True, alpha=0.3)

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

# ---------------------------------------------------------------------------
# Engineering application note: heat conduction in a cylinder
# ---------------------------------------------------------------------------
print()
print("Engineering note — heat conduction in a hollow cylinder:")
print("  Temperature profile: T(r) ∝ I0(α·r)  for the interior")
print("  Heat flux decays outward ∝ K0(α·r)  for the exterior")
r_vals = np.array([0.5, 1.0, 1.5, 2.0])
alpha = 1.0
for r in r_vals:
    print(f"  I0({alpha}·{r}) = {bessel_i0(alpha * r):.6f}   "
          f"K0({alpha}·{r}) = {bessel_k0(alpha * r):.6f}")

Plot Output

Modified Bessel functions of the first kind Iₙ(x), second kind Kₙ(x), and Bessel functions of the second kind Yₙ(x).

Console Output

bessel_i0(0)  = 1.0000000000  (expected: 1.0000000000)
bessel_i1(0)  = 0.0000000000  (expected: 0.0000000000)
bessel_i0(1)  = 1.2660658778
bessel_k0(1)  = 0.4210244382
bessel_k0(0.01) = 4.721245  (large — K0 diverges as x→0+)
bessel_yn(0, 1) = 0.0882569642
bessel_yn(1, 1) = -0.7812128213
Saved: test_output/example_imsl_bessel_modified.svg

Engineering note ΓÇö heat conduction in a hollow cylinder:
  Temperature profile: T(r) ∝ I0(α·r)  for the interior
  Heat flux decays outward ∝ K0(α·r)  for the exterior
  I0(1.0┬╖0.5) = 1.063483   K0(1.0┬╖0.5) = 0.924419
  I0(1.0┬╖1.0) = 1.266066   K0(1.0┬╖1.0) = 0.421024
  I0(1.0┬╖1.5) = 1.646723   K0(1.0┬╖1.5) = 0.213806
  I0(1.0┬╖2.0) = 2.279585   K0(1.0┬╖2.0) = 0.113894