# 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, 2017 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, 2017 at 18:50

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 Sep 29, 2017 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. Oct 4, 2017 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. Oct 5, 2017 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