Scipy quad integral of imaginary numbers

Question

I want to calculate integral of implicit function containing imaginary numbers

where f(iz) is something like:

and g(ix) is something like:

I want to calculate it numerically. Python scipy.quad doesn't calculate integrals of imaginary numbers (explained in Code 1 below). Quadpy isn't efficient also, because it passes entire numpy array instead of single values of integral (explained in code 2 below) and thus needs additional manipulation. So I am thinking about dividing integrals like in the way shown below (where Re is real part and Im is imaginary part):

and expanding above equation:

Can I do that?

And here are codes, where I show two ways to approach the problem.

Code 1:

First approach scipy.quad. According to: Use scipy.integrate.quad to integrate complex numbers I tried dividing into real and imaginary parts and calculating integral values of them separately:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def f_D(x, z):
    print("f_D:",z)
    return 1*np.exp(1j*x)*z

def phi_0(z):
    print("phi:",z)
    return 2/(z*M_r(z))

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D",z)
    return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]

def integral_P_Vm(z):
    print("Calling function, z=",z)
    return M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

def integral_P_Vm_real(z):
    print("Calling real function, z=",z)
    # return np.real(2*np.exp(1j*(z))*jv(m*B, (B*z))*2)
    return np.real(M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z))

def integral_P_Vm_imag(z):
    print("Calling imag function, z=",z)
    # return np.imag(2*np.exp(1j*(z))*jv(m*B, (B*z))*2)
    return np.imag(M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z))


Quad_integral = quad(integral_P_Vm, 0, 1, limit=100)
P_Vm_integral_real = quad(integral_P_Vm_real, 0, 1, limit=100)
P_Vm_integral_imag = quad(integral_P_Vm_imag, 0, 1, limit=100)

print("Quad",Quad_integral)
print("Real part:",P_Vm_integral_real)
print("Imaginary part",P_Vm_integral_imag)

The output of this approach is:

Quad (-5.63500792439988e-12, 1.4665732299181548e-21)
Real part: (-5.63500792439988e-12, 1.4665732299181548e-21)
Imaginary part (9.08758050641791e-12, 3.696603118432144e-22)

As you can see, normal quad omits imaginary part and thus I am not sure if integral in def psi_D(z): also needs to be divided into imaginary and real parts.

Code 2:

Second approach - Quadpy.quad the code is:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def f_D(x, z):
    print("f_D:",z)
    return 1*np.exp(1j*x)*z

def phi_0(z):
    print("phi:",z)
    return 2/(z*M_r(z))

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D",z)
    return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]

def integral_P_Vm(z):
    print("Calling function, z=",z)
    return M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

Quadpy_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

print("Quadpy",Quadpy_integral)

and the error is in line 28:

Exception has occurred: TypeError
only size-1 arrays can be converted to Python scalars

I tried to overcome this error by iterating numpy array z by modyfing:

def psi_D(z):
    print("psi_D",z)
    return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]

into:

def psi_D(z):
    print("psi_D",z)
    for e in z:
        print(e)
        return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]

But the iteration stops at first element in z array and I have an error:

Exception has occurred: TypeError
f_D() argument after * must be an iterable, not numpy.float64

I have no idea how to iterate over entire numpy array in this case.

It's been a long time since I took a complex math class. Just what does the complex integral mean? Same as treating the real and imaginary parts as 2dimensions, in effect a 2d integral on `x` and `y`, or are they more tightly bound? `np.quad` is a single value integral. There is a `dblquad` — hpaulj, May 27 '21 at 00:32
did you take a look at https://stackoverflow.com/questions/5965583/use-scipy-integrate-quad-to-integrate-complex-numbers? — yann ziselman, May 27 '21 at 05:44

score 1 · Accepted Answer · answered May 27 '21 at 09:21

This is my interpretation of your integral

The answer about how to generalize the integral can be extended to double integral (provided by dblquad function).

import numpy as np
from scipy.integrate import dblquad

def complex_dblquad(func, a, b, g, h, **kwargs):
    def real_func(z, x):
        return np.real(func(z, x))
    def imag_func(z, x):
        return np.imag(func(z, x))
    real_integral = dblquad(real_func, a, b, g, h, **kwargs)
    imag_integral = dblquad(imag_func, a, b, g, h, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

Then you compute your integral as follows

complex_dblquad(lambda z,x: np.exp(1j*x*z), 0, 1, 0.3, 0.9)

Thank you! Using your answer, I wrote the code and used following equation: `complex_integral = complex_dblquad(lambda z,x: np.sqrt(1+z*z)**2*(2/np.sqrt(1+z*z))**2*H_interpolation(x)*np.exp(1j*2/np.sqrt(1+z*z)*x), 0, 1, -0.5, 0.5)` Is this proper implementation of this set of equations? https://i.stack.imgur.com/8d4VF.png More details: https://math.stackexchange.com/questions/4166396/did-i-properly-implement-double-integral-equation-in-dblquad-python-helicoidal — Anna Majewska, Jun 07 '21 at 21:43

Scipy quad integral of imaginary numbers

1 Answers1

Linked