9

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.

Shuhao Cao
  • 2,552
  • 17
  • 30
  • 2
    If the origin is within the integration domain, may I suggest breaking up your integral and then transforming each one to spherical coordinates? – J. M. Dec 16 '11 at 04: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. –  Dec 16 '11 at 19:36

1 Answers1

8

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./3
quad(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.

J. M.
  • 3,155
  • 28
  • 37
David Ketcheson
  • 16,522
  • 4
  • 54
  • 105