Python numerical integration
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