Saya mengalami beberapa frustrasi atas cara matlab menangani integrasi numerik vs Scipy. Saya mengamati perbedaan berikut dalam kode pengujian saya di bawah ini:

1. Versi Matlab berjalan rata-rata 24 kali lebih cepat dari yang setara dengan python saya!
2. Versi Matlab mampu menghitung integral tanpa peringatan, sementara python kembali nan+nanj

Apa yang bisa saya lakukan untuk memastikan saya mendapatkan kinerja yang sama dalam python sehubungan dengan dua poin yang disebutkan? Menurut dokumentasi kedua metode harus menggunakan "quadrature adaptif global" untuk mendekati integral.

Di bawah ini adalah kode dalam dua versi (cukup mirip meskipun python mengharuskan fungsi integral dibuat sehingga dapat menangani integrand kompleks.)

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

Pertanyaannya memiliki dua subquestions yang sangat berbeda. Saya akan membahas yang pertama saja.

Versi Matlab berjalan rata-rata 24 kali lebih cepat dari yang setara dengan python saya!

Yang kedua adalah subyektif. Saya akan mengatakan bahwa memberi tahu pengguna bahwa ada beberapa masalah dengan integral adalah hal yang baik dan perilaku SciPy ini mengungguli yang Matlab `s untuk tetap diam dan entah bagaimana mencoba menghadapinya secara internal dengan cara yang hanya diketahui oleh insinyur Matlab yang memutuskan untuk menjadi yang terbaik.

Saya mengubah rentang integrasi dari 0 menjadi 30 (alih-alih dari 0 menjadi np.inf ) untuk menghindari waring NaN dan menambahkan kompilasi JIT. Untuk membandingkan solusi yang saya ulangi integrasi untuk 300 kali, hasilnya dari laptop saya.

Tanpa kompilasi JIT:

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

Dengan kompilasi JIT:

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

Dengan cara ini menambahkan dua baris kode mengarah ke faktor percepatan sekitar 40 kali kode Python dibandingkan dengan versi non-JIT. Saya tidak memiliki Matlab di laptop saya untuk memberikan perbandingan yang lebih baik, namun, jika skala dengan baik untuk PC Anda dari 24/40 = 0,6, jadi Python dengan JIT harus hampir dua kali lebih cepat dari Matlab untuk algoritma pengguna khusus ini. Kode lengkap:

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

Komentari baris @jit untuk melihat perbedaan untuk PC Anda.

Kadang-kadang fungsi untuk diintegrasikan tidak bisa JIT. Dalam hal ini, menggunakan metode integrasi lain akan menjadi solusinya.

Saya akan merekomendasikan scipy.integrate.romberg (ref) . rombergdapat mengintegrasikan fungsi yang kompleks, dan dapat mengevaluasi fungsi dengan array.