Numerical Integration Calculator (javascript, written by Claude Code)

Claude Code wrote this javascript 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 JavaScript project for a clean Ubuntu 22.04 machine with only Node.js 20+ (LTS) installed. Include:
- Source code
- package.json 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: mathjs 12.4.1, chalk 4.1.2
 */

const math = require('mathjs');
const chalk = require('chalk');

// ---------------------------------------------------------------------------
// Trapezoidal Rule
// ---------------------------------------------------------------------------
function trapezoidal(f, a, b, n) {
    const h = (b - a) / n;
    let sum = 0.5 * (f(a) + f(b));
    for (let i = 1; i < n; i++) {
        sum += f(a + i * h);
    }
    return sum * h;
}

// ---------------------------------------------------------------------------
// Simpson's 1/3 Rule (n must be even)
// ---------------------------------------------------------------------------
function simpsons(f, a, b, n) {
    if (n % 2 !== 0) n++;
    const h = (b - a) / n;
    let sum = f(a) + f(b);
    for (let i = 1; i < n; i += 2) {
        sum += 4 * f(a + i * h);
    }
    for (let i = 2; i < n; i += 2) {
        sum += 2 * f(a + i * h);
    }
    return (sum * h) / 3;
}

// ---------------------------------------------------------------------------
// Gauss-Legendre Quadrature (5-point)
// ---------------------------------------------------------------------------
function gaussLegendre(f, a, b, n) {
    // Standard Gauss-Legendre nodes and weights for 5-point quadrature
    const nodes5 = [
        -0.9061798459386640, -0.5384693101056831, 0.0,
         0.5384693101056831,  0.9061798459386640
    ];
    const weights5 = [
        0.2369268850561891, 0.4786286704993665, 0.5688888888888889,
        0.4786286704993665, 0.2369268850561891
    ];

    // Split interval into sub-intervals for higher n
    const numSub = Math.max(1, Math.floor(n / 5));
    const subLen = (b - a) / numSub;
    let total = 0;

    for (let s = 0; s < numSub; s++) {
        const subA = a + s * subLen;
        const subB = subA + subLen;
        const mid = 0.5 * (subA + subB);
        const half = 0.5 * (subB - subA);

        for (let i = 0; i < 5; i++) {
            const x = mid + half * nodes5[i];
            total += weights5[i] * f(x) * half;
        }
    }
    return total;
}

// ---------------------------------------------------------------------------
// Romberg Integration
// ---------------------------------------------------------------------------
function romberg(f, a, b, maxOrder = 10) {
    const R = Array.from({ length: maxOrder }, () => new Array(maxOrder).fill(0));
    const h = b - a;
    R[0][0] = 0.5 * h * (f(a) + f(b));

    for (let i = 1; i < maxOrder; i++) {
        const hi = h / Math.pow(2, i);
        const numNewPoints = Math.pow(2, i - 1);
        let sum = 0;
        for (let k = 1; k <= numNewPoints; k++) {
            sum += f(a + hi * (2 * k - 1));
        }
        R[i][0] = 0.5 * R[i - 1][0] + hi * sum;

        for (let j = 1; j <= i; j++) {
            const factor = Math.pow(4, j);
            R[i][j] = (factor * R[i][j - 1] - R[i - 1][j - 1]) / (factor - 1);
        }

        if (i > 0 && Math.abs(R[i][i] - R[i - 1][i - 1]) < 1e-12) {
            return { value: R[i][i], order: i + 1, table: R.slice(0, i + 1) };
        }
    }
    return { value: R[maxOrder - 1][maxOrder - 1], order: maxOrder, table: R };
}

// ---------------------------------------------------------------------------
// mathjs-based reference integrator (using numeric integration)
// ---------------------------------------------------------------------------
function mathjsReference(exprStr, a, b) {
    try {
        const node = math.parse(exprStr);
        const compiled = node.compile();
        const f = (x) => compiled.evaluate({ x });

        // Use Simpson's with high n as reference
        return simpsons(f, a, b, 10000);
    } catch (e) {
        return NaN;
    }
}

