Làm thế nào để thực hiện hai tích hợp số biến trong Julia?

Question

Tôi có thể thực hiện tích hợp số biến duy nhất trong Julia bằng cách sử dụng quadgk. Một số ví dụ đơn giản:

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 

Các documentation for quadgk dường như không bao hàm một sự hỗ trợ cho hội nhập đa chiều, và đủ chắc chắn tôi nhận được một lỗi nếu tôi cố gắng để lạm dụng nó cho một 2D không thể thiếu:

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

Tài liệu này cho thấy có một số gói bên ngoài để tích hợp, nhưng không đặt tên cho chúng. Tôi đoán rằng một trong những gói như vậy có thể làm tích phân hai chiều. Gói tốt nhất cho công việc này là gì?

Answer 1

Tôi nghĩ rằng bạn sẽ muốn kiểm tra các gói Cubature:

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

Có thể cho rằng, quadgk nên chỉ đơn giản là được lấy ra từ thư viện chuẩn vì nó hạn chế và chỉ lừa người vào không tìm kiếm một gói để làm tích hợp.

Answer 2

Ngoài Cubature.jl, còn có một gói Julia khác cho phép bạn tính toán tích phân số nhiều chiều: Cuba.jl (https://github.com/giordano/Cuba.jl). Bạn có thể cài đặt nó sử dụng quản lý gói:

Pkg.add("Cuba") 

Các tài liệu đầy đủ của gói có sẵn tại https://cubajl.readthedocs.org (cũng trong PDF version)

Disclaimer: Tôi là tác giả của gói.

---

Cuba.jl chỉ đơn giản là một wrapper Julia xung quanh Cuba Library, bởi Thomas Hahn, và cung cấp bốn thuật toán độc lập để tính tích phân: Vegas, Suave, Divonne, Cuhre.

Các thể thiếu của cos (x) trong lĩnh vực [0, 1] có thể được tính với một trong các lệnh sau:

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, tích phân

nơi Ω = [0, 1] ³ và

có thể được tính với Julia kịch bản sau đây:

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 

mang đến cho đầu ra sau đây

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