1

Here is an integration routine I found on another question see Use scipy.integrate.quad to integrate complex numbers. I am happy with its functionality but I wish to be able to perform the integration on my input function for different values of the parameter 't' and to do so using vectorization. 

import numpy as np
from scipy.integrate import quad

def complex_quadrature(func, a, b, t):
    def real_func(x,t):
        return np.real(func(x,t))
    def imag_func(x,t):
        return np.imag(func(x,t))
    real_integral = quad(real_func, a, b, args=(t))
    imag_integral = quad(imag_func, a, b, args=(t))
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


vcomplex_quadrature = np.vectorize(complex_quadrature)


t = np.linspace(1,4,4)

print(vcomplex_quadrature(lambda x, t: (np.exp(1j*x*t)), 0, np.pi/2, t))

Can anyone provide some help with this?

Community
  • 1
  • 1
Dipole
  • 1,840
  • 2
  • 24
  • 35

1 Answers1

0

You want to do a Fourier Transform. Use numpy.fft or scipy.fft instead.

For a more general function, there is no completely good way. np.vectorize is just a convenience function with a for loop under the hood, so you will get no performance gain here.

Any integral can be expressed as a ODE. You could then express all the independent integrals as a set of ODEs and feed them to the solver in Scipy, but this will damage accuracy because the solver uses an adaptive step. When you have several equations, the choice of the step will be the same for all of them, and estimated on the general behaviour, instead of a "personalised" individual step.

The good news are that this problem is embarrassingly parallelisable, which means you can write as a parallel map and use all the cores in your machine. This will give you a much better performance boost than any numpy vectorisation.

For simple parallelisation in Python, use joblib

Davidmh
  • 3,797
  • 18
  • 35