I am trying to do two-dimensional numerical integration. For this my first attempt was something like the following, using DoubleExponentialFormulas.jl. This example is significantly reduced and no longer does anything mathematically useful:
using DoubleExponentialFormulas #using QuadGK function f4(c, t, S, α) t^(c - 1) * S - (c - 1) * α end function f3(j, t, S, α) S * prod((f4(c, t, S, α) / c) for c = 0:j) end function f2(n, a, b, r, ξ, η) sum(f3(j, (r - b) / a, r, b / a) for j = 0:(n+2)) end function inner_integral(n, a, b, δ, ξ) f(η) = f2(n, a, b, a, ξ, η) I, E = quadde(f, 0, 2π, rtol = δ) I end function outer_integral(n, a, b, δ) f(ξ) = inner_integral(n, a, b, δ, ξ) * cos(b + ξ) I, E = quadde(f, -π / 2 - b, 0, π / 2 - b, rtol = δ) I end function sum_of_integrals(N, a, b, δ) sum(outer_integral(n, a, b, δ) for n = 0:N) end @time println(sum_of_integrals(3, 100, 0.2, 10^-6))
However, I cannot get this code to compile within a reasonable timeframe. I left a variant of this running for 3.5 hours before I ran out of system memory. I have verified that this time is actually spent in (some early step of) compilation: inserting a statement that should lead to a compilation error, e.g. a call to a non-existent function, does not cause an error during all this time.
Using QuadGK.jl instead (by switching out the
using statements and replacing the
quadde calls with
quadgk) makes the compilation work but it nevertheless takes around 90 seconds, which still seems extremely slow to me. (Regardless, for mathematical reasons I have to use double-exponential quadrature for integration, so I cannot use
quadgk) There are some other things that speed up the compilation some more, most importantly inlining the
f4 functions as loops within
f2 or making the integration intervals smaller, but it never reaches a usable point with
quadde for me.
I am running Julia 1.6.1 on Arch Linux; my CPU is an i5-8350U.
I am aware that nested 1D integration is not the best way to do multidimensional integration, but at the same time I feel the compiler should not behave in this way. In a sense, I am not really looking for a solution to make the compilation run faster (the real solution would probably be “use an actual multidimensional integrator”), but more asking whether this is expected behaviour from the compiler, and if not, whether it would make sense to investigate this further on my part so that the compiler may be improved in the future.