Numeric double integration I've made a simple program for numerically aproximating double integral, which accepts that the bounds of the inner integral are functions:

import numpy as np
import time

def double_integral(func, limits, res=1000):
        t = time.clock()
        t1 = time.clock()
        t2 = time.clock()
        s = 0
        a, b = limits[0], limits[1]
        outer_values = np.linspace(a, b, res)
        c_is_func = callable(limits[2])
        d_is_func = callable(limits[3])
        for y in outer_values:
            if c_is_func:
                c = limits[2](y)
            else:
                c = limits[2]
            if d_is_func:
                d = limits[3](y)
            else:
                d = limits[3]
            dA = ((b - a) / res) * ((d - c) / res)
            inner_values = np.linspace(c, d, res)
            for x in inner_values:
                t2 = time.clock() - t2
                s += func(x, y) * dA
            t1 = time.clock() - t1
         t = time.clock() - t
         return s, t, t1 / res, t2 / res**2


This is, however, terribly slow. When res=1000, such that the integral is a sum of a million parts, it takes about 5 seconds to run, but the answer is only correct to about the 3rd decimal in my experience. Is there any way to speed this up?

The code i am running to check the integral is

def f(x, y):
    if (4 - y**2 - x**2) &lt; 0:
        return 0              #This is to avoid taking the root of negarive #'s
    return np.sqrt(4 - y**2 - x**2)

def c(y):
    return np.sqrt(2 * y - y**2)

def d(y):
    return np.sqrt(4 - y**2)
# b d
# S S f(x,y) dx dy
# a c
a, b, = 0, 2
print(double_integral(f, [a, b, c, d]))


The integral is eaqual to 16/9
If you want to use numpy, use numpy properly. Inestead of 

for x in inner_values:
    s += func(x, y) * dA


use the more idiomatic, and much faster

s += dA * np.sum(func(inner_values, y))


Note: this requires rewriting f as 

return np.sqrt(np.maximum(0, 4 - y**2 - x**2))


so it can accept an array as an input. This does not reduce accuracy, but brings time down to a much more accaptable .04 seconds for a 100x improvement. The takeaway here is Numpy is not magic. It provides quick vectorization.