Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Given a two-variable function f[x,y], how to get the contour plots corresponding to several values of the integral of the original function. For example, I would like to plot the contours $C_j$ of the domains $S_j$, in the $xy$ plane, such that

$$ \int_{S_j} dxdy f(x,y) = \alpha_j $$

for $\alpha_j = \alpha_{j-1} + \epsilon$, with a fixed step $\epsilon$.

Additionally, I would like to be able to choose a value of $\alpha$ and plot the contour of the corresponding domain.

My difficulty is the following. I don't know how to make a condition inside ContourPlot such that it tells that the contours should correspond to the values of the integrals.

A possible application would be to determine the region corresponding to 68.3% of the events of a two-variable Gaussian distribution. This is called a 1-sigma region.

share|improve this question
Could you give an example? Any restrictions on f or the contours? –  ssch Mar 4 '13 at 23:01
It doesn't seem to me that you put much work into this before you asked for help here. We're here to help you get out of difficulties, not to do your work for you, –  m_goldberg Mar 4 '13 at 23:11
I am facing this kind of problem in my codes and I've tried to state it in a very general way, so that a possible answer could be useful to others. –  fcpenha Mar 4 '13 at 23:17
It looks like you may have over-generalized this question into something meaningless. These "contours" are not well defined. First, surely you mean to refer to regions, because these are not contour integrals; they are 2D integrals. Also, with few exceptions, for any given $\alpha_j$ there will be a large set of regions, having differing boundaries $C_j$, which yield the same value of the integral. Do you intend that $f$ be constant on each $C_j$? –  whuber Mar 4 '13 at 23:18
I have edited my question after these comments. I hope it is better, now. –  fcpenha Mar 5 '13 at 0:00
show 1 more comment

1 Answer

Here is a difficult running example:

f[x_, y_] := Abs[x + y]^(1/2) Cos[x^2 + y^2/2]^2;
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, Boxed -> False]


I interpret the question in this way: given a value $z$ in the range of $f$, we need to be able to compute the integral of $f$ over the set of $R_z=\{(x,y)\ |\ f(x,y) \le z\}$. Here is a general way:

i[f_, z_, xlim_: {-Infinity, Infinity}, ylim_: {-Infinity, Infinity}] := 
  Block[{x, y, g = Function[{x, y}, u = f[x, y]; u Boole[u <= z]]}, 
  Chop[NIntegrate[g[x, y], Evaluate@Prepend[xlim, x], Evaluate@Prepend[ylim, y]]]];

Its arguments are $f$, $z$, and the endpoints of a rectangle containing the domain of $f$ (defaulting to all of $\mathbb{R}^2$. For example,

i[f, 0.5, {-2, 2}, {-2, 2}]

returns $1.69395$, which is the numerical estimate of the integral $\int_{R_{0.5}} f(x,y) dx dy$.

The question asks how to invert $i$. To this end, let's compute the values of $i(f,z)$ for a carefully chosen range of values of $z$ matching the characteristics of $f$ (as gleaned by inspecting the preceding plot):

x = Range[2/10, 13/10, 1/10]^1.5;
data = ParallelMap[i[f, #, {-2, 2}, {-2, 2}] &, x]

Because $f$ is a somewhat difficult function to integrate, this takes some time and raises a lot of warnings. You can improve the accuracy according to needs and time available. But let's continue. Here is what we just computed:

Show[ListPlot[{x, data}\[Transpose], Joined -> True, PlotStyle -> Thick], 
  ListPlot[{x, data}\[Transpose], PlotStyle -> PointSize[0.015]]]

Plot 2

That's smooth enough for this illustration; here's the interpolating function of its inverse:

h = Interpolation[{data, x}\[Transpose]]

With it we can quickly produce any of the desired contours. Here, we show them for $z=1,2,\ldots, 7$:

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, Contours -> h /@ Range[1, 7]]

Plot 3

The contours themselves are labeled (as tooltips, not visible in this figure) with the $h(z_i)$, giving the height cutoff for $f$ corresponding to each contour.

Obviously there are some errors near the extreme values of $z$, related to the finite accuracy with which the integrals have been computed. Recompute them with higher accuracy if desired. For $f$ with sufficiently simple expressions, NIntegrate can be replaced by Integrate in i to achieve perfect accuracy and the interpolation in h can even be replaced by a call to Solve.

share|improve this answer
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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