尝试在 Python 中实现 Matlab 代码时出现问题

Question

我在 Matlab 中有一段运行良好的代码（见下文）：

clc; clear all;
f = 8500;
c0 = 343;
rho = 1.225;
omega = 2*pi*f;
k = omega/c0;
Z = -426;

lx = 0.1;
ly = 0.1;
nx = 50;
ny = nx/2;

integrand1 = @(x,y,kx) real(((exp(1i*(kx*x + sqrt(k.^2 - kx.^2).*abs(y))))./sqrt(k.^2 - kx.^2)));
integrand2 = @(x,y,kx) imag(((exp(1i*(kx*x + sqrt(k.^2 - kx.^2).*abs(y))))./sqrt(k.^2 - kx.^2)));
Gz = @(x,y) integral(@(kx)integrand1(x,y,kx), -100*k, 100*k, "AbsTol", 1e-4) + 1i*...
            integral(@(kx)integrand2(x,y,kx), -100*k, 100*k, "AbsTol", 1e-4);
Test = Gz(lx,ly)

但我需要在 Python 中实现它。我尝试使用与 lambda 表达式相同的方法将函数定义为函数句柄，但不幸的是，它不起作用并且出现错误。这些是我的 Python 代码：

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt


f = 8500       
rho = 1.225        
c0 = 343          
omega = 2*np.pi*f   
k = omega/c0       
Z = -426           

Lx = 0.1                       
Ly = 0.1                       
nx = 50                        
ny = int(nx/2)                  


integrandReal = lambda kx, x, y: np.real(((2*np.sqrt(k**2 - kx**2)*Z)/(np.sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + np.sqrt(k**2 - kx**2)*y)))/(np.sqrt(k**2 - kx**2))))
integrandImag = lambda kx, x, y: np.imag(((2*np.sqrt(k**2 - kx**2)*Z)/(np.sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + np.sqrt(k**2 - kx**2)*y)))/(np.sqrt(k**2 - kx**2))))

integral = lambda x, y: integrate.quad(integrandReal, -100*k, 100*k, args=(x,y)) + 1j*integrate.quad(integrandImag, -100*k, 100*k, args=(x,y))

G = integral(Lx,Ly)

这些是错误：

runfile(XYZ/'untitled0.py', wdir='/Users/XYZ ')
/Users/XYZ/untitled0.py:29: RuntimeWarning: invalid value encountered in sqrt
  integrandReal = lambda kx, x, y: np.real(((2*np.sqrt(k**2 - kx**2)*Z)/(np.sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + np.sqrt(k**2 - kx**2)*y)))/(np.sqrt(k**2 - kx**2))))
/Users/ /untitled0.py:32: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  integral = lambda x, y: integrate.quad(integrandReal, -100*k, 100*k, args=(x,y)) + 1j*integrate.quad(integrandImag, -100*k, 100*k, args=(x,y))
/Users/XYZ/untitled0.py:30: RuntimeWarning: invalid value encountered in sqrt
  integrandImag = lambda kx, x, y: np.imag(((2*np.sqrt(k**2 - kx**2)*Z)/(np.sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + np.sqrt(k**2 - kx**2)*y)))/(np.sqrt(k**2 - kx**2))))
Traceback (most recent call last):

  File "/Users/XYZ/untitled0.py", line 34, in <module>
    G = integral(0.1,0.1)

  File "/Users/XYZ/untitled0.py", line 32, in <lambda>
    integral = lambda x, y: integrate.quad(integrandReal, -100*k, 100*k, args=(x,y)) + 1j*integrate.quad(integrandImag, -100*k, 100*k, args=(x,y))

TypeError: can't multiply sequence by non-int of type 'complex'

非常欢迎任何帮助。 先感谢您！

Answer 1

scipy.integrate.quad 返回一个元组，不是数字，因此您应该选择其中的第一个返回值，即数值结果。

并且使用 np.sqrt 将在负输入上返回 nan ，您应该使用 numpy.lib.scimath.sqrt 因为它可以处理负输入并返回复数，如下所示。

from numpy.lib.scimath import sqrt
integrandReal = lambda kx, x, y: np.real(((2*sqrt(k**2 - kx**2)*Z)/(sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + sqrt(k**2 - kx**2)*y)))/(sqrt(k**2 - kx**2))))
integrandImag = lambda kx, x, y: np.imag(((2*sqrt(k**2 - kx**2)*Z)/(sqrt(k**2 - kx**2)*Z + omega*rho))*((np.exp(1j*(kx*x + sqrt(k**2 - kx**2)*y)))/(sqrt(k**2 - kx**2))))

integral = lambda x, y: integrate.quad(integrandReal, -100*k, 100*k, args=(x,y))[0] + 1j*integrate.quad(integrandImag, -100*k, 100*k, args=(x,y))[0]