Numerical Integration Calculator (python, written by Claude Code)

Claude Code wrote this python project from the task below. It installed and ran on a clean Ubuntu 22.04 machine as written.

Task given to the agent:

TASK: Numerical Integration Calculator

Write a program that computes definite integrals of mathematical functions using multiple numerical methods, comparing accuracy and convergence across methods.

FUNCTIONAL REQUIREMENTS:
- Accept a mathematical expression as a command-line argument via --function flag (e.g., --function "sin(x)*exp(-x)")
- Accept integration bounds via --lower and --upper flags
- Support multiple numerical integration methods selectable via --method flag: trapezoidal rule, Simpson's rule, Simpson's 3/8 rule, Gaussian quadrature, and Romberg integration
- Support a configurable number of subintervals via --intervals flag (default: 1000) for methods that use subdivision
- Run all methods and compare results via --compare flag, showing each method's result, estimated error, and computation time
- Support adaptive integration: automatically refine the interval count until the result converges within a specified tolerance via --tolerance flag (default: 1e-10)
- Parse mathematical expressions supporting: basic operators (+, -, *, /, ^), standard functions (sin, cos, tan, exp, log, sqrt, abs), constants (pi, e), and nested parentheses
- Support improper integrals with infinite bounds via --infinite flag using appropriate limit-based techniques
- Support tabulated data integration: read (x, y) pairs from a CSV file via --data flag and integrate using the available methods
- Print results to console: integral value, estimated error, method used, intervals used, and computation time
- Save results as JSON with --output flag (default: integration_result.json)
- If no input is given, compute several well-known integrals (integral of sin(x) from 0 to pi = 2, integral of exp(-x^2) from 0 to infinity = sqrt(pi)/2, integral of 1/x from 1 to e = 1) using all methods and display a comparison table with exact vs computed values and relative errors
- Handle errors: division by zero within the integration range, non-convergent integrals, invalid mathematical expressions, and bounds where lower > upper

Create a complete Python project for a clean Ubuntu 22.04 machine with only Python 3.10+ installed. Include:
- Source code
- requirements.txt with all dependencies (direct and transitive) pinned to exact versions
- README.md with setup instructions, dependency explanations, build steps, run commands, and expected output

[]

"""
Numerical Integration Calculator
Computes definite integrals using Trapezoidal, Simpson's, Gauss-Legendre,
and Romberg methods with convergence comparison.

Dependencies: scipy 1.12.0, numpy 1.26.4, matplotlib 3.8.2
"""

import numpy as np
from scipy import integrate
from scipy.special import roots_legendre
import matplotlib.pyplot as plt
import sys
import time


def trapezoidal(f, a, b, n):
    """Compute definite integral using the composite trapezoidal rule."""
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])


def simpsons(f, a, b, n):
    """Compute definite integral using Simpson's 1/3 rule.
    n must be even; if odd, it is incremented by 1.
    """
    if n % 2 != 0:
        n += 1
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return (h / 3) * (y[0] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]) + y[-1])


def gauss_legendre(f, a, b, n):
    """Compute definite integral using Gauss-Legendre quadrature with n nodes."""
    nodes, weights = roots_legendre(n)
    # Transform from [-1, 1] to [a, b]
    transformed_nodes = 0.5 * (b - a) * nodes + 0.5 * (a + b)
    transformed_weights = 0.5 * (b - a) * weights
    return np.sum(transformed_weights * f(transformed_nodes))


def romberg(f, a, b, max_order=10, tol=1e-12):
    """Compute definite integral using Romberg integration.
    Returns the result and the Romberg table for convergence analysis.
    """
    R = np.zeros((max_order, max_order))
    h = b - a
    R[0, 0] = 0.5 * h * (f(a) + f(b))

    for i in range(1, max_order):
        h_i = h / (2 ** i)
        # Composite trapezoidal with 2^i subintervals
        midpoints = a + h_i * (2 * np.arange(1, 2 ** (i - 1) + 1) - 1)
        R[i, 0] = 0.5 * R[i - 1, 0] + h_i * np.sum(f(midpoints))

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

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

    return R[max_order - 1, max_order - 1], R


def convergence_study(f, a, b, exact_value, max_n=256):
    """Perform a convergence study comparing all four methods."""
    ns = [2 ** k for k in range(1, int(np.log2(max_n)) + 1)]

    errors_trap = []
    errors_simp = []
    errors_gauss = []
    errors_romb = []

    for n in ns:
        # Trapezoidal
        result_trap = trapezoidal(f, a, b, n)
        errors_trap.append(abs(result_trap - exact_value))

        # Simpson's (needs even n)
        n_simp = n if n % 2 == 0 else n + 1
        result_simp = simpsons(f, a, b, n_simp)
        errors_simp.append(abs(result_simp - exact_value))

        # Gauss-Legendre (use n nodes, capped at reasonable value)
        n_gauss = min(n, 64)
        result_gauss = gauss_legendre(f, a, b, n_gauss)
        errors_gauss.append(abs(result_gauss - exact_value))

        # Romberg (use log2(n) order)
        order = max(2, int(np.log2(n)))
        result_romb, _ = romberg(f, a, b, max_order=order)
        errors_romb.append(abs(result_romb - exact_value))

    return ns, errors_trap, errors_simp, errors_gauss, errors_romb


