I have an algorithm for solving a certain problem that requires that I compute two-dimensional integrals as a subroutine, and I'd like to make some kind of statement about its running time. Suppose $S$ is a simply connected region in the plane whose boundary is composed of $n$ line segments and $m$ arcs of some type (say, arcs of a circle, ellipse, or parabola). I'd like to numerically approximate the integral of some function $f(x)$ over $S$. What's the "computational complexity", in big-O notation, of computing this within a given error $\epsilon$? Clearly this depends on lots of factors, such as the curvature of the arcs and the derivatives of the function $f(x)$ but I haven't had much success in pinning down anything concrete. By way of comparison, the wikipedia page for Simpson's rule
http://en.wikipedia.org/wiki/Simpson%27s_rule#Error
says that the error committed by the composite Simpson's rule is bounded (in absolute value) by
$\frac{h^4}{180} (b-a) M$
where $h$ is the step size and $M = \max_{\xi\in[a,b]}|f^{(4)}(\xi)|$, and therefore the complexity of computing such an integral within precision $\epsilon$ is $M^{1/4}(b-a)^{5/4}\epsilon^{-1/4}$. Can I make any similar statement for a two-dimensional integral over my region $S$? What if I were able to triangulate $S$ into $N$ small pieces?
Best Answer
Since you mentioned triangulation, Quarteroni, Sacco, and Saleri's book gives the following formula in the section on multi-dimensional integration: if your domain of integration is a convex polygon $\Omega$ with a triangulation $\mathcal{T}_h$ consisting of $N_T$ triangles, where $h$ is the maximum edge length in $\mathcal{T}_h$, and your quadrature rule has degree of exactness equal to $n$ with all nonnegative weights, then there exists a positive constant $K_n$, independent of $h$, such that the error $E$ is bounded by
$|E| \leq K_n h^{n+1} |\Omega| M_{n+1}$
where $M_{n+1}$ is the maximum value of the modules of the derivatives of order $n+1$ and $|\Omega|$ is the area of $\Omega$. The composite midpoint formula and the composite trapezoidal formula have nonnegative weights and degree of exactness $n=1$ and so we have
$|E| \leq K_1 h^2 |\Omega| M_2$
for some constant $K_1$. It follows that, if your desired error is $\epsilon$, then the maximum length of any edge of your triangulation $h$ must be at most $\epsilon^{1/2}(|\Omega| M_2 K_1)^{-1/2}$, i.e. $O(\epsilon^{1/2}(|\Omega| M_2)^{-1/2})$. If you break $\Omega$ into $N_T$ triangles as I mentioned before, the edge length will be $O(\sqrt{|\Omega|/N_T})$, so you'll need to break $\Omega$ into $O(\frac{|\Omega|^2 M_2}{\epsilon})$ triangles. So, quadratic in the area, linear in the maximum value of the second derivative, and inversely proportional to $\epsilon$. It's not clear to me where convexity enters the picture, so perhaps this will adapt to the general case you've described. The authors state that a proof can be found in Isaacson E. and Keller H. (1966), Analysis of Numerical Methods. Wiley, New York. (This is not my area of expertise so I welcome any corrections)
UPDATE: I checked the Isaacson book and the relevant section, "7.4 Composite Formulae for Multiple Integrals" doesn't make any mention of convexity, so you may be safe in using the $O(\frac{|\Omega|^2 M_2}{\epsilon})$ that I gave before.