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))

• Mast: Is there a reason you went for a lambda?
• Pii: No, but that's just an input that the user gives me. I really only care about the efficiency of integral

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). The memory usage will also be reduced.

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

Just for the record, I'm assuming that you aren't looking for general advice with numerical integration, and rather just using the midpoint rule. If you were, Simpson's rule or similar would be much more accurate, and would not take longer to compute.

It turns out that accumulating with a 128 bit float adds significant accuracy in some pathological cases, while only taking 10% longer.

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
 )