Inverse Laplace Transform Example — Three Classical Pairs

This example demonstrates numerical inversion of three classical Laplace transform pairs using transforms.inverse_laplace() and transforms.inverse_laplace_smooth() (Talbot’s fixed-contour method):

• \(F(s) = \tfrac{1}{s}\) → \(f(t) = 1\) (unit step)

• \(F(s) = \tfrac{1}{s+2}\) → \(f(t) = e^{-2t}\) (decaying exponential)

• \(F(s) = \tfrac{1}{s^2+1}\) → \(f(t) = \sin(t)\) (harmonic oscillation)

Maximum absolute errors relative to the analytical solutions are printed to the console, typically on the order of \(10^{-11}\).

Example Code

"""IMSL Inverse Laplace Transform example: three classical Laplace pairs.

Demonstrates numerical inversion of three well-known Laplace pairs using
:func:`transforms.inverse_laplace` and :func:`transforms.inverse_laplace_smooth`:

- F(s) = 1/s         → f(t) = 1            (unit step)
- F(s) = 1/(s+2)     → f(t) = exp(-2t)     (decaying exponential, a=2)
- F(s) = 1/(s²+1)    → f(t) = sin(t)       (harmonic oscillation)

Outputs:
- Table printed to stdout showing max absolute error for each pair
- SVG plot saved to test_output/example_imsl_laplace.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Callable, Dict

import matplotlib.pyplot as plt
import numpy as np

from transforms import inverse_laplace, inverse_laplace_smooth


def run_demo_imsl_laplace() -> Dict[str, object]:
    """Run IMSL inverse Laplace transform example.

    Numerically inverts three classical Laplace-domain functions and compares
    the results against their known analytical time-domain solutions.
    Uses :func:`transforms.inverse_laplace` for the first two cases and
    :func:`transforms.inverse_laplace_smooth` for the third.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``max_errors`` (dict mapping
            case name to float) and ``plot_path`` (str).
    """
    t = np.linspace(0.05, 5.0, 200)

    # --- Case 1: F(s) = 1/s  →  f(t) = 1 (unit step) ---
    def F_step(s: complex) -> complex:
        """Laplace transform of the unit step: 1/s.

        Args:
            s (complex): Complex Laplace variable.

        Returns:
            complex: Value of F(s) = 1/s.
        """
        return 1.0 / s

    f_step_num = inverse_laplace(F_step, t)
    f_step_ana = np.ones_like(t)
    err_step = float(np.max(np.abs(f_step_num - f_step_ana)))

    # --- Case 2: F(s) = 1/(s+2)  →  f(t) = exp(-2t) ---
    a = 2.0

    def F_exp(s: complex) -> complex:
        """Laplace transform of exp(-at): 1/(s+a), a=2.

        Args:
            s (complex): Complex Laplace variable.

        Returns:
            complex: Value of F(s) = 1/(s+a).
        """
        return 1.0 / (s + a)

    f_exp_num = inverse_laplace(F_exp, t)
    f_exp_ana = np.exp(-a * t)
    err_exp = float(np.max(np.abs(f_exp_num - f_exp_ana)))

    # --- Case 3: F(s) = 1/(s²+1)  →  f(t) = sin(t) ---
    def F_sin(s: complex) -> complex:
        """Laplace transform of sin(t): 1/(s²+1).

        Args:
            s (complex): Complex Laplace variable.

        Returns:
            complex: Value of F(s) = 1/(s^2+1).
        """
        return 1.0 / (s ** 2 + 1.0)

    f_sin_num = inverse_laplace_smooth(F_sin, t)
    f_sin_ana = np.sin(t)
    err_sin = float(np.max(np.abs(f_sin_num - f_sin_ana)))

    max_errors = {
        "1/s → 1": err_step,
        "1/(s+2) → exp(-2t)": err_exp,
        "1/(s²+1) → sin(t)": err_sin,
    }

    print("\nIMSL Inverse Laplace Transform Example")
    print(f"  Time range: t ∈ [{t[0]:.2f}, {t[-1]:.2f}],  N={len(t)} points")
    print("-" * 50)
    print(f"  {'Case':<28} {'Max |error|':>12}")
    print("-" * 50)
    for case, err in max_errors.items():
        print(f"  {case:<28} {err:>12.2e}")
    print("-" * 50)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_laplace.svg"

    fig, axes = plt.subplots(1, 3, figsize=(14, 5))

    _cases: list[tuple[str, np.ndarray, np.ndarray, str, str]] = [
        ("F(s)=1/s\nf(t)=1", f_step_num, f_step_ana, "#6d28d9", f"Max err={err_step:.2e}"),
        ("F(s)=1/(s+2)\nf(t)=exp(−2t)", f_exp_num, f_exp_ana, "#0369a1", f"Max err={err_exp:.2e}"),
        ("F(s)=1/(s²+1)\nf(t)=sin(t)", f_sin_num, f_sin_ana, "#059669", f"Max err={err_sin:.2e}"),
    ]

    for ax, (title, numerical, analytical, color, err_label) in zip(axes, _cases):
        ax.plot(t, analytical, color="#94a3b8", linewidth=2.5, linestyle="-",
                label="Analytical", zorder=1)
        ax.plot(t, numerical, color=color, linewidth=1.5, linestyle="--",
                label=f"Numerical\n{err_label}", zorder=2)
        ax.set_xlabel("t")
        ax.set_ylabel("f(t)")
        ax.set_title(title)
        ax.legend(fontsize=8)
        ax.grid(True, alpha=0.3)

    fig.suptitle("Inverse Laplace Transform: Numerical vs Analytical", fontsize=13)
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"max_errors": max_errors, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_laplace()

Plot Output

Three-panel comparison of numerical (dashed) vs. analytical (solid) inverse Laplace transform results for the unit step, decaying exponential, and harmonic sine functions. Maximum errors are of order \(10^{-11}\).

Console Output

IMSL Inverse Laplace Transform Example
  Time range: t Γêê [0.05, 5.00],  N=200 points
--------------------------------------------------
  Case                          Max |error|
--------------------------------------------------
  1/s → 1                          2.21e-11
  1/(s+2) → exp(-2t)               2.02e-11
  1/(s²+1) → sin(t)                7.01e-12
--------------------------------------------------