integration.tex

\documentclass[a4paper, 11pt]{report}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{xcolor}
\usepackage[makestderr]{pythontex}

\definecolor{trap}{RGB}{231, 76, 60}
\definecolor{simp}{RGB}{46, 204, 113}
\definecolor{gauss}{RGB}{52, 152, 219}
\definecolor{romberg}{RGB}{155, 89, 182}

\title{Numerical Integration Methods:\\
Quadrature Algorithms and Error Analysis}
\author{Department of Computational Mathematics\\Technical Report NM-2024-001}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This report presents a comprehensive analysis of numerical integration (quadrature) methods. We implement and compare the trapezoidal rule, Simpson's rule, Romberg integration, and Gaussian quadrature. Error analysis demonstrates convergence rates, and we verify results against analytically known integrals. All computations use PythonTeX for reproducibility.
\end{abstract}

\tableofcontents

\chapter{Introduction}

Numerical integration approximates definite integrals when analytical solutions are unavailable or impractical. We seek to compute:
\begin{equation}
I = \int_a^b f(x) \, dx
\end{equation}

\section{Newton-Cotes Formulas}
These methods interpolate $f(x)$ with polynomials at equally spaced nodes:
\begin{itemize}
    \item Trapezoidal rule: Linear interpolation
    \item Simpson's rule: Quadratic interpolation
    \item Higher-order rules: Cubic, quartic, etc.
\end{itemize}

\chapter{Trapezoidal Rule}

\section{Formula}
The composite trapezoidal rule with $n$ subintervals:
\begin{equation}
T_n = h\left[\frac{f(a) + f(b)}{2} + \sum_{i=1}^{n-1} f(a + ih)\right], \quad h = \frac{b-a}{n}
\end{equation}

Error: $E_T = -\frac{(b-a)h^2}{12}f''(\xi) = O(h^2)$

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.special import roots_legendre
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

def trapezoidal(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h * (0.5*y[0] + np.sum(y[1:-1]) + 0.5*y[-1])

def simpson(f, a, b, n):
    if n % 2 == 1:
        n += 1
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1])

def romberg(f, a, b, max_iter=10, tol=1e-10):
    R = np.zeros((max_iter, max_iter))
    h = b - a
    R[0, 0] = 0.5 * h * (f(a) + f(b))

    for i in range(1, max_iter):
        h = h / 2
        sum_new = sum(f(a + (2*k-1)*h) for k in range(1, 2**(i-1) + 1))
        R[i, 0] = 0.5 * R[i-1, 0] + h * sum_new

        for j in range(1, i + 1):
            R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1)

        if i > 0 and abs(R[i, i] - R[i-1, i-1]) < tol:
            return R[i, i], i, R[:i+1, :i+1]

    return R[max_iter-1, max_iter-1], max_iter, R

def gauss_legendre(f, a, b, n):
    nodes, weights = roots_legendre(n)
    # Transform from [-1, 1] to [a, b]
    x = 0.5 * (b - a) * nodes + 0.5 * (a + b)
    return 0.5 * (b - a) * np.sum(weights * f(x))
\end{pycode}

\section{Visualization}

\begin{pycode}
# Test function
def f_test(x):
    return np.exp(-x**2)

a, b = 0, 2
exact = 0.5 * np.sqrt(np.pi) * (1 - 2/np.sqrt(np.pi) * integrate.quad(lambda t: np.exp(-t**2), 2, np.inf)[0])
exact_value = integrate.quad(f_test, a, b)[0]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Trapezoidal visualization
ax = axes[0, 0]
x_plot = np.linspace(a, b, 200)
ax.plot(x_plot, f_test(x_plot), 'b-', linewidth=2, label='$f(x) = e^{-x^2}$')
n_trap = 6
x_trap = np.linspace(a, b, n_trap + 1)
y_trap = f_test(x_trap)
for i in range(n_trap):
    ax.fill([x_trap[i], x_trap[i], x_trap[i+1], x_trap[i+1]],
            [0, y_trap[i], y_trap[i+1], 0], alpha=0.3, color='red')
