我现在正在使用scipy.integrate.quad成功地积分一些实值积分。现在出现了需要积分复合函数的情况。quad不能够处理它，其他的scipy.integrate例程也不能处理它。因此，我想问：是否有一种方法使用scipy.integrate积分一个复合函数，而不必将其分离成实部和虚部？
2个回答
60
为什么不能将其分解为实部和虚部？scipy.integrate.quad需要被积函数返回浮点数（也就是实数），以便使用它的算法。
import scipy
from scipy.integrate import quad
def complex_quadrature(func, a, b, **kwargs):
def real_func(x):
return scipy.real(func(x))
def imag_func(x):
return scipy.imag(func(x))
real_integral = quad(real_func, a, b, **kwargs)
imag_integral = quad(imag_func, a, b, **kwargs)
return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
E.g.,
>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2)
((0.99999999999999989+0.99999999999999989j),
(1.1102230246251564e-14,),
(1.1102230246251564e-14,))
这是您期望的舍入误差 - 从0到pi/2积分exp(i x)等于(1/i)(e^i pi/2 - e^0) = -i(i - 1) = 1 + i ~ (0.99999999999999989+0.99999999999999989j)。
顺便提一下，如果对每个人来说不是100％清楚，那么积分是一个线性函数，意味着∫ {f(x) + k g(x)} dx = ∫ f(x) dx + k ∫ g(x) dx（其中k是相对于x的常数）。或者对于我们的特定情况，∫ z(x) dx = ∫ Re z(x) dx + i ∫ Im z(x) dx，因为z(x) = Re z(x) + i Im z(x)。
如果您要在复平面上的路径（除了沿着实轴）或复平面上的区域上进行积分，则需要更复杂的算法。
注意：Scipy.integrate不会直接处理复合积分。为什么？它在FORTRAN QUADPACK库中完成了繁重的工作，具体而言，在qagse.f中，它明确要求在执行其“基于21点Gauss-Kronrod积分的全局自适应积分”的计算之前将函数/变量设为实数。除非您想尝试修改底层FORTRAN以处理复数，将其编译成新库，否则您将无法使其工作。
如果你真的想在一次积分中使用复数执行高斯-克朗罗德方法，请查看维基百科页面并直接实现以下内容（使用15点和7点规则）。请注意，我将函数记忆化以重复常见调用到公共变量（假设函数调用很慢，就像函数非常复杂一样）。由于我不想自己计算节点/权重，并且那些是维基百科上列出的节点/权重，但对于测试用例（~1e-14）获得了合理的误差，因此仅执行了7点和15点规则。import scipy
from scipy import array
def quad_routine(func, a, b, x_list, w_list):
c_1 = (b-a)/2.0
c_2 = (b+a)/2.0
eval_points = map(lambda x: c_1*x+c_2, x_list)
func_evals = map(func, eval_points)
return c_1 * sum(array(func_evals) * array(w_list))
def quad_gauss_7(func, a, b):
x_gauss = [-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759]
w_gauss = array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
return quad_routine(func,a,b,x_gauss, w_gauss)
def quad_kronrod_15(func, a, b):
x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]
return quad_routine(func,a,b,x_kr, w_kr)
class Memoize(object):
def __init__(self, func):
self.func = func
self.eval_points = {}
def __call__(self, *args):
if args not in self.eval_points:
self.eval_points[args] = self.func(*args)
return self.eval_points[args]
def quad(func,a,b):
''' Output is the 15 point estimate; and the estimated error '''
func = Memoize(func) # Memoize function to skip repeated function calls.
g7 = quad_gauss_7(func,a,b)
k15 = quad_kronrod_15(func,a,b)
# I don't have much faith in this error estimate taken from wikipedia
# without incorporating how it should scale with changing limits
return [k15, (200*scipy.absolute(g7-k15))**1.5]
测试用例：
>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0)
[(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19]
我不信任误差估计 - 我从维基百科中获取了从[-1到1]积分时的推荐误差估计，但这些值对我来说似乎不合理。例如，与真实值相比，误差约为5e-15而不是1e-19。我相信如果有人参考num recipes，就可以得到更准确的估计。（可能需要将一个减去另一个除以2的结果乘以某个幂或类似的东西）。
请记住，Python版本的准确性不如调用scipy的QUADPACK两次进行积分。（如果需要，您可以改进它）。
- dr jimbob
4
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
**作为前缀可以将其变成一个字典，包含其他可选的命名参数，你可以一并传递。请参见https://docs.python.org/2/faq/programming.html#how-can-i-pass-optional-or-keyword-parameters-from-one-function-to-another因此，您可以使用可选参数调用complex_quadrature以传递给quad。some_sliceable [1：]切片可切片对象的零号元素之后的所有内容。通常，它只包含估计的绝对误差。但是，如果您调用complex_quadrature（f，a，b，full_output = True），则会返回额外的信息。 - dr jimbob