def plot_convergence(ns, errors_trap, errors_simp, errors_gauss, errors_romb,
                     title="Convergence Comparison", filename="convergence.png"):
    """Plot convergence comparison of all four methods."""
    plt.figure(figsize=(10, 7))

    # Replace zeros with minimum float for log plot
    min_err = 1e-16
    errors_trap = [max(e, min_err) for e in errors_trap]
    errors_simp = [max(e, min_err) for e in errors_simp]
    errors_gauss = [max(e, min_err) for e in errors_gauss]
    errors_romb = [max(e, min_err) for e in errors_romb]

    plt.loglog(ns, errors_trap, 'o-', label='Trapezoidal', linewidth=2, markersize=6)
    plt.loglog(ns, errors_simp, 's-', label="Simpson's 1/3", linewidth=2, markersize=6)
    plt.loglog(ns, errors_gauss, '^-', label='Gauss-Legendre', linewidth=2, markersize=6)
    plt.loglog(ns, errors_romb, 'D-', label='Romberg', linewidth=2, markersize=6)

    # Reference slopes
    n_arr = np.array(ns, dtype=float)
    plt.loglog(ns, (n_arr[0] / n_arr) ** 2 * errors_trap[0], '--', color='gray',
               alpha=0.5, label='O(n^-2) reference')
    plt.loglog(ns, (n_arr[0] / n_arr) ** 4 * errors_simp[0], ':', color='gray',
               alpha=0.5, label='O(n^-4) reference')

    plt.xlabel('Number of points (n)', fontsize=12)
    plt.ylabel('Absolute Error', fontsize=12)
    plt.title(title, fontsize=14)
    plt.legend(fontsize=10)
    plt.grid(True, which='both', alpha=0.3)
    plt.tight_layout()
    plt.savefig(filename, dpi=150)
    plt.close()
    print(f"Convergence plot saved to {filename}")


def scipy_reference(f, a, b):
    """Compute reference value using scipy.integrate.quad."""
    result, error = integrate.quad(f, a, b)
    return result, error


# --- Test functions with known exact integrals ---

TEST_FUNCTIONS = {
    "polynomial": {
        "f": lambda x: 3 * x ** 4 - 2 * x ** 3 + x ** 2 - 5 * x + 7,
        "a": 0, "b": 2,
        "exact": 3 * (2 ** 5) / 5 - 2 * (2 ** 4) / 4 + (2 ** 3) / 3 - 5 * (2 ** 2) / 2 + 7 * 2,
        "desc": "3x^4 - 2x^3 + x^2 - 5x + 7 on [0, 2]"
    },
    "trigonometric": {
        "f": lambda x: np.sin(x) * np.cos(x),
        "a": 0, "b": np.pi / 2,
        "exact": 0.5,  # sin^2(pi/2)/2 - sin^2(0)/2
        "desc": "sin(x)*cos(x) on [0, pi/2]"
    },
    "exponential": {
        "f": lambda x: np.exp(-x ** 2),
        "a": 0, "b": 1,
        "exact": 0.7468241328124271,  # erf(1)*sqrt(pi)/2
        "desc": "exp(-x^2) on [0, 1] (Gaussian)"
    },
    "oscillatory": {
        "f": lambda x: np.sin(10 * x) * np.exp(-x),
        "a": 0, "b": np.pi,
        "exact": None,  # Will compute via scipy
        "desc": "sin(10x)*exp(-x) on [0, pi] (oscillatory)"
    },
    "singular_endpoint": {
        "f": lambda x: np.where(x > 0, 1.0 / np.sqrt(x + 1e-15), 0.0),
        "a": 1e-10, "b": 1,
        "exact": 2.0 * (1.0 - np.sqrt(1e-10)),
        "desc": "1/sqrt(x) on [~0, 1] (near-singular)"
    }
}


