"""
ch07_composite_trapezoidal_simpson.py

Composite trapezoidal and Simpson rules for numerical integration.

The program evaluates

    integral sin(x) dx,   0 <= x <= pi,

for successively finer uniform grids. The exact value is 2.

The example illustrates the different convergence rates of the
composite trapezoidal and Simpson rules as the number of subintervals
is doubled.

The Simpson routine is based on code written by Alexander Godunov,
October 2009. Python version prepared for the companion website, 2026.

The numerical routines use only the Python standard library.
Matplotlib is used only for visualization.
"""

import math
import matplotlib.pyplot as plt


def f(x):
    """Test function for the integration example."""
    return math.sin(x)


def composite_trapezoidal(fun, a, b, n):
    """Composite trapezoidal rule on [a,b]."""
    if n < 1:
        raise ValueError("composite_trapezoidal: n must be positive.")

    h = (b - a) / n
    total = 0.5 * (fun(a) + fun(b))

    for i in range(1, n):
        total += fun(a + i * h)

    return h * total


def composite_simpson(fun, a, b, n):
    """Composite Simpson rule on [a,b]; n must be even."""
    if n < 2 or n % 2 != 0:
        raise ValueError(
            "composite_simpson: n must be a positive even integer."
        )

    h = (b - a) / n

    sum_odd = 0.0
    for i in range(1, n, 2):
        sum_odd += fun(a + i * h)

    sum_even = 0.0
    for i in range(2, n - 1, 2):
        sum_even += fun(a + i * h)

    return h * (fun(a) + fun(b) + 4.0 * sum_odd + 2.0 * sum_even) / 3.0


def plot_integrand(a, b):
    """Plot the integrand before carrying out the numerical integration."""
    nplot = 1000
    xplot = [a + i * (b - a) / nplot for i in range(nplot + 1)]
    yplot = [f(x) for x in xplot]

    plt.figure()
    plt.plot(xplot, yplot, linewidth=1.5)
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title("Integrand for composite trapezoidal and Simpson rules")
    plt.grid(True)


def main():
    a = 0.0
    b = math.pi
    exact = 2.0

    plot_integrand(a, b)

    print("Composite trapezoidal and Simpson rules")
    print("Integral of sin(x) from 0 to pi")
    print(f"Exact value          = {exact:15.7e}\n")

    print("        n       Trapezoidal          Simpson")
    print("------------------------------------------------")

    n = 2
    for _ in range(16):
        trap = composite_trapezoidal(f, a, b, n)
        simp = composite_simpson(f, a, b, n)
        print(f"{n:9d}{trap:19.7e}{simp:19.7e}")
        n *= 2

    plt.show()


if __name__ == "__main__":
    main()