Numerical Integration Calculator (python, written by Codex)

Codex 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 (Python)

## Requirements
- Ubuntu 22.04
- Python 3.10+

## Install
```bash
python3 -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt
```

## Run
```bash
python src/main.py --function "sin(x)*exp(-x)" --lower 0 --upper 10 --method simpson
python src/main.py --function "exp(-x^2)" --lower 0 --upper inf --method romberg --infinite
python src/main.py --function "1/x" --lower 1 --upper 2 --compare
python src/main.py --data points.csv --method trapezoidal
```

If no input is provided, well-known integrals are computed and compared.

numpy==2.1.1
scipy==1.14.1
sympy==1.13.3

#!/usr/bin/env python3
import argparse
import csv
import json
import math
import sys
import time
from pathlib import Path

import numpy as np
from scipy import integrate
import sympy as sp


x = sp.symbols("x")


def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(description="Numerical Integration Calculator")
    p.add_argument("--function", dest="fn")
    p.add_argument("--lower")
    p.add_argument("--upper")
    p.add_argument("--method", choices=["trapezoidal", "simpson", "simpson38", "gaussian", "romberg"], default="trapezoidal")
    p.add_argument("--intervals", type=int, default=1000)
    p.add_argument("--compare", action="store_true")
    p.add_argument("--tolerance", type=float, default=1e-10)
    p.add_argument("--infinite", action="store_true")
    p.add_argument("--data")
    p.add_argument("--output", default="integration_result.json")
    return p.parse_args()


def parse_bound(v):
    if v is None:
        return None
    s = str(v).strip().lower()
    if s in {"inf", "+inf", "infinity", "+infinity"}:
        return math.inf
    if s in {"-inf", "-infinity"}:
        return -math.inf
    return float(v)


def compile_fn(expr: str):
    e = sp.sympify(expr)
    f = sp.lambdify(x, e, modules=["numpy", "math"])
    return lambda xv: float(f(xv))


def trap(f, a, b, n):
    xs = np.linspace(a, b, n + 1)
    ys = np.array([f(v) for v in xs])
    return float(np.trapz(ys, xs))


def simp(f, a, b, n):
    if n % 2:
        n += 1
    xs = np.linspace(a, b, n + 1)
    ys = np.array([f(v) for v in xs])
    return float(integrate.simpson(ys, x=xs))


def simp38(f, a, b, n):
    while n % 3:
        n += 1
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        s += (2 if i % 3 == 0 else 3) * f(a + i * h)
    return float(3 * h * s / 8)


def gauss5(f, a, b):
    xs = np.array([0, -0.5384693101056831, 0.5384693101056831, -0.9061798459386640, 0.9061798459386640])
    ws = np.array([0.5688888888888889, 0.4786286704993665, 0.4786286704993665, 0.2369268850561891, 0.2369268850561891])
    c1, c2 = (b - a) / 2, (b + a) / 2
    return float(c1 * np.sum(ws * np.array([f(c1 * t + c2) for t in xs])))


def romberg(f, a, b):
    return float(integrate.romberg(f, a, b, divmax=8, tol=1e-12, rtol=1e-12, show=False))


def integrate_method(name, f, a, b, n):
    if name == "trapezoidal":
        return trap(f, a, b, n)
    if name == "simpson":
        return simp(f, a, b, n)
    if name == "simpson38":
        return simp38(f, a, b, n)
    if name == "gaussian":
        return gauss5(f, a, b)
    if name == "romberg":
        return romberg(f, a, b)
    raise ValueError(name)


def adaptive(name, f, a, b, n, tol):
    prev = integrate_method(name, f, a, b, n)
    cur_n = n
    for _ in range(20):
        cur_n *= 2
        cur = integrate_method(name, f, a, b, cur_n)
        if abs(cur - prev) < tol:
            return {"value": cur, "intervals": cur_n, "estimated_error": abs(cur - prev)}
        prev = cur
    return {"value": prev, "intervals": cur_n, "estimated_error": float("nan")}


