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.