我对 matlab 处理数值积分与 Scipy 的方式感到有些沮丧。我在下面的测试代码中观察到以下差异:
- Matlab 版本的运行速度平均比我的 python版本快 24 倍!
- Matlab 的版本能够在没有警告的情况下计算积分,而 python 返回
nan+nanj
对于提到的两点,我能做些什么来确保我在 python 中获得相同的性能?根据文档,这两种方法都应该使用“全局自适应正交”来近似积分。
下面是两个版本的代码(虽然python要求创建一个积分函数以便它可以处理复杂的被积函数,但非常相似。)
Python
import numpy as np
from scipy import integrate
import time
def integral(integrand, a, b, arg):
def real_func(x,arg):
return np.real(integrand(x,arg))
def imag_func(x,arg):
return np.imag(integrand(x,arg))
real_integral = integrate.quad(real_func, a, b, args=(arg))
imag_integral = integrate.quad(imag_func, a, b, args=(arg))
return real_integral[0] + 1j*imag_integral[0]
vintegral = np.vectorize(integral)
def f_integrand(s, omega):
sigma = np.pi/(np.pi+2)
xs = np.exp(-np.pi*s/(2*sigma))
x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+np.sqrt(1-xs**2))
x2 = 1-2*sigma/np.pi*(1-xs)
zeta = x2+x1*1j
Vc = 1/(2*sigma)
theta = -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
t1 = 1/np.sqrt(1+np.tan(theta)**2)
t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
return np.real((t1-1j*t2)/np.sqrt(zeta**2-1))*np.exp(1j*omega*s/Vc);
t0 = time.time()
omega = 10
result = integral(f_integrand, 0, np.inf, omega)
print time.time()-t0
print result
MATLAB
function [ out ] = f_integrand( s, omega )
sigma = pi/(pi+2);
xs = exp(-pi.*s./(2*sigma));
x1 = -2*sigma./pi.*(log(xs./(1+sqrt(1-xs.^2)))+sqrt(1-xs.^2));
x2 = 1-2*sigma./pi.*(1-xs);
zeta = x2+x1*1j;
Vc = 1/(2*sigma);
theta = -1*asin(exp(-pi./(2.0.*sigma).*s));
t1 = 1./sqrt(1+tan(theta).^2);
t2 = -1./sqrt(1+1./tan(theta).^2);
out = real((t1-1j.*t2)./sqrt(zeta.^2-1)).*exp(1j.*omega.*s./Vc);
end
t=cputime;
omega = 10;
result = integral(@(s) f_integrand(s,omega),0,Inf)
time_taken = cputime-t