def run_single_test(name, test_info, verbose=True):
    """Run all integration methods on a single test function."""
    f = test_info["f"]
    a, b = test_info["a"], test_info["b"]

    # Get exact value (use scipy if not provided)
    if test_info["exact"] is not None:
        exact = test_info["exact"]
    else:
        exact, _ = scipy_reference(f, a, b)

    if verbose:
        print(f"\n{'=' * 70}")
        print(f"Test: {test_info['desc']}")
        print(f"Exact value: {exact:.15f}")
        print(f"{'=' * 70}")

    results = {}
    n_points = 100

    # Trapezoidal
    t0 = time.perf_counter()
    val = trapezoidal(f, a, b, n_points)
    elapsed = time.perf_counter() - t0
    results["Trapezoidal"] = {"value": val, "error": abs(val - exact), "time": elapsed}

    # Simpson's
    t0 = time.perf_counter()
    val = simpsons(f, a, b, n_points)
    elapsed = time.perf_counter() - t0
    results["Simpson's"] = {"value": val, "error": abs(val - exact), "time": elapsed}

    # Gauss-Legendre
    t0 = time.perf_counter()
    val = gauss_legendre(f, a, b, 20)
    elapsed = time.perf_counter() - t0
    results["Gauss-Legendre"] = {"value": val, "error": abs(val - exact), "time": elapsed}

    # Romberg
    t0 = time.perf_counter()
    val, table = romberg(f, a, b, max_order=10)
    elapsed = time.perf_counter() - t0
    results["Romberg"] = {"value": val, "error": abs(val - exact), "time": elapsed}

    # scipy.integrate.quad reference
    t0 = time.perf_counter()
    val, err = scipy_reference(f, a, b)
    elapsed = time.perf_counter() - t0
    results["scipy.quad"] = {"value": val, "error": abs(val - exact), "time": elapsed}

    if verbose:
        print(f"\n{'Method':<20} {'Result':<22} {'Error':<15} {'Time (ms)':<12}")
        print("-" * 70)
        for method, data in results.items():
            print(f"{method:<20} {data['value']:<22.15f} {data['error']:<15.2e} "
                  f"{data['time'] * 1000:<12.4f}")

    return results, exact


def main():
    """Main function - run all tests and generate convergence plots."""
    print("=" * 70)
    print("  Numerical Integration Calculator")
    print("  Methods: Trapezoidal, Simpson's, Gauss-Legendre, Romberg")
    print("=" * 70)

    all_results = {}

    for name, test_info in TEST_FUNCTIONS.items():
        results, exact = run_single_test(name, test_info)
        all_results[name] = {"results": results, "exact": exact}

    # Convergence study on a smooth function
    print(f"\n\n{'=' * 70}")
    print("  Convergence Study: exp(-x^2) on [0, 1]")
    print(f"{'=' * 70}")

    f = TEST_FUNCTIONS["exponential"]["f"]
    exact = TEST_FUNCTIONS["exponential"]["exact"]
    ns, e_trap, e_simp, e_gauss, e_romb = convergence_study(f, 0, 1, exact, max_n=512)

    print(f"\n{'n':<8} {'Trapezoidal':<15} {'Simpsons':<15} {'Gauss-Legendre':<15} {'Romberg':<15}")
    print("-" * 68)
    for i, n in enumerate(ns):
        print(f"{n:<8} {e_trap[i]:<15.2e} {e_simp[i]:<15.2e} "
              f"{e_gauss[i]:<15.2e} {e_romb[i]:<15.2e}")

    plot_convergence(ns, e_trap, e_simp, e_gauss, e_romb,
                     title="Convergence: exp(-x^2) on [0, 1]",
                     filename="convergence_gaussian.png")

    # Convergence study on oscillatory function
    f_osc = TEST_FUNCTIONS["oscillatory"]["f"]
    exact_osc = all_results["oscillatory"]["exact"]
    ns2, e_trap2, e_simp2, e_gauss2, e_romb2 = convergence_study(
        f_osc, 0, np.pi, exact_osc, max_n=512)

    plot_convergence(ns2, e_trap2, e_simp2, e_gauss2, e_romb2,
                     title="Convergence: sin(10x)*exp(-x) on [0, pi]",
                     filename="convergence_oscillatory.png")

    # Summary
    print(f"\n\n{'=' * 70}")
    print("  Summary: Best method for each test function")
    print(f"{'=' * 70}")
    for name, data in all_results.items():
        best_method = min(
            [(m, d["error"]) for m, d in data["results"].items()],
            key=lambda x: x[1]
        )
        print(f"  {name:<25} -> {best_method[0]:<20} (error: {best_method[1]:.2e})")

    print("\nDone. Convergence plots saved.")
    return 0


if __name__ == "__main__":
    sys.exit(main())

# Numerical Integration Calculator (Python - Trial 1)

Computes definite integrals using four numerical methods with convergence comparison.

## Methods
- **Trapezoidal Rule**: Composite trapezoidal integration
- **Simpson's 1/3 Rule**: Composite Simpson's rule (requires even subintervals)
- **Gauss-Legendre Quadrature**: High-accuracy quadrature using optimal nodes/weights
- **Romberg Integration**: Richardson extrapolation on trapezoidal estimates

## Dependencies
- scipy 1.12.0 - Gauss-Legendre roots/weights, reference quad integration
- numpy 1.26.4 - Array operations and mathematical functions
- matplotlib 3.8.2 - Convergence comparison plots

## Setup
```bash
pip install -r requirements.txt
```

## Usage
```bash
python integrator.py
```

## Output
- Console output with integration results and errors for each test function
- Convergence comparison tables showing error vs number of points
- PNG plots comparing convergence rates of all four methods

scipy==1.12.0
numpy==1.26.4
matplotlib==3.8.2