How to do two variable numeric integration in Julia?

Question

I can do single variable numeric integration in Julia using quadgk. Some simple examples:

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

The documentation for quadgk doesn't seem to imply an support for multidimensional integration, and sure enough I get an error if I attempt to misuse it for a 2D integral:

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

The documentation does suggest there are some external packages for integration, but doesn't name them. I'm guessing that one such package can do two dimensional integrals. What's the best such package for this task?

Answer 1

I think you'll want to check out the Cubature package:

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

Arguably, quadgk should simply be removed from the standard library because it's limited and just misleads people into not looking for a package to do integration.

Answer 2

In addition to Cubature.jl, there is another Julia package that allows you to compute multidimensional numerical integrals: Cuba.jl (https://github.com/giordano/Cuba.jl). You can install it using package manager:

Pkg.add("Cuba")

The complete documentation of the package is available at https://cubajl.readthedocs.org (also in PDF version)

Disclaimer: I'm the author of the package.

---

Cuba.jl is simply a Julia wrapper around Cuba Library, by Thomas Hahn, and provides four independent algorithms to calculate integrals: Vegas, Suave, Divonne, Cuhre.

The integral of cos(x) in the domain [0, 1] can be computed with one of the following commands:

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, the integral

where Ω = [0, 1]³ and

can be computed with the following Julia script:

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

which gives the following output

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