# What numerical quadrature to choose to integrate a function with singularities?

For example, I would like to numerically compute the $L^2$-norm of $\displaystyle u = \frac{1}{(x^2+y^2+z^2)^{1/3}}$ in some domain that includes zero, I tried Gauss quadrature and it fails, it is kinda far from the real $L^2$-norm on the unit ball using spherical coordinates to integrate, is there any good way to do this? This problem is often seen in the finite element computing toy problems for domains with re-entrant corners. Thanks.

• If the origin is within the integration domain, may I suggest breaking up your integral and then transforming each one to spherical coordinates? Dec 16 '11 at 4:50
• I agree with JM -- if you know the location and structure of the singularities beforehand, you're better off using that structural information in writing the calls to your quadrature routines intelligently vs. feeding it to a numerical package and hoping that (a) it finds the singularities and (b) does the right thing with them.
– user389
Dec 16 '11 at 19:36

from mpmath import *

You may need to increase the precision (e.g. mp.dps=30) and it will likely be slow, but should be quite accurate.
You could also try nesting calls to MATLAB's quadgk(), which uses adaptive Gauss-Kronrod quadrature in 1D.