Python numerical integration

Question

Could the time complexity of this definite integral algorithm be improved?

import numpy as np

def integral(f, left, right, epsilon):
    w = int((right-left)/epsilon)
    domain = np.linspace(left, right, w)

    heights = [f(((a+b)*0.5)) for a, b in zip(domain[:-1], domain[1:])]

    return sum(heights)*epsilon

# Test an example
print('Result:', integral(lambda x: np.sin(x), 1, 10, 0.01))

Answer 1

By changing the start and end points to start+epsilon/2 and end-epsilon/2, we can simplify the calculation tremendously. This means that the domain is what we want to evaluate the function at, removing the zip and averaging logic, secondly, because we now have heights of the same length as domain, we can just call f(domain), to have it update in place. Finally, we can use np.sum to return the answer. Put together, these changes result in simpler code, and roughly 100x faster (especially for small epsilon)

def integral1(f, left, right, epsilon):
    left += epsilon / 2
    right -= epsilon / 2
    domain = np.linspace(left, right, (right-left) // epsilon)
    return np.sum(f(domain), dtype = np.float128) * epsilon

Answer 2

You might want to check out quadpy (a project of mine). It does the heavy lifting for all sorts of integration problems.

import numpy
import quadpy

val, error_estimate = quadpy.line_segment.integrate_adaptive(
       numpy.sin,
       [0.0, 1.0],
       1.0e-10
       )