Je ressens une certaine frustration quant à la façon dont matlab gère l'intégration numérique par rapport à Scipy. J'observe les différences suivantes dans mon code de test ci-dessous:

1. La version de Matlab tourne en moyenne 24 fois plus vite que mon équivalent python!
2. La version de Matlab est capable de calculer l'intégrale sans avertissements, tandis que python retourne nan+nanj

Que puis-je faire pour m'assurer d'obtenir les mêmes performances en python par rapport aux deux points mentionnés? Selon la documentation, les deux méthodes devraient utiliser une "quadrature adaptative globale" pour approximer l'intégrale.

Vous trouverez ci-dessous le code dans les deux versions (assez similaire, bien que python nécessite la création d'une fonction intégrale pour pouvoir gérer des intégrales complexes.)

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

La question comporte deux sous-questions très différentes. Je ne traiterai que du premier.

La version de Matlab tourne en moyenne 24 fois plus vite que mon équivalent python!

Le second est subjectif. Je dirais que faire savoir à l'utilisateur qu'il y a un problème avec l'intégrale est une bonne chose et ce comportement SciPy surpasse celui de Matlab pour le garder silencieux et essayer de le traiter en interne de la manière connue uniquement par les ingénieurs de Matlab qui a décidé que c'était le meilleur.

J'ai changé la durée d'intégration pour passer de 0 à 30 (au lieu de 0 à np.inf ) pour éviter les avertissements NaN et j'ai ajouté une compilation JIT. Pour comparer la solution, j'ai répété l'intégration 300 fois, les résultats proviennent de mon ordinateur portable.

Sans compilation JIT:

$ ./test_integrate.py
34.20992112159729
(0.2618828053067563+0.24474506983644717j)

Avec la compilation JIT:

$ ./test_integrate.py
0.8560323715209961
(0.261882805306756+0.24474506983644712j)

De cette façon, l'ajout de deux lignes de code entraîne un facteur d'accélération d'environ 40 fois le code Python par rapport à une version non JIT. Je n'ai pas de Matlab sur mon ordinateur portable pour fournir une meilleure comparaison, cependant, s'il évolue bien sur votre PC que 24/40 = 0,6, donc Python avec JIT devrait être presque deux fois plus rapide que Matlab pour cet algorithme utilisateur particulier. Code complet:

#!/usr/bin/env python3
import numpy as np
from scipy import integrate
from numba import complex128,float64,jit
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)


@jit(complex128(float64, float64), nopython=True, cache=True)
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
for i in range(300): 
    #result = integral(f_integrand, 0, np.inf, omega)
    result = integral(f_integrand, 0, 30, omega)
print (time.time()-t0)
print (result)

Mettez en commentaire la ligne @jit pour voir la différence pour votre PC.

Parfois, la fonction à intégrer ne peut pas être JIT. Dans ce cas, l'utilisation d'une autre méthode d'intégration serait la solution.

Je recommanderais scipy.integrate.romberg (ref) . rombergpeut intégrer des fonctions complexes et peut évaluer la fonction avec un tableau.