Hello, I was trying to integrate a hemisphere using scipy dblquad function, and its not working

Question

my code is

import numpy as np
from scipy import integrate
from math import *
import cmath


f = lambda y, x: cmath.sqrt(1 - x**2 - y**2)
hemisphere = integrate.dblquad(f, -1, 1, lambda x: -1, lambda x: 1)

print(hemisphere)

and the error i get is

TypeError: can't convert complex to float

It is because the root is negative, so it contains complex numbers.

Is there something I can do so it works properly?

Thanks a lot.

Answer 1

Use math.sqrt or np.sqrt instead of cmath.sqrt: integrate.dblquad can't handle complex numbers, even if the imaginary part of those complex numbers is always zero, and mathematically it makes little sense to be using complex numbers here.

To avoid the potential issue of taking the square root of a negative number, you can either:

• adjust the integrand so that it's zero outside the region of interest, or
• adjust the limits so that you're integrating over the unit disk instead of over a square

Here's a version of your code that does the former. Note the use of np.maximum and np.sqrt instead of the Python builtin max and math.sqrt: this ensures that f works for array arguments as well as scalar arguments, which may allow dblquad to compute its result faster.

import numpy as np
from scipy import integrate

f = lambda y, x: np.sqrt(np.maximum(0.0, 1 - x**2 - y**2))
hemisphere = integrate.dblquad(f, -1, 1, lambda x: -1, lambda x: 1)

print(hemisphere)

For me, this produces the expected approximation to 2/3 π:

(2.0943951045290796, 1.855738140932317e-08)

And here's a version of your code that adjusts the limits instead:

import numpy as np
from scipy import integrate

f = lambda y, x: np.sqrt(1 - x**2 - y**2)
hemisphere = integrate.dblquad(f, -1, 1, lambda x: -np.sqrt(1-x*x),
                                         lambda x: np.sqrt(1-x*x))

print(hemisphere)