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)
                c = limits[2]
            if d_is_func:
                d = limits[3](y)
                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) < 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

  • 3
    \$\begingroup\$ (Welcome to CR!) Don't even think of the clock in inner loops. Have a look at timeit. Do you know SciPy integrate? \$\endgroup\$
    – greybeard
    Commented Feb 22, 2018 at 20:53
  • 2
    \$\begingroup\$ You really ought to include all your import statements; I'm assuming at the very least you have import numpy as np right? \$\endgroup\$
    – Snowbody
    Commented Feb 22, 2018 at 21:09

3 Answers 3


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.


As noted by greybeard in comments, double integration is available in SciPy as scipy.integrate.dblquad. This has a similar interface to the code in the post:

scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

Compute a double integral.

Return the double (definite) integral of func(y, x) from x = a..b and y = gfun(x)..hfun(x).

The differences are (i) func takes its arguments in the other order; (ii) the lower and upper boundaries must be specified by callables (but this is not a limitation because you can specify a constant boundary \$y=c\$ with the function lambda x:c); (iii) there are arguments for absolute and relative tolerance of the result.

So for the example problem you'd write:

import numpy as np

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

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

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

and then:

>>> import scipy.integrate
>>> scipy.integrate.dblquad(f, 0, 2, c, d)
(1.7777777777777706, 1.3374816809630374e-09)

(The second result is an estimate of the absolute error in the first result, though an over-estimate in this case.)

  • \$\begingroup\$ I thought that might be the problem, but as you might notice, this integral is symmetrical. f(x,y)=f(y,x). I sent my code to a friend, he ran it as is and got the right answer. I am starting to believe an exorcism might be my only way out. \$\endgroup\$ Commented Feb 23, 2018 at 10:27
  • \$\begingroup\$ So i have been using Anaconda, but tried running it in vanilla python (the exact same code) and got the right answer... \$\endgroup\$ Commented Feb 23, 2018 at 10:34

For a function strickly growing. Between a and b Recursively Search at x interssection point between the line a and b and compare the delta y if too much you you add the point c and check the intersection. After that you have to sum the trapezoid sirfaces


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.