scipy.integrate.quad를 사용하여 복소수를 통합하십시오.

Question

저는 실제 적분을 성공적으로 통합하기 위해 scipy.integrate.quad를 지금 사용하고 있습니다. 복잡한 통합을 통합해야하는 상황이 나타났습니다. 쿼드는 다른 scipy.integrate 루틴처럼 그렇게 할 수 없을 것 같습니다. scipy.integrate를 사용하여 복잡한 integrand를 통합 할 수있는 방법이 있습니까? 실제와 가상 부분의 적분을 분리하지 않아도됩니까?

Answer 1

실재와 허수 부로 구분하는 것이 잘못된 이유는 무엇입니까? scipy.integrate.quad에는 사용하는 알고리즘에 대해 통합 함수 return float (실수로도 함)가 필요합니다. 당신이 오류를 반올림 기대하는 것입니다

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

예,

>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2) 
((0.99999999999999989+0.99999999999999989j), 
(1.1102230246251564e-14,), 
(1.1102230246251564e-14,)) 

- 특급 (IX)의 통합 0이, 파이/2 (1/I) (예^나는 파이/2 - e^0) = -i (i - 1) = 1 + i - (0.99999999999999989 + 0.99999999999999989j).

그리고 모든 사람에게 100 % 명확하지 않은 경우의 기록은 일차 함수 적입니다. 즉, ∫ {f (x) + kg (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 구적법을 기반으로 한 "전역 적응 형 직교 (global adaptive quadrature)"를 수행하기 전에 함수/변수를 명시 적으로 요구하며, Peter Wynn의 엡실론 연산." 따라서, 기본 FORTRAN을 수정하여 복소수를 처리하고 새로운 라이브러리로 컴파일하면 작동하지 않을 것입니다.

정확히 하나의 통합으로 복소수를 가진 Gauss-Kronrod 방법을 수행하려면 wikipedias page을보고 아래에 설명 된대로 직접 구현하십시오 (15pt, 7pt 규칙 사용). 참고, 나는 일반적인 변수에 대한 일반적인 호출을 반복하는 함수를 memoize했습니다 (함수 호출이 매우 복잡한 것처럼 느린 경우). 또한 노드와 가중치를 직접 계산하지 않고 위키 피 디아에 나열된 것들 이었지만 테스트 케이스 (~ 1e-14)에 대해 합당한 오류가 발생했다고 생각했기 때문에 7-pt 및 15-pt 규칙 만 수행했습니다.

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를 참고하면 더 정확한 견적을 얻을 수있을 것입니다. (어쩌면 (a-b)/2에 의해 몇 가지 힘이나 비슷한 것에 배수해야합니다).

파이썬 버전은 scipy의 QUADPACK 기반 통합을 두 번 호출하는 것보다 덜 정확합니다. (원하는 경우 개선 할 수 있습니다.)

Answer 2

나는 파티에 늦었다는 것을 알고 있지만 아마도 quadpy (내 프로젝트)이 도움이 될 수 있습니다.이

import numpy 
import quadpy 
import scipy 

val = quadpy.line_segment.integrate(
     lambda x: scipy.exp(1j*x), 
     [0, 1], 
     quadpy.line_segment.GaussKronrod(3) 
     ) 

print(val) 

이 제대로 당신은 다른 방식을 사용할 수

(0.841470984808+0.459697694132j) 

대신 GaussKronrod(3)

을 제공합니다.