How do I obtain the enclosed area of this particular parametric plot? - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-18T14:22:11Z https://mathematica.stackexchange.com/feeds/question/22543 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/22543 18 How do I obtain the enclosed area of this particular parametric plot? Black Milk https://mathematica.stackexchange.com/users/6632 2013-04-02T21:36:33Z 2018-10-12T16:27:58Z <p>I'm trying to find a way to obtain the enclosed area of this particular plot. Can someone show me how?</p> <pre><code>curveplot = ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}, {t, 0, 2 π}, PlotStyle -&gt; {{Thickness[0.01], Darker[Purple]}}, AspectRatio -&gt; Automatic, PlotRange -&gt; All, AxesLabel -&gt; {"x", "y"}] </code></pre> <p><img src="https://i.stack.imgur.com/UWb4w.png" alt="Parametric Plot"></p> https://mathematica.stackexchange.com/questions/22543/-/22545#22545 19 Answer by Daniel Lichtblau for How do I obtain the enclosed area of this particular parametric plot? Daniel Lichtblau https://mathematica.stackexchange.com/users/51 2013-04-02T21:49:11Z 2015-06-30T03:21:46Z <p>You can get the curve in polynomial implicit form as below.</p> <pre><code>poly = GroebnerBasis[{x^2 - ct, y^2 - st, ct^2 + st^2 - 1}, {x, y}, {ct, st}][] (* Out= -1 + x^4 + y^4 *) </code></pre> <p>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).</p> <pre><code>area = Integrate[Boole[poly &lt;= 0], {x, -2, 2}, {y, -2, 2}] (* Out= (2 Gamma[1/4] Gamma[5/4])/Sqrt[π] *) N[area] (* Out= 3.7081493546 *) </code></pre> <p>There are other ways to do this if you cannot find an implicit form, but this seems most direct in this case.</p> <p>--- edit ---</p> <p>If you can just solve separately for <code>x</code> and <code>y</code> in terms of the parameter <code>t</code> 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.</p> <pre><code>reg = Function[{x, y}, If[And @@ {0 &lt;= x &lt;= 1, 0 &lt;= y &lt;= 1, y &lt;= Sqrt[Sin[ArcCos[x^2]]]}, 1, 0]]; approxarea = 4*NIntegrate[reg[x, y], {x, 0, 1}, {y, 0, 1}] (* Out= 3.70814937167 *) </code></pre> <p>One can actually recover the exact area from this by using <code>Integrate</code> instead of <code>NIntegrate</code>. But this seems like a viable approach in situations where the exact value might not be readily computed.</p> <p>--- end edit ---</p> <p>--- edit 2 ---</p> <p>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.</p> <pre><code>segs = Cases[curveplot, _Line, Infinity][[1, 1]]; segs = {Join[segs, N[Table[{j, 1 - j}, {j, 0, 1, 1/100}]]]}; </code></pre> <p>I added an extra level of <code>List</code> due to requirements of some further code. First let's reform a line to take a look at this region.</p> <pre><code>Graphics[Apply[Line, segs]] </code></pre> <p><img src="https://i.stack.imgur.com/oWG6t.png" alt="segment of the curve"></p> <p>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.)</p> <p>To create the in-out function I use code directly from <a href="https://mathematica.stackexchange.com/a/9489">here</a>.</p> <pre><code>nbins = 100; Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} = polyToSegmentList[{segs[]}, nbins];] (* Out= {0.040000, Null} *) len = 100000; pts = RandomReal[1, {len, 2}]; Timing[ inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &amp;, pts];] approxarea = 4.*Length[Cases[inout, True]]/len + 2. (* Out= {2.750000, Null} *) (* Out= 3.7092 *) </code></pre> <p>I would imagine one could do a bit better by integrating just the unit square in the first quadrant with <code>Method -&gt; "QuasiMonteCarlo"</code>, sowing the points via the <code>EvaluationMonitor</code> 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.</p> <p>-- end edit 2 ---</p> https://mathematica.stackexchange.com/questions/22543/-/22552#22552 9 Answer by rm -rf for How do I obtain the enclosed area of this particular parametric plot? rm -rf https://mathematica.stackexchange.com/users/5 2013-04-02T22:45:47Z 2013-04-02T22:45:47Z <p>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.</p> <pre><code>curveplot = ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}, {t, 0, 2 π}, PlotStyle -&gt; Thick, AspectRatio -&gt; 1, Axes -&gt; False, PlotRange -&gt; {-1, 1}] </code></pre> <p><img src="https://i.stack.imgur.com/T4bIL.png" alt=""></p> <p>Next, close the holes in the curve using morphological operations:</p> <pre><code>im = Binarize@curveplot ~Opening~ 7 // DeleteSmallComponents </code></pre> <p><img src="https://i.stack.imgur.com/jh1I8.png" alt=""></p> <p>Depending on the curve, you'll have to tweak the parameters and be careful with <code>DeleteSmallComponents</code>, in case you have disconnected components in your plot (always do a visual check or <code>ImageAdd[curveplot, im]</code> to see if it's correct). Sometimes, this may even not be necessary if your curve is "nice". </p> <p>Finally, use <code>ComponentMeasurements</code> 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})$)</p> <pre><code>4 (1 /. ComponentMeasurements[im, "Area"]) / Times @@ ImageDimensions@im (* 3.66738 *) </code></pre> <p>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 <code>ImageSize</code> in the original plot. Using <code>ImageSize -&gt; 2000</code> gives me 3.70129 as the area (but note that the image processing steps will take longer to compute).</p> https://mathematica.stackexchange.com/questions/22543/-/22553#22553 11 Answer by Michael E2 for How do I obtain the enclosed area of this particular parametric plot? Michael E2 https://mathematica.stackexchange.com/users/4999 2013-04-02T22:48:00Z 2013-04-11T12:29:42Z <p>One can use one of the line integral forms of the <a href="http://mathworld.wolfram.com/Area.html">area</a>, derived from <a href="http://mathworld.wolfram.com/GreensTheorem.html">Green's Theorem</a>: $$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.</p> <pre><code>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 *) </code></pre> <p>One problem with this parametrization are the derivatives of <code>Abs</code> and <code>Sign</code>. 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 <code>1</code> for <code>Sign[x]^2</code>, since <code>x</code> will be <code>0</code> only at few isolated points. Thus the differential is</p> <pre><code>dA = dA /. {Sign'[x_] :&gt; 0, Abs'[x_] :&gt; Sign[x]} /. {Sign[_]^2 :&gt; 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]]]) *) </code></pre> <p><code>NIntegrate</code> returns a small imaginary component</p> <pre><code>NIntegrate[dA, {t, 0, 2 π}] (* 3.70815 - 1.97076*10^-10 I *) </code></pre> <p>Oddly, setting <code>WorkingPrecision</code> reduces the imaginary error, even if it is set to less than <code>MachinePrecision</code> (<code>15.9546</code>)</p> <pre><code>NIntegrate[dA, {t, 0, 2 π}, WorkingPrecision -&gt; 10] (* 3.708149355 + 0.*10^-21 I *) </code></pre> <p><code>Integrate</code> returns an exact answer in this case:</p> <pre><code>Integrate[dA, {t, 0, 2 π}] (* (3 Sqrt π 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 *) </code></pre> https://mathematica.stackexchange.com/questions/22543/-/22587#22587 17 Answer by J. M. will be back soon for How do I obtain the enclosed area of this particular parametric plot? J. M. will be back soon https://mathematica.stackexchange.com/users/50 2013-04-03T11:21:15Z 2015-10-16T16:37:38Z <p>Since it seems to have not been mentioned yet: yet another way to obtain an approximation of the area of your <a href="http://mathworld.wolfram.com/Superellipse.html" rel="nofollow">Lamé curve</a> is to use the <a href="http://en.wikipedia.org/wiki/Shoelace_formula" rel="nofollow">shoelace method</a> for computing the area. Here's a <em>Mathematica</em> demonstration:</p> <pre><code>pts = First[Cases[ ParametricPlot[{Sqrt[Abs[Cos[t]]] Sign[Cos[t]], Sqrt[Abs[Sin[t]]] Sign[Sin[t]]}, {t, 0, 2 π}, Exclusions -&gt; None, Method -&gt; {MaxBend -&gt; 1.}, PlotPoints -&gt; 100] // Normal, Line[l_] :&gt; l, ∞]]; PolygonSignedArea[pts_?MatrixQ] := Total[Det /@ Partition[pts, 2, 1, 1]]/2 PolygonSignedArea[pts] 3.7081447086368127 </code></pre> <p>The value thus obtained is pretty close to the results in the other answers.</p> <hr> <p>Surprisingly, there is in fact an <strong>undocumented</strong> built-in function for computing the area of a polygon:</p> <pre><code>(* $VersionNumber &lt; 10. *) GraphicsMeshMeshInit[]; PolygonArea[pts] 3.708144708636812 (*$VersionNumber &gt;= 10. *) GraphicsPolygonUtils`PolygonArea[pts] 3.708144708636812 </code></pre> <hr> <p>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</p> <p>$$\left|\frac{x}{a}\right|^r+\left|\frac{y}{b}\right|^r=1$$</p> <p>the formula for the area of a Lamé curve (formula 5 <a href="http://mathworld.wolfram.com/Superellipse.html" rel="nofollow">here</a>) is</p> <p>$$A=\frac{4^{1-\tfrac1{r}}ab\sqrt\pi\;\Gamma\left(1+\tfrac1{r}\right)}{\Gamma\left(\tfrac1{r}+\tfrac12\right)}$$</p> <p>In particular, for the OP's specific case, $a=b=1$, and $r=4$. Thus,</p> <pre><code>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 </code></pre> <p>This is the same as the answer Daniel obtained through more general methods.</p> https://mathematica.stackexchange.com/questions/22543/-/22603#22603 6 Answer by BoLe for How do I obtain the enclosed area of this particular parametric plot? BoLe https://mathematica.stackexchange.com/users/6555 2013-04-03T16:22:18Z 2013-04-03T16:22:18Z <p>Integrating <code>InterpolatingFunction</code>.</p> <pre><code>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 </code></pre> https://mathematica.stackexchange.com/questions/22543/-/88484#88484 4 Answer by LCFactorization for How do I obtain the enclosed area of this particular parametric plot? LCFactorization https://mathematica.stackexchange.com/users/9851 2015-07-17T12:24:48Z 2018-10-12T16:27:58Z <p>As a modified version of Michael E2's answer:</p> <p>I tried to rewrite your original curve as below for <span class="math-container">$t\in \left[0,\tfrac{\pi}{2}\right]$</span>, in order to make sure the derivatives of the parametric form can be obtained easily by avoiding <code>Abs</code> or <code>Sign</code>:</p> <pre><code>ncurve = {(Cos[t]^2 )^(1/4), (Sin[t]^2)^(1/4)} </code></pre> <p>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:</p> <pre><code>area = 4*Integrate[ncurve[] D[ncurve[], t], {t, 0, Pi/2}] </code></pre> <p>which gives:</p> <blockquote> <p><span class="math-container">$$\dfrac{8\; \Gamma^2 \left(\tfrac{5}{4}\right)}{\sqrt{\pi }}$$</span></p> </blockquote> <p>and the numerical result is therefore:</p> <pre><code>area // N[#, 20] &amp; </code></pre> <blockquote> <p>3.7081493546027438369</p> </blockquote>