How to make two numerical integration variables in Julia?

Question

I can do a single variable numerical 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 does not seem to imply support for multidimensional integration, and of course I get an error if I try to abuse it for a two-dimensional integral:

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

The documentation does have some external integration packages, but they are not called. I assume that one such package can perform two-dimensional integrals. What is the best such package for this task?

Answer 1

I think you will want to check out the Cubature package:

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

Perhaps quadgk should simply be removed from the standard library, as it is limited and simply misleads people without looking for an integration package.

Answer 2

In addition to Cubature.jl there is another Julia package that allows you to calculate multidimensional numerical integrals: Cuba.jl ( https://github.com/giordano/Cuba.jl ). It can be installed using the package manager:

 Pkg.add("Cuba") 

Full package documentation is available at https://cubajl.readthedocs.org (also in PDF version )

Disclaimer: I am the author of the package.

---

Cuba.jl is just a wrapper for Julia around the Cuban Library from Thomas Hahn and provides four independent algorithms for calculating the integrals: Vegas, Suave, Divonne, Kure.

The integral cos (x) in the region [0, 1] can be calculated using 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 calculated using the following version of the 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 next exit

 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