Finding Volume Of A Segment - Mathematica Stack Exchange

Sorry, we no longer support your browser Please upgrade to Microsoft Edge, Google Chrome, or Firefox. Learn more about our browser support.

1. 1. Home
   2. Questions
   3. Tags
   4. Users
   5. Unanswered
2. Teams

   Ask questions, find answers and collaborate at work with Stack Overflow for Teams.

   Try Teams for free Explore Teams
3. Teams

4. Ask questions, find answers and collaborate at work with Stack Overflow for Teams. Explore Teams

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

Learn more about Teams Finding volume of a segment Ask Question Asked 10 years, 1 month ago Modified 8 years, 3 months ago Viewed 940 times 8 $\begingroup$

I'm still pretty new to Mathematica, so I would like to seek advice regarding a geometrical problem.

I am currently trying to define that as an extra condition in the Mathematica code below.

reg = ImplicitRegion[x^2/a^2 + y^2/b^2 + z^2/c^2 <= 1 {z, y, x}]; Volume[reg, Assumptions -> a > 0 && b > 0 && c > 0]

Any one has any idea how to incorporate it into the extra conditions in defining the implicit region?

• calculus-and-analysis
• geometry
• regions

Share Improve this question Follow edited Aug 12, 2016 at 3:40 J. M.'s missing motivation♦ 125k1111 gold badges407407 silver badges581581 bronze badges asked Oct 10, 2014 at 4:42 CorseCorse 45922 silver badges77 bronze badges $\endgroup$ 2

• $\begingroup$ What is the equation of the plane ....? $\endgroup$ – Dr. belisarius Commented Oct 10, 2014 at 4:53
• $\begingroup$ In this case, the equation of the plane (XZ) is y = 0 right? I would need to rotate this plane about the point X and use the rotated plane as a condition in evaluating the volume. Not sure how do I do that though. $\endgroup$ – Corse Commented Oct 10, 2014 at 5:01

Add a comment |

3 Answers 3

Sorted by: Reset to default Highest score (default) Date modified (newest first) Date created (oldest first) 13 $\begingroup$

Let you have a vector ${\bf p}$, which is perpendicular to the plane and an ellipsoid with axes $(a,b,c)$. The illustration (2D for simplicity):

Mathematica can calculate the numeric value of the clipped volume easily

Nvolume[p_, abc_] := Volume[RegionIntersection[ ImplicitRegion[{x, y, z}.N[p] > 0, {x, y, z}], Ellipsoid[N[abc] {1, 0, 0}, N[abc]]]] p = RandomReal[{-1, 1}, 3]; abc = RandomReal[{1, 2}, 3]; Nvolume[p, abc] (* 16.2584 *)

Mathematica cannot derive the general formula, but it isn't difficult to derive manually. Let us introduce new coordinates

$$ x' = x/a, \quad y' = y/b, \quad z' = z/c. $$

In these coordinates the ellipsoid becomes the unit ball

The Jacobian of this transformation is $J=abc$. In the new coordinates the normalized perpendicular vector is

$$ {\bf n} = \frac{(ap_x,bp_y,cp_z)}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}. $$

Now it is simple to integrate the volume along the axis $\xi$ because the cross section is a circle

$$ V=abc\int_{-n_x}^1\pi(1-\xi^2)d\xi = \pi abc \left(\frac{2}{3} + n_x - \frac{n_x^3}{3}\right) $$

