Mam frustrację związaną ze sposobem, w jaki Matlab obsługuje integrację numeryczną vs. Scipy. Poniżej obserwuję następujące różnice w moim kodzie testowym:

1. Wersja Matlaba działa średnio 24 razy szybciej niż mój odpowiednik w Pythonie!
2. Wersja Matlaba jest w stanie obliczyć całkę bez ostrzeżeń, podczas gdy Python powraca nan+nanj

Co mogę zrobić, aby uzyskać taką samą wydajność w Pythonie w odniesieniu do dwóch wymienionych punktów? Zgodnie z dokumentacją obie metody powinny wykorzystywać „globalną kwadraturę adaptacyjną” w celu przybliżenia całki.

Poniżej znajduje się kod w dwóch wersjach (dość podobny, chociaż python wymaga utworzenia funkcji całkowej, aby mogła obsługiwać złożone całki).

Pyton

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

Pytanie ma dwa bardzo różne pytania. Zajmę się tylko pierwszym.

Wersja Matlaba działa średnio 24 razy szybciej niż mój odpowiednik w Pythonie!

Drugi jest subiektywny. Powiedziałbym, że poinformowanie użytkownika, że ​​jest jakiś problem z całką, jest dobrą rzeczą i to zachowanie SciPy przewyższa Matlaba, jeśli chodzi o milczenie i jakoś spróbować poradzić sobie z tym wewnętrznie w sposób znany tylko inżynierom Matlaba, którzy zdecydowałem, że będzie najlepszy.

Zmieniłem zakres integracji na z 0 na 30 (zamiast z 0 na np.inf ), aby uniknąć ostrzeżenia NaN i dodałem kompilację JIT. Aby przetestować rozwiązanie, powtórzyłem integrację 300 razy, wyniki pochodzą z mojego laptopa.

Bez kompilacji JIT:

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

Z kompilacją JIT:

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

W ten sposób dodanie dwóch wierszy kodu prowadzi do przyspieszenia około 40 razy kodu Pythona w porównaniu z wersją inną niż JIT. Nie mam Matlaba na moim laptopie, aby zapewnić lepsze porównanie, jednak jeśli skaluje się dobrze na twoim komputerze niż 24/40 = 0,6, więc Python z JIT powinien być prawie dwa razy szybszy niż Matlab dla tego konkretnego algorytmu użytkownika. Pełny kod:

#!/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)

Skomentuj linię @jit, aby zobaczyć różnicę na swoim komputerze.

Czasami funkcja integracji nie może być JIT. W takim przypadku rozwiązaniem może być zastosowanie innej metody integracji.

Polecam scipy.integrate.romberg (ref) . rombergpotrafi integrować złożone funkcje i oceniać funkcję z tablicą.