Compute indefinite integrals by the trapezoidal rule

trapz.py

import numpy as np


def running_trapz(f, mesh):
    """Return the 1D indefinite integral of the function *f* given
    at the points in *mesh*. The integral is returned at the points
    of *mesh*, and the leftmost value is always zero. The integral
    is approximated using the trapezoidal rule.
    """

    dx = np.diff(mesh)
    weights = np.zeros_like(mesh)
    weights[0:-1] += .5*dx
    weights[1:] += .5*dx

    result = np.empty_like(f)
    result[0] = 0
    result[1:] = 0.5*np.cumsum(dx*f[:-1] + dx*f[1:])
    return result


def test_running_trapz(n):
    x = np.linspace(0, 5, n, endpoint=True)
    h = x[1] - x[0]

    soln = running_trapz(np.sin(x), x)

    err = soln - (1-np.cos(x))

    if 0:
        import matplotlib.pyplot as pt
        pt.plot(x, soln, label="num")
        pt.plot(x, 1-np.cos(x), label="true")
        pt.legend()
        pt.show()

    return h, np.sqrt(np.trapz(err**2, x))


def estimate_order(f, point_counts):
    n1, n2 = point_counts
    h1, err1 = f(n1)
    h2, err2 = f(n2)

    print "h=%g err=%g" % (h1, err1)
    print "h=%g err=%g" % (h2, err2)

    from math import log
    est_order = (log(err2/err1) / log(h2/h1))
    print "EOC: %g" % est_order

    return est_order


if __name__ == "__main__":
    eoc = estimate_order(test_running_trapz, [20, 40])
    assert eoc > 1.85