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


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?


2 Answers 2


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.


This paper seems to be a good reference:

Average Case Complexity of Weighted Integration and Approximation over $\mathbb{R}^d$ with isotropic weight, by Plaskota, Ritter, Wasilkowski.

Of course, there are a million different models, of which the above paper addresses only one...


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.