python - Utilice scipy.integrate.quad para integrar números complejos

numpy complex-numbers

Question

Estoy usando en este momento el scipy.integrate.quad para integrar con éxito algunos integrands reales. Ahora apareció una situación que necesito para integrar un integrand complejo. Quad parece no poder hacerlo, como las otras rutinas scipy.integrate, así que pregunto: ¿hay alguna forma de integrar un integrando complejo usando scipy.integrate, sin tener que separar la integral en las partes real e imaginaria?

Answer 1

Me doy cuenta de que llego tarde a la fiesta, pero quizás quadpy (un proyecto mío) pueda ayudar. Esta

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)

correctamente da

(0.841470984808+0.459697694132j)

En lugar de GaussKronrod(3) , puedes usar cualquier otro esquema.

Answer 2

¿Qué hay de malo en separarlo en partes reales e imaginarias? scipy.integrate.quad requiere la función integrada de devoluciones flotantes (también conocidos como números reales) para el algoritmo que utiliza.

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

P.ej,

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

que es lo que esperas al error de redondeo - integral de exp (ix) de 0, pi / 2 es (1 / i) (e ^ i pi / 2 - e ^ 0) = -i (i - 1) = 1 + i ~ (0.99999999999999989 + 0.999999999999999998j).

Y para el registro en caso de que no sea 100% claro para todos, la integración es lineal funcional, lo que significa que ∫ {f (x) + kg (x)} dx = ∫ f (x) dx + k ∫ g (x ) dx (donde k es una constante con respecto a x). O para nuestro caso específico ∫ z (x) dx = ∫ Re z (x) dx + i ∫ Im z (x) dx como z (x) = Re z (x) + i Im z (x).

Si está tratando de hacer una integración sobre una ruta en el plano complejo (que no sea a lo largo del eje real) o la región en el plano complejo, necesitará un algoritmo más sofisticado.

Nota: Scipy.integrate no manejará directamente la integración compleja. ¿Por qué? Realiza el trabajo pesado en la biblioteca FORTRAN QUADPACK , específicamente en qagse.f que requiere explícitamente que las funciones / variables sean reales antes de hacer su "cuadratura adaptativa global basada en la cuadratura de Gauss-Kronrod de 21 puntos dentro de cada subintervalo, con la aceleración de Peter". El algoritmo épsilon de Wynn ". Por lo tanto, a menos que desee probar y modificar el FORTRAN subyacente para que maneje números complejos, compílelo en una nueva biblioteca, no logrará que funcione.

Si realmente desea hacer el método de Gauss-Kronrod con números complejos en exactamente una integración, consulte la página de wikipedias e implemente directamente como se hace a continuación (usando la regla de 15 puntos, 7 puntos). Tenga en cuenta que la función memoize''d para repetir llamadas comunes a las variables comunes (suponiendo que las llamadas a funciones sean lentas como si la función fuera muy compleja). También cumplí con las reglas de 7 y 15 puntos, ya que no tenía ganas de calcular los nodos / pesos yo mismo y esos eran los que figuran en wikipedia, pero obtenía errores razonables para los casos de prueba (~ 1e-14)

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]

Caso de prueba:

>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0)
[(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19]

No confío en la estimación de error. Tomé algo de wiki para la estimación de error recomendada al integrarme de [-1 a 1] y los valores no me parecen razonables. Por ejemplo, el error anterior comparado con la verdad es ~ 5e-15 no ~ 1e-19. Estoy seguro de que si alguien consultara varias recetas, podría obtener una estimación más precisa. (Probablemente tenga que multiplicar por (ab)/2 a algún poder o algo similar).

Recuerde, la versión de python es menos precisa que simplemente llamar dos veces a la integración basada en QUADPACK de scipy. (Podrías mejorarlo si lo deseas).