Exponential and Trigonometric Integrals Example

Demonstrates exponential and trigonometric integral functions using specialfunctions.e1(), specialfunctions.ei(), specialfunctions.si(), and specialfunctions.ci(). The application section computes \(\mathrm{Si}(\pi) \approx 1.8519\), which governs the Gibbs phenomenon: the maximum of the Fourier partial sum of a unit step function converges to \(0.5 + \mathrm{Si}(\pi)/\pi \approx 1.0895\).

Example Code

"""Example: IMSL Special Functions — Exponential and Trigonometric Integrals."""
from __future__ import annotations

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import e1, ei, si, ci

# ---------------------------------------------------------------------------
# Known values
# ---------------------------------------------------------------------------
si_pi = si(np.pi)
ci_1  = ci(1.0)
e1_1  = e1(1.0)

print(f"Si(π)  = {si_pi:.10f}  (expected ≈ 1.8519272584)")
print(f"Ci(1)  = {ci_1:.10f}  (expected ≈ 0.3374039229)")
print(f"E1(1)  = {e1_1:.10f}  (expected ≈ 0.2193839344)")

# ---------------------------------------------------------------------------
# Application: Si(π) in signal processing
#   The normalised sinc function sinc(x) = sin(πx)/(πx) has the property
#   that its integral from -∞ to ∞ equals 1.  Si(π) appears in the partial
#   sums of Fourier series and the Gibbs phenomenon overshoot ≈ 1 + Si(π)/π.
# ---------------------------------------------------------------------------
gibbs_overshoot = 0.5 + si_pi / np.pi
print(f"\nGibbs max of unit-step Fourier sum = 0.5 + Si(π)/π = {gibbs_overshoot:.6f}  (≈ 1.0895)")

# ---------------------------------------------------------------------------
# Plots
# ---------------------------------------------------------------------------
x = np.linspace(0.1, 10.0, 600)
E1_vals = e1(x)
Ei_vals = ei(x)
Si_vals = si(x)
Ci_vals = ci(x)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle("Exponential and Trigonometric Integrals", fontsize=14)

# --- Panel 1: E1 and Ei ---
axes[0].plot(x, E1_vals, "b-",  lw=1.8, label=r"$E_1(x)$")
axes[0].plot(x, np.clip(Ei_vals, -5, 10), "r-", lw=1.8, label=r"$\mathrm{Ei}(x)$")
axes[0].axhline(0, color="k", lw=0.5, ls="--")
axes[0].set_xlim(0.1, 10)
axes[0].set_ylim(-5, 10)
axes[0].set_title(r"$E_1(x)$ and $\mathrm{Ei}(x)$")
axes[0].set_xlabel("x")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# --- Panel 2: Si and Ci ---
axes[1].plot(x, Si_vals, "g-",  lw=1.8, label=r"$\mathrm{Si}(x)$")
axes[1].plot(x, Ci_vals, "m-",  lw=1.8, label=r"$\mathrm{Ci}(x)$")
axes[1].axhline(np.pi / 2, color="g", ls=":", lw=0.8, label=r"$\pi/2$")
axes[1].axhline(0, color="k", lw=0.5, ls="--")
axes[1].axvline(np.pi, color="gray", ls=":", lw=0.8)
axes[1].annotate(r"$x=\pi$", xy=(np.pi + 0.1, -1.5), fontsize=9, color="gray")
axes[1].set_xlim(0.1, 10)
axes[1].set_ylim(-2, 2.5)
axes[1].set_title(r"$\mathrm{Si}(x)$ and $\mathrm{Ci}(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_exponential_integrals.svg", bbox_inches="tight")
print("\nSaved: test_output/example_imsl_exponential_integrals.svg")

Plot Output

Left: \(E_1(x)\) (decays monotonically) and \(\mathrm{Ei}(x)\) (grows, with a logarithmic singularity as \(x\to 0^+\)). Right: \(\mathrm{Si}(x)\) converges to \(\pi/2\); the vertical dotted line marks \(x=\pi\) where \(\mathrm{Si}(\pi)\) appears in the Gibbs formula.

Console Output

Si(π)  = 1.8519370520  (expected ≈ 1.8519272584)
Ci(1)  = 0.3374039229  (expected Γëê 0.3374039229)
E1(1)  = 0.2193839344  (expected Γëê 0.2193839344)

Gibbs max of unit-step Fourier sum = 0.5 + Si(π)/π = 1.089490  (≈ 1.0895)

Saved: test_output/example_imsl_exponential_integrals.svg