Scipy 虚数四重积分

如何解决Scipy 虚数四重积分

我想计算包含虚数的隐函数的积分

其中 f(iz) 类似于：

和 g(ix) 类似于：

我想用数字来计算它。 Python scipy.quad 不计算虚数的积分（在下面的代码 1 中解释）。 Quadpy 也效率不高，因为它传递整个 numpy 数组而不是单个整数值（在下面的代码 2 中解释），因此需要额外的操作。所以我正在考虑以如下所示的方式划分积分（其中 Re 是实部，Im 是虚部）：

并展开上面的等式：

我可以这样做吗？

这里是代码，我展示了两种解决问题的方法。

代码 1：

第一种方法scipy.quad。根据： Use scipy.integrate.quad to integrate complex numbers 我尝试将实部和虚部分开并分别计算它们的整数值：

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)
    # 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)
    # 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,1,limit=100)
P_Vm_integral_real = quad(integral_P_Vm_real,limit=100)
P_Vm_integral_imag = quad(integral_P_Vm_imag,limit=100)

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

这种方法的输出是：

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

如您所见，普通四边形省略了虚部，因此我不确定 def psi_D(z): 中的积分是否也需要分为虚部和实部。

代码 2：

第二种方法 - Quadpy.quad 代码是：

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,((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

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

print("Quadpy",Quadpy_integral)

错误在第 28 行：

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

我试图通过 modyfing 迭代 numpy 数组 z 来克服这个错误：

def psi_D(z):
    print("psi_D",args=(z))[0]

进入：

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

但是迭代在 z 数组中的第一个元素处停止并且我有一个错误：

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

在这种情况下，我不知道如何遍历整个 numpy 数组。

解决方法

这是我对你积分的解释

关于如何推广积分的答案可以推广到二重积分（由dblquad函数提供）。

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,**kwargs)
    imag_integral = dblquad(imag_func,**kwargs)
    return (real_integral[0] + 1j*imag_integral[0],real_integral[1:],imag_integral[1:])

然后你计算你的积分如下

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