尝试在 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 return 是一个元组，不是数字，所以你应该 select 它的第一个 return 是数字结果。

并使用 np.sqrt 将 return nan 用于负输入，您应该改用 numpy.lib.scimath.sqrt 因为它可以处理负输入和 return 复数，如下所示。

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]

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

Problem trying to implement Matlab code in Python

python

numpy

matlab

scipy

quad