ax.scatter(x_trap, y_trap, color='red', s=50, zorder=5)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(f'Trapezoidal Rule ($n = {n_trap}$)')
ax.legend()
ax.grid(True, alpha=0.3)

# Simpson visualization
ax = axes[0, 1]
ax.plot(x_plot, f_test(x_plot), 'b-', linewidth=2, label='$f(x) = e^{-x^2}$')
n_simp = 4
h_simp = (b - a) / n_simp
for i in range(0, n_simp, 2):
    x0, x1, x2 = a + i*h_simp, a + (i+1)*h_simp, a + (i+2)*h_simp
    y0, y1, y2 = f_test(x0), f_test(x1), f_test(x2)
    # Quadratic fit
    coeffs = np.polyfit([x0, x1, x2], [y0, y1, y2], 2)
    x_para = np.linspace(x0, x2, 50)
    y_para = np.polyval(coeffs, x_para)
    ax.fill_between(x_para, 0, y_para, alpha=0.3, color='green')
x_simp = np.linspace(a, b, n_simp + 1)
ax.scatter(x_simp, f_test(x_simp), color='green', s=50, zorder=5)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(f"Simpson's Rule ($n = {n_simp}$)")
ax.legend()
ax.grid(True, alpha=0.3)

# Convergence comparison
ax = axes[1, 0]
n_values = np.array([4, 8, 16, 32, 64, 128, 256])
trap_errors = []
simp_errors = []
gauss_errors = []

for n in n_values:
    trap_errors.append(abs(trapezoidal(f_test, a, b, n) - exact_value))
    simp_errors.append(abs(simpson(f_test, a, b, n) - exact_value))
    if n <= 20:
        gauss_errors.append(abs(gauss_legendre(f_test, a, b, n) - exact_value))

ax.loglog(n_values, trap_errors, 'ro-', label='Trapezoidal $O(h^2)$')
ax.loglog(n_values, simp_errors, 'gs-', label="Simpson's $O(h^4)$")
ax.loglog(n_values[:len(gauss_errors)], gauss_errors, 'b^-', label='Gauss-Legendre')
ax.set_xlabel('Number of subintervals $n$')
ax.set_ylabel('Absolute Error')
ax.set_title('Convergence Comparison')
ax.legend()
ax.grid(True, alpha=0.3)

# Romberg tableau
ax = axes[1, 1]
result, iters, R = romberg(f_test, a, b)
# Show Romberg tableau as heatmap
R_display = R.copy()
R_display[R_display == 0] = np.nan
im = ax.imshow(np.abs(R_display - exact_value), cmap='viridis_r', aspect='auto')
ax.set_xlabel('Richardson Extrapolation Level $j$')
ax.set_ylabel('Refinement Level $i$')
ax.set_title('Romberg Tableau (Error)')
plt.colorbar(im, ax=ax, label='$|R_{i,j} - I|$')