def handle_infinite(f, a, b, method, n, tol):
    if math.isfinite(a) and math.isfinite(b):
        return adaptive(method, f, a, b, n, tol)
    if math.isfinite(a) and not math.isfinite(b):
        g = lambda t: f(a + t / (1 - t)) / ((1 - t) ** 2)
        return adaptive(method, g, 0, 1 - 1e-8, n, tol)
    if not math.isfinite(a) and math.isfinite(b):
        g = lambda t: f(b - t / (1 - t)) / ((1 - t) ** 2)
        return adaptive(method, g, 0, 1 - 1e-8, n, tol)
    g = lambda t: f(math.tan(math.pi * (t - 0.5))) * math.pi / (math.cos(math.pi * (t - 0.5)) ** 2)
    return adaptive(method, g, 1e-8, 1 - 1e-8, n, tol)


def integrate_data(path: Path, method: str):
    rows = []
    with path.open("r", encoding="utf-8") as f:
        r = csv.reader(f)
        for row in r:
            if not row:
                continue
            rows.append((float(row[0]), float(row[1])))
    rows.sort(key=lambda t: t[0])
    x = np.array([r[0] for r in rows])
    y = np.array([r[1] for r in rows])

    if method == "trapezoidal":
        return float(np.trapz(y, x))

    # linear interpolation + simpson fallback
    f = lambda xp: float(np.interp(xp, x, y))
    return simp(f, float(x[0]), float(x[-1]), 1000)


def compare_all(f, a, b, n, tol):
    methods = ["trapezoidal", "simpson", "simpson38", "gaussian", "romberg"]
    out = []
    for m in methods:
        t0 = time.perf_counter()
        r = handle_infinite(f, a, b, m, n, tol)
        out.append({"method": m, "value": r["value"], "estimated_error": r["estimated_error"], "intervals": r["intervals"], "time_ms": (time.perf_counter() - t0) * 1000})
    return out


def known_demo():
    cases = [
        ("integral sin(x) [0,pi]", "sin(x)", 0, math.pi, 2.0),
        ("integral exp(-x^2) [0,inf]", "exp(-x^2)", 0, math.inf, math.sqrt(math.pi) / 2),
        ("integral 1/x [1,e]", "1/x", 1, math.e, 1.0),
    ]
    out = []
    for name, expr, a, b, exact in cases:
        f = compile_fn(expr)
        for m in ["trapezoidal", "simpson", "simpson38", "gaussian", "romberg"]:
            r = handle_infinite(f, a, b, m, 1000, 1e-10)
            out.append({"case": name, "method": m, "computed": r["value"], "exact": exact, "relative_error": abs((r["value"] - exact) / exact)})
    return out


def main() -> int:
    args = parse_args()

    try:
        if not args.fn and not args.data:
            comp = known_demo()
            print(json.dumps(comp, indent=2))
            Path(args.output).write_text(json.dumps({"config": vars(args), "comparison": comp}, indent=2), encoding="utf-8")
            return 0

        if args.data:
            t0 = time.perf_counter()
            val = integrate_data(Path(args.data), args.method)
            result = {"mode": "tabulated", "method": args.method, "value": val, "time_ms": (time.perf_counter() - t0) * 1000}
            print(json.dumps(result, indent=2))
            Path(args.output).write_text(json.dumps({"config": vars(args), "result": result}, indent=2), encoding="utf-8")
            return 0

        a = parse_bound(args.lower)
        b = parse_bound(args.upper)
        if a is None or b is None:
            print("Error: --lower and --upper are required", file=sys.stderr)
            return 1
        if math.isfinite(a) and math.isfinite(b) and a > b:
            print("Error: lower bound cannot exceed upper bound", file=sys.stderr)
            return 1

        f = compile_fn(args.fn)

        if args.compare:
            comp = compare_all(f, a, b, args.intervals, args.tolerance)
            print(json.dumps(comp, indent=2))
            Path(args.output).write_text(json.dumps({"config": vars(args), "comparison": comp}, indent=2), encoding="utf-8")
            return 0

        t0 = time.perf_counter()
        r = handle_infinite(f, a, b, args.method, args.intervals, args.tolerance)
        result = {"method": args.method, "value": r["value"], "estimated_error": r["estimated_error"], "intervals": r["intervals"], "time_ms": (time.perf_counter() - t0) * 1000}
        print(json.dumps(result, indent=2))
        Path(args.output).write_text(json.dumps({"config": vars(args), "result": result}, indent=2), encoding="utf-8")
        return 0
    except Exception as exc:
        print(f"Error: {exc}", file=sys.stderr)
        return 1


if __name__ == "__main__":
    raise SystemExit(main())