Jak wykonać dwie integracje numeryczne zmiennych w Julia?

Question

Potrafię wykonywać integrację numeryczną jednej zmiennej w Julii przy użyciu quadgk. Kilka prostych przykładów:

julia> f(x) = cos(x) 
f (generic function with 1 method) 

julia> quadgk(f, 0, pi) 
(8.326672684688674e-17,0.0) 

julia> quadgk(f, 0, pi/2) 
(1.0,1.1102230246251565e-16) 

julia> g(x) = cos(x)^2 
g (generic function with 1 method) 

julia> quadgk(g, 0, pi/2) 
(0.7853981633974483,0.0) 

julia> pi/4 
0.7853981633974483 

documentation for quadgk nie wydaje się sugerować poparcie dla wielowymiarowej integracji, a na pewno wystarczy otrzymuję komunikat o błędzie, jeśli próbują nadużywać go za integralną 2D:

julia> quadgk(h, 0, pi/2, 0, pi/2) 
ERROR: `h` has no method matching h(::Float64) 

Dokumentacja sugeruje, że istnieją pewne pakiety zewnętrzne do integracji, ale ich nie wymienia. Zgaduję, że jeden taki pakiet może zrobić dwuwymiarowe całki. Jaki jest najlepszy taki pakiet do tego zadania?

Answer 1

myślę, że chcesz sprawdzić pakiet Kubatura:

https://github.com/stevengj/Cubature.jl

Prawdopodobnie quadgk powinny być po prostu usunięte z biblioteki standardowej, ponieważ jest ograniczona, a tylko wprowadza w błąd ludzi do nie szuka pakietu zrobić integrację.

Answer 2

Oprócz Cubature.jl jest jeszcze inny pakiet Julia, która pozwala na obliczenie całki wielowymiarowe numeryczne: Cuba.jl (https://github.com/giordano/Cuba.jl). Można ją zainstalować za pomocą menedżera pakietów:

Pkg.add("Cuba") 

Kompletna dokumentacja pakietu jest dostępna w https://cubajl.readthedocs.org (również w PDF version)

Disclaimer: Jestem autor pakietu.

---

Cuba.jl jest po prostu Julia otoki wokół Cuba Library, Thomas Hahn i oferuje cztery niezależne algorytmy do obliczania całek: Vegas, Suave, pomiędzy polem Cuhre.

Całka cos (x) w obszarze [0, 1], może być obliczony z jednym z następujących poleceń:

Vegas((x,f)->f[1]=cos(x[1]), 1, 1) 
Suave((x,f)->f[1]=cos(x[1]), 1, 1) 
Divonne((x,f)->f[1]=cos(x[1]), 1, 1) 
Cuhre((x,f)->f[1]=cos(x[1]), 1, 1) 

As a more advanced example całka

gdzie Ω = [0, 1] ³ i

można obliczyć za pomocą następującego scenariusza: Julia

using Cuba 

function integrand(x, f) 
    f[1] = sin(x[1])*cos(x[2])*exp(x[3]) 
    f[2] = exp(-(x[1]^2 + x[2]^2 + x[3]^2)) 
    f[3] = 1/(1 - x[1]*x[2]*x[3]) 
end 

result = Cuhre(integrand, 3, 3, epsabs=1e-12, epsrel=1e-10) 
answer = [(e-1)*(1-cos(1))*sin(1), (sqrt(pi)*erf(1)/2)^3, zeta(3)] 
for i = 1:3 
    println("Component $i") 
    println(" Result of Cuba: ", result[1][i], " ± ", result[2][i]) 
    println(" Exact result: ", answer[i]) 
    println(" Actual error: ", abs(result[1][i] - answer[i])) 
end 

co daje następujący wynik

Component 1 
Result of Cuba: 0.6646696797813739 ± 1.0050367631018485e-13 
Exact result: 0.6646696797813771 
Actual error: 3.219646771412954e-15 
Component 2 
Result of Cuba: 0.4165383858806454 ± 2.932866749838454e-11 
Exact result: 0.41653838588663805 
Actual error: 5.9926508200192075e-12 
Component 3 
Result of Cuba: 1.2020569031649702 ± 1.1958522385908214e-10 
Exact result: 1.2020569031595951 
Actual error: 5.375033751420233e-12