volume[p_, abc_] := π Times @@ abc (2/3 + # - #^3/3) &@@ Normalize[abc p] volume[p, abc] (* 16.2584 *)

The result is the same.

Update: OP asks also about the area of the intersection. It is also an interesting question.

Mathematica region functionality is very powerful for numerical computations:

Narea[p_, abc_] := Area[RegionIntersection[ImplicitRegion[{x, y, z}.N[p] == 0, {x, y, z}], Ellipsoid[N[abc] {1, 0, 0}, N[abc]]]] Narea[p, abc] (* 6.20243 *)

The analytic formula can be derived using the Dirac $\delta$-function \begin{multline} A = \int_\text{ellipse} \delta \left({\bf r}\cdot\frac{{\bf p}}{p}\right)d{\bf r} = abcp \int_\text{unit ball} \delta \left(x'ap_x+y'bp_y+z'cp_z\right)d{\bf r}' = \\ \frac{abcp}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}\int_\text{unit ball} \delta \left({\bf r}'\cdot{\bf n}\right)d{\bf r}'. \end{multline} It is the cross section of the unit ball. Hence \begin{equation} A = \frac{\pi abcp (1-n_x^2)}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}. \end{equation}

area[p_, abc_] := π Times @@ abc (1 - #^2) & @@ Normalize[abc p] Norm[p]/Norm[abc p]; area[p, abc] (* 6.20243 *) Share Improve this answer Follow edited Oct 21, 2014 at 11:48 answered Oct 10, 2014 at 20:00 ybeltukovybeltukov 43.9k55 gold badges109109 silver badges215215 bronze badges $\endgroup$ 15

• $\begingroup$ this is beautiful method and exposition $\endgroup$ – ubpdqn Commented Oct 10, 2014 at 23:33
• $\begingroup$ Hi there @ybeltukov, thats a concise analytical method nicely complementing ubpdqn's proposed implementation! I understand that a plane takes {x,y,z}.N[p] = 0, but what does it mean when {x, y, z}.N[p] > 0 ? Also, how could I get a rotation angle of the plane w.r.t horizontal using the normal vector p? $\endgroup$ – Corse Commented Oct 11, 2014 at 6:32
• $\begingroup$ @Corse, {x, y, z}.N[p] > 0 means the half of the space on the one side of the plane. As in ubpdqn's you can use p = {Sin[a], 0, Cos[a]}. $\endgroup$ – ybeltukov Commented Oct 11, 2014 at 12:51
• $\begingroup$ @ybeltukov thanks for the clarification, i realise that the value of 'a' has to be in terms of a negative angle to get the volume under the rotated plane. Sorry to bother, could there be a simple modification to the code to get the corresponding area of the intersecting plane through the ellipsoid? $\endgroup$ – Corse Commented Oct 12, 2014 at 2:29
• $\begingroup$ @Corse, you can change the sign of p (e.g. p = -{Sin[a], 0, Cos[a]}) or change > 0 to < 0 in Nvolume and signs of # and #^3 in volume $\endgroup$ – ybeltukov Commented Oct 12, 2014 at 2:34

| Show 10 more comments 8 $\begingroup$

Visualization of the @ybeltukov's answer

This is not an answer, just visualization of ybeltukov's answer.

Step 1 direction3D for direction handling.

opt1 = Sequence[Orange, Thick, Arrowheads[.15]]; opt2 = Sequence[PlotRange -> 1, Ticks -> None, Boxed -> False, ViewPoint -> {1, 1, 4}, ViewVertical -> {0, 1, 0}, PlotRegion -> {{-.2, 1.05}, {-.2, 1.05}}, ImageSize -> 60]; direction3D[Dynamic[d_]] := DynamicModule[ {$p = {0, π/2}, dd}, LocatorPane[Dynamic[$p, ($p = #; dd @@ #) &], Graphics3D[ {{opt1, Dynamic@Arrow[Tube[{{0, 0, 0}, d}, .04]]}, {Dynamic@Cylinder[{{0, 0, 0}, d/100}]}}, opt2 ], {{-π, π}, {π, 0}},Appearance -> None, ImageMargins -> 0 ], Initialization :> (d = {1, 0, 0}; dd[t_, s_] := (d = {Sin[t] Sin[s], Cos[s], Cos[t] Sin[s]})) ]

Step 2 ybeltukov's result.

volume[p_, abc_] := π Times @@ abc (2/3 + # - #^3/3) & @@Normalize[abc p]

step3 visualization

opts = Sequence[Boxed -> False, Axes -> True, ViewPoint -> {1, 1, 4}, ViewVertical -> {0, 1, 0}, AxesOrigin -> {0, 0, 0}, PlotRange -> {{-3, 3}, {-2, 2}, {-2, 2}}]; Manipulate[ elp = Graphics3D[{Opacity[.7], Orange, Ellipsoid[{0, 0, 0}, {a, b, c}], Opacity[1], Black, Text[Style[#, 16, Bold], #2] & @@@ {{"x", {3.1, 0, 0}}, {"y", {0, 2.2, 0}}, {"z", {0, 0, 2.2}}} }, opts]; pln = ContourPlot3D[ p.{x + a, y, z} == 0, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, Mesh -> None, ContourStyle -> {Opacity[.5], Darker@Green}]; Show[elp, pln, Graphics3D[{ Text[Style[ StringJoin["V= ", ToString[N@volume[p, {a, b, c}]/Pi], "\[Pi]"], 20], {2, 2, 0}]}]], Pane[Row[{ Pane[direction3D[Dynamic[p]], ImageMargins -> {{0, 40}, {0, 0}}], Pane[Column[{Control@{{a, 2}, 1, 2}, Control@{b, 1, 2}, Control@{c, 1, 2}}]]}]]]

Share Improve this answer Follow edited Nov 9, 2014 at 0:07 answered Oct 12, 2014 at 16:08 Junho LeeJunho Lee 5,17511 gold badge1616 silver badges3333 bronze badges $\endgroup$ 1

• $\begingroup$ thats a nice visualization! thank you $\endgroup$ – Corse Commented Oct 13, 2014 at 1:44

Add a comment | 6 $\begingroup$

Functions for generating ellipsoid $x^2/a^2+y^2/b^2+z^2/c^2=1$ and plane through point $\vec{p}$ with normal $\vec{n}$, i.e $\vec{n}\cdot(\vec{r}-\vec{p})=0$

el[a_, b_, c_, x_, y_, z_] := x^2/a^2 + y^2/b^2 + z^2/c^2 pl[n_, p_, x_, y_, z_] := z /. First@Solve[n.({x, y, z} - p) == 0, z]

Using as an example: a=1, b=2,c=3 and normal {1,0,par}:

Manipulate[ Column[{Show[ Plot3D[Evaluate[pl[{1, 0, par}, {-1, 0, 0}, x, y, z]], {x, -3, 3}, {y, -3, 3}, Mesh -> None, PlotStyle -> Opacity[0.5]], RegionPlot3D[ Evaluate[ reg = (el[1, 2, 3, x, y, z] < 1 && pl[{1, 0, par}, {-1, 0, 0}, x, y, z] > z)], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, BoxRatios -> {1, 1, 1}, Mesh -> False, PlotStyle -> Red, PerformanceGoal -> "Quality"], RegionPlot3D[ Evaluate[el[1, 2, 3, x, y, z] < 0.9], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, PlotStyle -> Directive[Blue, Opacity[0.2]], BoxRatios -> {1, 1, 1}, Mesh -> None], PlotRange -> Table[{-3, 3}, {3}], ImageSize -> 400], StringForm["Volume of red region:`1`", NumberForm[RegionMeasure[ImplicitRegion[reg, {x, y, z}]], 3]], StringForm["Volume of ellipsoid:`1`", NumberForm[N@4 Pi 1 2 3/3, 3]] }] , {par, 0.5, 4}]

UPDATE

In relation to comment:

rot[a_, p_, x_, y_, z_] := pl[{Sin[a], 0, Cos[a]}, p, x, y, z]

This interactive graphic shows relation of normal to plane and what I understand is the desired angle from horizontal plane:

Manipulate[ Show[RegionPlot3D[ Evaluate[el[1, 2, 3, x, y, z] < 1], {x, -4, 4}, {y, -4, 4}, {z, -4, 4}, Mesh -> False, PlotStyle -> Opacity[0.3]], Plot3D[rot[a Degree, {-1, 0, 0}, x, y, z], {x, -4, 4}, {y, -4, 4}, Mesh -> False, PlotStyle -> {Blue, Opacity[0.5]}], Graphics3D[{{Arrow[2 {{0, 0, 0}, {0, 0, 1}}]}, {Red, Arrow[2 {{0, 0, 0}, {Sin[a Degree], 0, Cos[a Degree]}}] }, {Arrow[{{-1, 0, 0}, {1, 0, 0}}]}, {Arrow[{{-1, 0, 0}, {-1 + 2 Cos[a Degree], 0, -2 Sin[a Degree]}}]} }]], {a, 0, 90, Appearance -> "Labeled"}]

Share Improve this answer Follow edited Oct 10, 2014 at 9:15 answered Oct 10, 2014 at 8:04 ubpdqnubpdqn 64.4k33 gold badges6565 silver badges153153 bronze badges $\endgroup$ 10

• $\begingroup$ hi there ubpdqn, thank you for the great implementation! Just wondering, how would I be able to obtain the volume of the red region in a generalized equation based on a specified rotation angle? I would still need to go through and fully understand your code there. I'm still not too clear as to how you managed to defined the plane (using vector product) and rotation axis. $\endgroup$ – Corse Commented Oct 10, 2014 at 8:43
• $\begingroup$ @Corse see update in relation to angles $\endgroup$ – ubpdqn Commented Oct 10, 2014 at 9:16
• $\begingroup$ you have addressed my exact problem!! but I have a couple of queries I hope you could advise: 1) What is the term 'par' supposed to represent and how do I convert it to the rotation angle 'a' ? 2) When I completely set the value of 'par' to the extreme right, which gives a value of par=4, the corresponding volume of the red region is 11 which is less than half of the volume of ellipsoid. How do I edit par to obtain the correct limits? (i.e. volume for a horizontal and a vertical plane) $\endgroup$ – Corse Commented Oct 10, 2014 at 11:28
• $\begingroup$ Also, to edit the geometry of the ellipse, I would have to edit the values of a_,b_,c_ for each occurence in the el[] function right? $\endgroup$ – Corse Commented Oct 10, 2014 at 11:29
• $\begingroup$ @Corse par just varies the normal vector that defines the plane (1,0,par). You could use the second or updated code rot to define the inclined plane and the relevant inequalities. Yes, I merely chose a=1,b=2,c=3 for illustrative purposes. You can choose what you like. Finally please also see my answer to your question re: rotation about arbitrary axis:mathematica.stackexchange.com/a/61782/1997 $\endgroup$ – ubpdqn Commented Oct 10, 2014 at 12:09

| Show 5 more comments

Your Answer

Thanks for contributing an answer to Mathematica Stack Exchange!

• Please be sure to answer the question. Provide details and share your research!

But avoid …

• Asking for help, clarification, or responding to other answers.
• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

Draft saved Draft discarded

Sign up or log in

Sign up using Google Sign up using Email and Password Submit

Post as a guest

Name Email

Required, but never shown

Post Your Answer Discard

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

• calculus-and-analysis
• geometry
• regions

or ask your own question.

• The Overflow Blog
• We'll Be In Touch - A New Podcast From Stack Overflow!
• The app that fights for your data privacy rights
• Featured on Meta
• More network sites to see advertising test
• We’re (finally!) going to the cloud!

Linked

1 Find volume of a translated and rotated ellipsoid below xy (horizontal) plane 6 Equation of a transformed (rotated) ellipsoid 1 Rotating a plane about an arbitrary axis 3 A triple integral involving Abs over an ellipsoidal region

Related

3 Finding the area bounded by a logspiral curve and two straight lines 74 Can Mathematica solve Plateau's problem (finding a minimal surface with specified boundary)? 5 Another Volume by ImplicitRegion and RegionPlot3D question 2 Finding the volume enclosed by two surfaces of revolution 1 RegionProduct to extrude a 3d volume with thickness variation 2 Finding Ellipse Axes 5 Creating a Domain of Regions Defined by Random, Non-Intersecting Curves 6 Techniques for NIntegrate on a 3D surface region

Hot Network Questions

• What happens if a client involved in active litigation disappears?
• Why aren't there square astronomical units or square light years?
• Sums and Products of Adjacent Numbers 2
• USB drives in space?
• eLife-like publications and Tenure Decisions
• Tate's thesis paper
• Why aren't there any large Ultraviolet (UV) space telescopes?
• Topology for the complex numbers endowed with different notions of convergence
• How can I tell if commercial packaging is suitable for Sous Vide cooking?
• Set Position Node
• Do rediscoveries justify publication?
• Future inextendible curve
• Draw a 3D plot of beams of an antenna using Pgfplots or any other package
• Is partial correctness decidable?
• How would Merfolk make a solar oven
• Can I protect my EV car charger's cable with aluminum tape and a stainless steel hose helix?
• What is it called when you have a hobby where you're good enough at to impress others but you yourself know you're only beginning?
• Can I use integrated cable stop from clamp-on downtube shifter in thumb shifter conversion?
• Suggestion for catching a flight with short layover in Amsterdam
• What is the probability that there will be at least one empty box?
• Why doesn't Hotelling's law seem to apply to the stances taken by political parties?
• Top loading an Aircraft vs bottom Loading
• Can I pretend that Spearman's = Pearson's correlation coefficients for meta analysis?
• Where in the Gospels does Jesus explain what 'The Truth' is?

more hot questions Question feed Subscribe to RSS Question feed

To subscribe to this RSS feed, copy and paste this URL into your RSS reader.

lang-mma