Python quad2和MATLAB integral2返回不同的结果

我一直在用Python重写一些MATLAB代码，在尝试用Python执行二重积分时遇到了一些问题。 下面的MATLAB代码为Prad返回一个值2.5133

               ant.Frequency = 3e8;
               N = 5;
               lambda = 3e8/ant.Frequency;
               d = 0.5*lambda;
               k = 2*pi/lambda;

               theta = linspace(0,180,180)*pi/180;
               phi = linspace(0,360,180) * pi/180;

               theta0 = 0*pi/180;
               phi0 = 0*pi/180;

               PsiX = (k.*d.*sin(theta).*cos(phi) - k*d*sin(theta0)*cos(phi0));
               PsiY = (k.*d.*sin(theta).*sin(phi) - k*d*sin(theta0)*sin(phi0));
               AF = ((sin(N.*PsiX./2)./(N.*sin(PsiX./2)))).*((sin(N.*PsiY./2)./(N.*sin(PsiY./2))));

               Pfun = @(x,y) ((1/N.*(sin(N.*(k.*d.*sin(x).*cos(y) - k.*d.*sin(theta0).*cos(phi0))./2)...
                   ./sin((k.*d.*sin(x).*cos(y) - k.*d.*sin(theta0).*cos(phi0))./2))).^2).*sin(x);

               Prad = integral2(Pfun,0,pi,0,2*pi);

我尝试在Python中创建相同的行为，结果Prad值为3.209E-16

    theta = np.linspace(0,180,180)*np.pi/180
    phi = np.linspace(0,360,180)*np.pi/180
    theta0 = np.array([0])
    phi0 = np.array([0])

    wavelength = 3E8/3E8
    k = 2 * np.pi/wavelength
    d = 0.5 * wavelength
    N = 5

    PsiX = (k*d*np.sin(theta)*np.cos(phi) - k*d*np.sin(theta0)*np.cos(phi0))
    PsiY = (k*d*np.sin(theta)*np.sin(phi) - k*d*np.sin(theta0)*np.sin(phi0))

    AF = ((np.sin(N*PsiX/2)/(N*np.sin(PsiX/2)))) * ((np.sin(N*PsiY/2)/(N*np.sin(PsiY/2))))

    Pfun = lambda x,y: ((1/N*(np.sin(N*(k*d*np.sin(x)*np.cos(y) - k*d*np.sin(theta0)*np.cos(phi0))/2)/np.sin((k*d*np.sin(x)*np.cos(y) - k*d*np.sin(theta0)*np.cos(phi0))/2)))**2)*np.sin(x)
    Prad = dblquad(Pfun,0,np.pi,lambda x: 0, lambda x: 2*np.pi)

即使使用MATLAB dblquad函数也可以得到期望的结果。我在Python代码中有语法错误吗?还是我误解了某些功能? 问题来源StackOverflow 地址：/questions/59381403/python-quad2-and-matlab-integral2-returning-different-results