Scipy 虚数的四次积分

Question

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

其中 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=",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)

这种方法的输出是：

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, 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)

错误在第 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",z)
    return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]

进入：

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]

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

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

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

Answer 1

这是我对你积分的解读

关于如何推广积分的答案可以扩展到二重积分（由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, 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:])

然后你计算你的积分如下

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

Scipy 虚数的四次积分

Scipy quad integral of imaginary numbers

python

numpy

scipy

quad