Numeric double integration

Question

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

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

Answer 1

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.

Answer 2

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

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

Answer 3

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