# Python numerical integration

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

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

## 2 Answers

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

• oh, and also the memory usage will be smaller – Oscar Smith Sep 29 '17 at 4:10
• 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. – Oscar Smith Oct 4 '17 at 22:21
• It turns out that accumulating with a 128 bit float adds significant accuracy in some pathalogical cases, while only taking 10% longer. Updated answer. – Oscar Smith Oct 5 '17 at 18:23

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
)