你好，我试图使用scipy dblqad函数来集成一个半球，但它不起作用

Question

我的代码是

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)

我得到的错误是

TypeError: can't convert complex to float

这是因为根是负数，所以它包含复数。

我能做些什么使它正常工作吗？

非常感谢。

Answer 1

使用math.sqrt或np.sqrt而不是cmath.sqrt：即使这些复数的虚部总是零，integrate.dblquad也无法处理复数，在数学上，在这里使用复数没有什么意义。

为了避免取负数的平方根的潜在问题，您可以：

• 调整被积函数，使其在感兴趣区域之外为零，或者
• 调整限制，以便在单位磁盘上积分，而不是在正方形上积分

以下是执行前者的代码版本。注意使用np.maximum和np.sqrt，而不是Python内置的max和math.sqrt：这确保了f适用于数组参数和标量参数，这可能允许dblquad更快地计算其结果。

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)

对我来说，这产生了预期的2/3 π:近似值

(2.0943951045290796, 1.855738140932317e-08)

这里有一个调整限制的代码版本：

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)