I constructed a tiny example on which my (naive?) use of the quadgk() fails. The example is a first order polynomial in x and y on the triangular domain with vertices (0,0), (1,0) and (0,1). Below a working and failing MWE. What do I overlook? Thx.
fooworks(x,y) = -8.0*x-8.0*y+4.0
gnuworks(y) = quadgk(x->fooworks(x,y),0,1-y)[1]
fooint = quadgk_count(y->gnuworks(y),0,1)
produces -0.66… as expected output.
# observe change in coefficient (-4 instead of -8) of the term with y
foofails(x,y) = -8.0*x-4.0*y+4.0
gnufails(y) = quadgk(x->foofails(x,y),0,1-y)[1]
fooint = quadgk(y->gnufails(y),0,1)
hangs in performing computations.
You aren’t setting any tolerance in these integrals, so it is using the default relative tolerance of \approx 10^{-8}. This is a problem here because the integral fooint is actually zero, so a relative tolerance |err| ≤ rtol * |int| can never be satisfied (unless the error is exactly zero), so it will keep refining the integrals until it hits the limits of machine precision (which will take a long time in 2d).
Passing an absolute tolerance atol=1e-8 to both quadgk calls causes the calculation to run very quickly.
(You might also look into IteratedIntegration.jl for performing nested quadgk calls more efficiently. Or HCubature.jl via a change of variables to map to a rectangular domain. Or, if you are really just integrating low-degree polynomials and similar very smooth functions, use a tensor-product of fixed-order Gaussian quadrature rules, or SimplexQuad.jl, or similar.)
3 Likes
Sincere thanks for the valuable input. Much appreciated.