plt.tight_layout()
plt.savefig('integration_methods.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{integration_methods.pdf}
\caption{Numerical integration methods: (a) trapezoidal approximation, (b) Simpson's parabolic approximation, (c) error convergence, (d) Romberg tableau.}
\end{figure}

\chapter{Simpson's Rule}

\section{Formula}
Composite Simpson's 1/3 rule (requires even $n$):
\begin{equation}
S_n = \frac{h}{3}\left[f(a) + 4\sum_{i=1,3,5,...}^{n-1} f(a+ih) + 2\sum_{i=2,4,6,...}^{n-2} f(a+ih) + f(b)\right]
\end{equation}

Error: $E_S = -\frac{(b-a)h^4}{180}f^{(4)}(\xi) = O(h^4)$

\chapter{Romberg Integration}

\section{Richardson Extrapolation}
Romberg integration uses the trapezoidal rule with Richardson extrapolation:
\begin{equation}
R_{i,j} = R_{i,j-1} + \frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}
\end{equation}

\begin{pycode}
# Test multiple integrals
test_functions = [
    (lambda x: np.exp(-x**2), 0, 2, "e^{-x^2}"),
    (lambda x: np.sin(x)/x if x != 0 else 1.0, 0.001, np.pi, "\\sin(x)/x"),
    (lambda x: 1/(1 + x**2), 0, 1, "1/(1+x^2)"),
]

results_table = []
for func, a, b, name in test_functions:
    f_vec = np.vectorize(func)
    exact = integrate.quad(f_vec, a, b)[0]

    trap_32 = trapezoidal(f_vec, a, b, 32)
    simp_32 = simpson(f_vec, a, b, 32)
    gauss_8 = gauss_legendre(f_vec, a, b, 8)
    romb_val, _, _ = romberg(f_vec, a, b)

    results_table.append({
        'name': name,
        'exact': exact,
        'trap_err': abs(trap_32 - exact),
        'simp_err': abs(simp_32 - exact),
        'gauss_err': abs(gauss_8 - exact),
        'romb_err': abs(romb_val - exact)
    })
\end{pycode}

\begin{table}[htbp]
\centering
\caption{Error comparison for different integrals}
\begin{tabular}{@{}lcccc@{}}
\toprule
Function & Trapezoidal & Simpson & Gauss-8 & Romberg \\
\midrule
$\py{results_table[0]['name']}$ & \py{f"{results_table[0]['trap_err']:.2e}"} & \py{f"{results_table[0]['simp_err']:.2e}"} & \py{f"{results_table[0]['gauss_err']:.2e}"} & \py{f"{results_table[0]['romb_err']:.2e}"} \\
$\py{results_table[1]['name']}$ & \py{f"{results_table[1]['trap_err']:.2e}"} & \py{f"{results_table[1]['simp_err']:.2e}"} & \py{f"{results_table[1]['gauss_err']:.2e}"} & \py{f"{results_table[1]['romb_err']:.2e}"} \\
$\py{results_table[2]['name']}$ & \py{f"{results_table[2]['trap_err']:.2e}"} & \py{f"{results_table[2]['simp_err']:.2e}"} & \py{f"{results_table[2]['gauss_err']:.2e}"} & \py{f"{results_table[2]['romb_err']:.2e}"} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{Gaussian Quadrature}

\section{Theory}
Gaussian quadrature chooses nodes and weights optimally:
\begin{equation}
\int_{-1}^{1} f(x) \, dx \approx \sum_{i=1}^{n} w_i f(x_i)
\end{equation}

Nodes $x_i$ are roots of Legendre polynomials. An $n$-point formula is exact for polynomials of degree $\leq 2n-1$.

\begin{pycode}
# Gaussian quadrature analysis
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Gauss nodes and weights
ax = axes[0]
for n in [2, 3, 4, 5]:
    nodes, weights = roots_legendre(n)
    ax.scatter(nodes, [n]*len(nodes), s=weights*500, alpha=0.7, label=f'$n = {n}$')
ax.set_xlabel('Node position $x_i$')
ax.set_ylabel('Number of points $n$')
ax.set_title('Gauss-Legendre Nodes')
ax.set_xlim(-1.1, 1.1)
ax.legend()
ax.grid(True, alpha=0.3)

# Polynomial exactness
ax = axes[1]
n_points = range(2, 11)
degrees_exact = []
for n in n_points:
    # Test polynomial of degree 2n-1
    deg = 2*n - 1
    poly = lambda x, d=deg: x**d
    exact_poly = 2 / (deg + 1) if deg % 2 == 0 else 0
    gauss_result = gauss_legendre(poly, -1, 1, n)
    error = abs(gauss_result - exact_poly)
    degrees_exact.append(error < 1e-10)
ax.bar(n_points, [2*n-1 for n in n_points], color='blue', alpha=0.7)
ax.set_xlabel('Number of points $n$')
ax.set_ylabel('Exact polynomial degree')
ax.set_title('Polynomial Exactness: $2n-1$')
ax.grid(True, alpha=0.3)

# Error vs number of Gauss points
ax = axes[2]
n_gauss = range(2, 16)
gauss_errors_exp = [abs(gauss_legendre(f_test, a, b, n) - exact_value) for n in n_gauss]
ax.semilogy(n_gauss, gauss_errors_exp, 'b^-')
ax.set_xlabel('Number of Gauss points $n$')
ax.set_ylabel('Absolute Error')
ax.set_title('Gauss-Legendre Convergence')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('gaussian_quadrature.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{gaussian_quadrature.pdf}
\caption{Gaussian quadrature: node positions (left), polynomial exactness (center), convergence for $e^{-x^2}$ (right).}
\end{figure}

\chapter{Adaptive Quadrature}

\begin{pycode}
# Adaptive Simpson's rule
def adaptive_simpson(f, a, b, tol=1e-8, max_depth=50):
    def simpson_recursive(a, b, fa, fm, fb, S, tol, depth):
        c = (a + b) / 2
        d = (a + c) / 2
        e = (c + b) / 2
        fd = f(d)
        fe = f(e)
        h = b - a

        S_left = h/12 * (fa + 4*fd + fm)
        S_right = h/12 * (fm + 4*fe + fb)
        S_new = S_left + S_right

        if depth >= max_depth or abs(S_new - S) < 15*tol:
            return S_new + (S_new - S)/15

        return (simpson_recursive(a, c, fa, fd, fm, S_left, tol/2, depth+1) +
                simpson_recursive(c, b, fm, fe, fb, S_right, tol/2, depth+1))

    c = (a + b) / 2
    fa, fm, fb = f(a), f(c), f(b)
    S = (b - a)/6 * (fa + 4*fm + fb)
    return simpson_recursive(a, b, fa, fm, fb, S, tol, 0)

# Test on oscillatory function
def f_osc(x):
    return np.sin(10*x) * np.exp(-x)

exact_osc = integrate.quad(f_osc, 0, 2*np.pi)[0]
adaptive_result = adaptive_simpson(f_osc, 0, 2*np.pi)
fixed_result = simpson(f_osc, 0, 2*np.pi, 100)

fig, ax = plt.subplots(figsize=(10, 4))
x = np.linspace(0, 2*np.pi, 500)
ax.plot(x, f_osc(x), 'b-', linewidth=1)
ax.fill_between(x, 0, f_osc(x), alpha=0.3)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title('$f(x) = \\sin(10x) e^{-x}$')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('adaptive_function.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{adaptive_function.pdf}
\caption{Oscillatory test function for adaptive quadrature.}
\end{figure}

Adaptive Simpson error: $\py{f'{abs(adaptive_result - exact_osc):.2e}'}$ \\
Fixed Simpson (n=100) error: $\py{f'{abs(fixed_result - exact_osc):.2e}'}$

\chapter{Summary}

\begin{pycode}
summary_data = [
    ('Trapezoidal', '$O(h^2)$', 'Simple, low accuracy'),
    ("Simpson's", '$O(h^4)$', 'Good accuracy, even $n$'),
    ('Romberg', '$O(h^{2k})$', 'High accuracy, extrapolation'),
    ('Gauss-Legendre', 'Exponential', 'Optimal for smooth functions'),
]
\end{pycode}

\begin{table}[htbp]
\centering
\caption{Summary of quadrature methods}
\begin{tabular}{@{}lll@{}}
\toprule
Method & Convergence & Notes \\
\midrule
\py{summary_data[0][0]} & \py{summary_data[0][1]} & \py{summary_data[0][2]} \\
\py{summary_data[1][0]} & \py{summary_data[1][1]} & \py{summary_data[1][2]} \\
\py{summary_data[2][0]} & \py{summary_data[2][1]} & \py{summary_data[2][2]} \\
\py{summary_data[3][0]} & \py{summary_data[3][1]} & \py{summary_data[3][2]} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{Conclusions}

\begin{enumerate}
    \item Trapezoidal rule: $O(h^2)$, simple but low accuracy
    \item Simpson's rule: $O(h^4)$, good balance of simplicity and accuracy
    \item Romberg: achieves high accuracy through Richardson extrapolation
    \item Gaussian quadrature: optimal for smooth functions, exponential convergence
    \item Adaptive methods: efficient for oscillatory or peaked integrands
\end{enumerate}

\end{document}