// ---------------------------------------------------------------------------
// Test functions
// ---------------------------------------------------------------------------
const TEST_FUNCTIONS = {
    polynomial: {
        f: x => 3 * Math.pow(x, 4) - 2 * Math.pow(x, 3) + x * x - 5 * x + 7,
        a: 0, b: 2,
        exact: 3 * Math.pow(2, 5) / 5 - 2 * Math.pow(2, 4) / 4 + Math.pow(2, 3) / 3 - 5 * 4 / 2 + 14,
        desc: '3x^4 - 2x^3 + x^2 - 5x + 7 on [0, 2]'
    },
    trigonometric: {
        f: x => Math.sin(x) * Math.cos(x),
        a: 0, b: Math.PI / 2,
        exact: 0.5,
        desc: 'sin(x)*cos(x) on [0, pi/2]'
    },
    exponential: {
        f: x => Math.exp(-x * x),
        a: 0, b: 1,
        exact: 0.7468241328124271,
        desc: 'exp(-x^2) on [0, 1] (Gaussian)'
    },
    oscillatory: {
        f: x => Math.sin(10 * x) * Math.exp(-x),
        a: 0, b: Math.PI,
        exact: null,
        desc: 'sin(10x)*exp(-x) on [0, pi] (oscillatory)'
    },
    singular_endpoint: {
        f: x => x > 0 ? 1.0 / Math.sqrt(x) : 0,
        a: 1e-10, b: 1,
        exact: 2.0 * (1.0 - Math.sqrt(1e-10)),
        desc: '1/sqrt(x) on [~0, 1] (near-singular)'
    }
};

// ---------------------------------------------------------------------------
// Convergence study
// ---------------------------------------------------------------------------
function convergenceStudy(f, a, b, exact, maxN = 256) {
    const ns = [];
    for (let k = 1; (1 << k) <= maxN; k++) ns.push(1 << k);

    const errTrap = [], errSimp = [], errGauss = [], errRomb = [];

    for (const n of ns) {
        errTrap.push(Math.abs(trapezoidal(f, a, b, n) - exact));
        const nSimp = n % 2 === 0 ? n : n + 1;
        errSimp.push(Math.abs(simpsons(f, a, b, nSimp) - exact));
        const nGauss = Math.min(n, 64);
        errGauss.push(Math.abs(gaussLegendre(f, a, b, nGauss) - exact));
        const order = Math.max(2, Math.floor(Math.log2(n)));
        errRomb.push(Math.abs(romberg(f, a, b, order).value - exact));
    }

    return { ns, errTrap, errSimp, errGauss, errRomb };
}

// ---------------------------------------------------------------------------
// Run a single test
// ---------------------------------------------------------------------------
function runSingleTest(name, testInfo) {
    let exact = testInfo.exact;
    if (exact === null) {
        exact = simpsons(testInfo.f, testInfo.a, testInfo.b, 10000);
    }

    console.log('\n' + chalk.cyan('='.repeat(70)));
    console.log(chalk.cyan(`Test: ${testInfo.desc}`));
    console.log(chalk.cyan(`Exact value: ${exact.toFixed(15)}`));
    console.log(chalk.cyan('='.repeat(70)));

    const results = {};
    const nPoints = 100;

    // Trapezoidal
    let t0 = process.hrtime.bigint();
    let val = trapezoidal(testInfo.f, testInfo.a, testInfo.b, nPoints);
    let elapsed = Number(process.hrtime.bigint() - t0) / 1e6;
    results['Trapezoidal'] = { value: val, error: Math.abs(val - exact), time_ms: elapsed };

    // Simpson's
    t0 = process.hrtime.bigint();
    val = simpsons(testInfo.f, testInfo.a, testInfo.b, nPoints);
    elapsed = Number(process.hrtime.bigint() - t0) / 1e6;
    results["Simpson's"] = { value: val, error: Math.abs(val - exact), time_ms: elapsed };

    // Gauss-Legendre
    t0 = process.hrtime.bigint();
    val = gaussLegendre(testInfo.f, testInfo.a, testInfo.b, 20);
    elapsed = Number(process.hrtime.bigint() - t0) / 1e6;
    results['Gauss-Legendre'] = { value: val, error: Math.abs(val - exact), time_ms: elapsed };

    // Romberg
    t0 = process.hrtime.bigint();
    const rResult = romberg(testInfo.f, testInfo.a, testInfo.b, 10);
    val = rResult.value;
    elapsed = Number(process.hrtime.bigint() - t0) / 1e6;
    results['Romberg'] = { value: val, error: Math.abs(val - exact), time_ms: elapsed };

    // Print table
    console.log(`\n${'Method'.padEnd(20)} ${'Result'.padEnd(22)} ${'Error'.padEnd(15)} ${'Time (ms)'.padEnd(12)}`);
    console.log('-'.repeat(70));
    for (const [method, data] of Object.entries(results)) {
        const errStr = data.error.toExponential(2);
        console.log(
            `${method.padEnd(20)} ${data.value.toFixed(15).padEnd(22)} ${errStr.padEnd(15)} ${data.time_ms.toFixed(4).padEnd(12)}`
        );
    }

    return { results, exact };
}

