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.
What numerical quadrature to choose to integrate a function with singularities?
numerics;finite element;quadrature
You should be able to get accurate results with mpmath, a Python module for arbitrary-precision floating-point computations.  There are examples of integration with singularities in the documentation.  You'll want to explicitly tell it to break up the interval:from mpmath import *f = lambda x,y,z: 1./(x**2+y**2+z**2)**1./3quad(f,[-1,0,1],[-1,0,1],[-1,0,1])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.