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

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

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

• 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. – Martin Johnsrud Feb 23 '18 at 10:27
• So i have been using Anaconda, but tried running it in vanilla python (the exact same code) and got the right answer... – Martin Johnsrud Feb 23 '18 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