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 f3
and 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.