// ---------------------------------------------------------------------------
// Main
// ---------------------------------------------------------------------------
function main() {
    console.log(chalk.bold('='.repeat(70)));
    console.log(chalk.bold('  Numerical Integration Calculator'));
    console.log(chalk.bold('  Methods: Trapezoidal, Simpson\'s, Gauss-Legendre, Romberg'));
    console.log(chalk.bold('='.repeat(70)));

    const allResults = {};
    for (const [name, testInfo] of Object.entries(TEST_FUNCTIONS)) {
        const { results, exact } = runSingleTest(name, testInfo);
        allResults[name] = { results, exact };
    }

    // Convergence study
    console.log('\n\n' + chalk.yellow('='.repeat(70)));
    console.log(chalk.yellow('  Convergence Study: exp(-x^2) on [0, 1]'));
    console.log(chalk.yellow('='.repeat(70)));

    const fGauss = TEST_FUNCTIONS.exponential.f;
    const exactGauss = TEST_FUNCTIONS.exponential.exact;
    const conv = convergenceStudy(fGauss, 0, 1, exactGauss, 512);

    console.log(`\n${'n'.padEnd(8)} ${'Trapezoidal'.padEnd(15)} ${"Simpson's".padEnd(15)} ${'Gauss-Legendre'.padEnd(15)} ${'Romberg'.padEnd(15)}`);
    console.log('-'.repeat(68));
    for (let i = 0; i < conv.ns.length; i++) {
        console.log(
            `${String(conv.ns[i]).padEnd(8)} ${conv.errTrap[i].toExponential(2).padEnd(15)} ` +
            `${conv.errSimp[i].toExponential(2).padEnd(15)} ${conv.errGauss[i].toExponential(2).padEnd(15)} ` +
            `${conv.errRomb[i].toExponential(2).padEnd(15)}`
        );
    }

    // Summary
    console.log('\n' + chalk.green('='.repeat(70)));
    console.log(chalk.green('  Summary: Best method for each test function'));
    console.log(chalk.green('='.repeat(70)));
    for (const [name, data] of Object.entries(allResults)) {
        let bestMethod = '';
        let bestError = Infinity;
        for (const [method, d] of Object.entries(data.results)) {
            if (d.error < bestError) {
                bestError = d.error;
                bestMethod = method;
            }
        }
        console.log(`  ${name.padEnd(25)} -> ${bestMethod.padEnd(20)} (error: ${bestError.toExponential(2)})`);
    }

    console.log('\nDone.');
}

main();

{
  "name": "numerical-integration-trial1",
  "version": "1.0.0",
  "description": "Numerical Integration Calculator using mathjs and chalk",
  "main": "integrator.js",
  "scripts": {
    "start": "node integrator.js"
  },
  "dependencies": {
    "mathjs": "12.4.1",
    "chalk": "4.1.2"
  }
}