Problem trying to implement Matlab code in Python-CodePudding

I have a piece of code in Matlab (see below) that works well:

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)

But I need to implement it in Python. I tried to use the same approach with lambda expressions to define my functions as functions handle, but unfortunately, it's not working and I am getting errors. These is my Python code:

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)

These are the errors:

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'

Any help is more than welcomed. Thank you in advance!

CodePudding user response：

scipy.integrate.quad returns a tuple, not a number, so you should select the first return from it which is the numerical result.

and using np.sqrt will return nan on negative input, you should instead use numpy.lib.scimath.sqrt as it can handle negative inputs and return complex numbers as shown below.

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]