19
$\begingroup$

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

Improve this question
$\endgroup$
0

6 Answers 6

Reset to default
19
$\begingroup$

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[π] *)

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]]

segment of the curve

Now we create an in-out function, generate a bunch of random points in the 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 ---

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks, I wanted to see if this checked out with a monte carlo estimate. $\endgroup$
    – Black Milk
    Apr 2, 2013 at 21:55
18
$\begingroup$

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:

(* $VersionNumber < 10. *)
Graphics`Mesh`MeshInit[];
PolygonArea[pts]    
   3.708144708636812

(* $VersionNumber >= 10. *)
Graphics`PolygonUtils`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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Never heard of that method before. Nice. (Yes, I upvoted..) $\endgroup$
    – Daniel Lichtblau
    Apr 3, 2013 at 15:42
  • $\begingroup$ I've taught the "shoelace" method but never heard of the name. (+1) $\endgroup$
    – Michael E2
    Apr 3, 2013 at 20:32
  • $\begingroup$ @J.M. Thank you for the additional insight. $\endgroup$
    – Black Milk
    Apr 3, 2013 at 20:55
  • $\begingroup$ @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 ... $\endgroup$
    – J. M.'s eventual burnout ♦
    Apr 4, 2013 at 2:37
  • $\begingroup$ You could just do PolygonArea@pts... no need to wrap it in Polygon $\endgroup$
    – rm -rf ♦
    Apr 25, 2013 at 21:19
12
$\begingroup$

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 *)
Improve this answer
$\endgroup$
9
$\begingroup$

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).

Improve this answer
$\endgroup$
1
  • 4
    $\begingroup$ Exclusions -> None in ParametricPlot[] should close up the artifact holes nicely. $\endgroup$
    – J. M.'s eventual burnout ♦
    Apr 3, 2013 at 15:27
6
$\begingroup$

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
Improve this answer
$\endgroup$
4
$\begingroup$

As a modified version of Michael E2's answer:

I tried to rewrite your original curve as below for $t\in \left[0,\tfrac{\pi}{2}\right]$, in order to make sure the derivatives of the parametric form can be obtained easily by avoiding Abs or Sign:

ncurve = {(Cos[t]^2 )^(1/4), (Sin[t]^2)^(1/4)}

Then the result of the closed area can be obtained by applying Green's theorem to the curve and by using the symmetry of the curve:

area = 4*Integrate[ncurve[[1]] D[ncurve[[2]], t], {t, 0, Pi/2}]

which gives:

$$\dfrac{8\; \Gamma^2 \left(\tfrac{5}{4}\right)}{\sqrt{\pi }}$$

and the numerical result is therefore:

area // N[#, 20] &

3.7081493546027438369

Improve this answer
$\endgroup$

Your Answer

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

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