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.

I'm trying to find a way to obtain the enclosed area of this particular plot. Can someone show me how?

curveplot = ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], 
Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}, {t, 0, 2 π}, 
PlotStyle -> {{Thickness[0.01`], Darker[Purple]}}, 
AspectRatio -> Automatic, PlotRange -> All, AxesLabel -> {"x", "y"}]

Parametric Plot

share|improve this question

5 Answers 5

up vote 11 down vote accepted

You can get the curve in polynomial implicit form as below.

poly = 
 GroebnerBasis[{x^2 - ct, y^2 - st, ct^2 + st^2 - 1}, {x, y}, {ct, 
    st}][[1]]

(* Out[290]= -1 + x^4 + y^4 *)

To get the area, integrate the characteristic function for the interior of the region. That that's where the polynomial is nonpositive (just notice that it is negative at the origin, say).

area = Integrate[Boole[poly <= 0], {x, -2, 2}, {y, -2, 2}]

(* Out[292]= (2 Gamma[1/4] Gamma[5/4])/Sqrt[\[Pi]] *)

N[area]

(* Out[293]= 3.7081493546 *)

There are other ways to do this if you cannot find an implicit form, but this seems most direct in this case.

--- edit ---

If you can just solve separately for x and y in terms of the parameter t then you can set up a region function. I do this below for the positive quadrant, and take advantage of symmetry to get the full area in approximate form.

reg = Function[{x, y}, 
   If[And @@ {0 <= x <= 1, 0 <= y <= 1, y <= Sqrt[Sin[ArcCos[x^2]]]}, 
    1, 0]];

approxarea = 4*NIntegrate[reg[x, y], {x, 0, 1}, {y, 0, 1}]

(* Out[321]= 3.70814937167 *)

One can actually recover the exact area from this by using Integrate instead of NIntegrate. But this seems like a viable approach in situations where the exact value might not be readily computed.

--- end edit ---

--- edit 2 ---

Here is a Monte Carlo method that does not rely on solving for anything. We extract the line segments, augment with a diagonal, and do some magic.

segs = Cases[curveplot, _Line, Infinity][[1, 1]];
segs = {Join[segs, N[Table[{j, 1 - j}, {j, 0, 1, 1/100}]]]};

I added an extra level of List due to requirements of some further code. First let's reform a line to take a look at this region.

Graphics[Apply[Line, segs]]

enter image description here

Now we create an in-out function, generate a bunch of random points in teh unit square of the first quadrant, take a Monte-Carlo approximation of this area. Then multiply by 4 and add 2. Why? because that's what one always does-- it's like selecting "c" when we don't know the multiple choice answer. (Okay, we multiply by 4 to account for all quadrants, and add 2 because we have in effect excised a square of side length sqrt(2) from the full region.)

To create the in-out function I use code directly from here.

nbins = 100;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} = 
   polyToSegmentList[{segs[[1]]}, nbins];]

(* Out[414]= {0.040000, Null} *)

len = 100000;
pts = RandomReal[1, {len, 2}];
Timing[
 inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &,
     pts];]
approxarea = 4.*Length[Cases[inout, True]]/len + 2.

(* Out[419]= {2.750000, Null} *)

(* Out[420]= 3.7092 *)

I would imagine one could do a bit better by integrating just the unit square in the first quadrant with method->"QuasiMonteCarlo", sowing the points via the EvaluationMonitor option, and using those instead of the random set above. This will give a low-discrepancy sequence. Or generate such a set directly; bit offhand I don't know how to do that.

-- end edit 2 ---

share|improve this answer
    
Thanks, I wanted to see if this checked out with a monte carlo estimate. –  Black Milk Apr 2 '13 at 21:55

Since it seems to have not been mentioned yet: yet another way to obtain an approximation of the area of your Lamé curve is to use the shoelace method for computing the area. Here's a Mathematica demonstration:

pts = First[Cases[
          ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]},
                         {t, 0, 2 π}, Exclusions -> None, Method -> {MaxBend -> 1.},
                         PlotPoints -> 100] // Normal, Line[l_] :> l, ∞]];

PolygonSignedArea[pts_?MatrixQ] := Total[Det /@ Partition[pts, 2, 1, 1]]/2

PolygonSignedArea[pts]
   3.7081447086368127

The value thus obtained is pretty close to the results in the other answers.


Surprisingly, there is in fact an undocumented built-in function for computing the area of a polygon:

Graphics`Mesh`MeshInit[];
PolygonArea[pts]
   3.708144708636812

But, what you really should know is that Lamé curves have been well studied, and there is in fact a closed form expression for the area of a Lamé curve. Given the Cartesian equation

$$\left|\frac{x}{a}\right|^r+\left|\frac{y}{b}\right|^r=1$$

the formula for the area of a Lamé curve (formula 5 here) is

$$A=\frac{4^{1-\tfrac1{r}}ab\sqrt\pi\;\Gamma\left(1+\tfrac1{r}\right)}{\Gamma\left(\tfrac1{r}+\tfrac12\right)}$$

In particular, for the OP's specific case, $a=b=1$, and $r=4$. Thus,

With[{a = 1, b = 1, r = 4}, 
 N[(4^(1 - 1/r) a b Sqrt[π] Gamma[1 + 1/r])/Gamma[1/r + 1/2], 20]]
   3.7081493546027438369

This is the same as the answer Daniel obtained through more general methods.

share|improve this answer
    
Never heard of that method before. Nice. (Yes, I upvoted..) –  Daniel Lichtblau Apr 3 '13 at 15:42
    
I've taught the "shoelace" method but never heard of the name. (+1) –  Michael E2 Apr 3 '13 at 20:32
    
@J.M. Thank you for the additional insight. –  Black Milk Apr 3 '13 at 20:55
    
@Michael, I must confess that I too have been using that method for a while now, but only learned it had a name on math.SE ... –  J. M. Apr 4 '13 at 2:37
    
You could just do PolygonArea@pts... no need to wrap it in Polygon –  rm -rf Apr 25 '13 at 21:19

One can use one of the line integral forms of the area, derived from Green's Theorem: $$A = \frac12 \int_C x \; dy - y \; dx = \int_C x \; dy = - \int_C y \; dx$$ The first one is symmetric, which sometimes is an advantage.

c[t_] := {Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}

dA = 1/2 c'[t].Cross[c[t]]
(* complicated output *)

One problem with this parametrization are the derivatives of Abs and Sign. They are discontinuous at isolated points, and as far as the integral is concerned, it does not matter what value we assign them at the discontinuities. So we can simplify matters by substituting for them. We can also substitute 1 for Sign[x]^2, since x will be 0 only at few isolated points. Thus the differential is

dA = dA /. {Sign'[x_] :> 0, Abs'[x_] :> Sign[x]} /. {Sign[_]^2 :> 1} // Simplify
(* (Abs[Cos[t]] Cos[t] Sign[Cos[t]] + Abs[Sin[t]] Sign[Sin[t]] Sin[t]) /
   (4 Sqrt[Abs[Cos[t]]] Sqrt[Abs[Sin[t]]]) *)

NIntegrate returns a small imaginary component

NIntegrate[dA, {t, 0, 2 π}]
(* 3.70815 - 1.97076*10^-10 I *)

Oddly, setting WorkingPrecision reduces the imaginary error, even if it is set to less than MachinePrecision (15.9546)

NIntegrate[dA, {t, 0, 2 π}, WorkingPrecision -> 10]
(* 3.708149355 + 0.*10^-21 I *)

Integrate returns an exact answer in this case:

Integrate[dA, {t, 0, 2 π}]
(* (3 Sqrt[2] π Gamma[5/4] + 4 Gamma[3/4] Gamma[5/4]^2) /
   (2 Sqrt[π] Gamma[3/4]) *)

FullSimplify @ %
(* (Sqrt[π/2] Gamma[1/4])/Gamma[3/4] *)

N[%, 20]
(* 3.7081493546027438369 *)
share|improve this answer

In cases where it is not possible to get a closed form for the curve, you can use the image processing functions to get an approximation for the enclosed area. First, some small changes to the plot — get rid of the axes, labels and everything else that isn't needed, and set the aspect ratio to 1.

curveplot = ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}, 
    {t, 0, 2 π}, PlotStyle -> Thick, AspectRatio -> 1, Axes -> False, PlotRange -> {-1, 1}]

Next, close the holes in the curve using morphological operations:

im = Binarize@curveplot ~Opening~ 7 // DeleteSmallComponents

Depending on the curve, you'll have to tweak the parameters and be careful with DeleteSmallComponents, in case you have disconnected components in your plot (always do a visual check or ImageAdd[curveplot, im] to see if it's correct). Sometimes, this may even not be necessary if your curve is "nice".

Finally, use ComponentMeasurements to get the area of the enclosed space. The result is in sq. pixels, so I normalize by the total area (in sq. pixels) and multiply by the area of the plot range rectangle in the original plot (i.e., $(x_{max}-x_{min})(y_{max}-y_{min})$)

4 (1 /. ComponentMeasurements[im, "Area"]) / Times @@ ImageDimensions@im
(* 3.66738 *)

which looks about right, since your curve fits in a square of side 2 and is close to Daniel's answer of 3.708. You can get a closer approximation if you increase the ImageSize in the original plot. Using ImageSize -> 2000 gives me 3.70129 as the area (but note that the image processing steps will take longer to compute).

share|improve this answer
2  
Exclusions -> None in ParametricPlot[] should close up the artifact holes nicely. –  J. M. Apr 3 '13 at 15:27

Integrating InterpolatingFunction.

p = Table[{
   Sign[Cos@t] Sqrt[Abs@Cos@t],
   Sign[Sin@t] Sqrt[Abs@Sin@t]}, {t, 0, Pi, Pi/100.}];
f = Interpolation[p];
2*NIntegrate[f@t, {t, -1, 1}]
3.70771
share|improve this answer

Your Answer

 
discard

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.