3
$\begingroup$

I'd like to calculate the volume of the intersection and the union of

R1 = ImplicitRegion[x^2 + y^2 + z^2 <= 9, {x, y, z}];
R2 = ImplicitRegion[(-1/2 + x)^2 + (y + 1/3)^2 + (z - 8)^2 <= 49, {x, y, z}];

In version 13 on Windows 10

Volume[RegionIntersection[R1, R2]]

22.5624

Volume[RegionUnion[R1, R2]]

1527.29

The above numerical results are sufficient in most of the cases, but not in all the cases. Not only for sportive interest (e.g. see that thread), I want to find exact results. Of course, I am aware of a spherical cap.

I find the only related question in this forum. However, the advice of @J. M.'s slightly less busy there to use

CylindricalDecomposition[x^2+y^2+z^2<=9&&(-1/2+x)^2+(y+1/3)^2+(z-8)^2<=49,{x,y,z}]

(x == (2631 - 2 Sqrt[1288846631])/27804 && y == (-877 + 36 x)/13848 - 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] && z == Sqrt[9 - x^2 - y^2]) || ((2631 - 2 Sqrt[1288846631])/27804 < x < (2631 + 2 Sqrt[1288846631])/ 27804 && ((y == (-877 + 36 x)/13848 - 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] && z == Sqrt[ 9 - x^2 - y^2]) || ((-877 + 36 x)/13848 - 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] < y < (-877 + 36 x)/13848 + 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] && 8 - 1/6 Sqrt[1751 + 36 x - 36 x^2 - 24 y - 36 y^2] <= z <= Sqrt[9 - x^2 - y^2]) || (y == (-877 + 36 x)/13848 + 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] && z == Sqrt[9 - x^2 - y^2]))) || (x == ( 2631 + 2 Sqrt[1288846631])/27804 && y == (-877 + 36 x)/13848 - 1/577 Sqrt[2222039 + 63144 x - 333648 x^2] && z == Sqrt[9 - x^2 - y^2])

seems to be only good wishes: I have strong doubts whether Mathematica is able to calculate those triple integrals in cartesian coordinates.

So the question is, how to exactly calculate these volumes?

Addition. In view of a valuable comment of @J. M.'s slightly less busy

R3 = ImplicitRegion[(1/2 + x)^2 + (y - 1/3)^2 + (z - 8)^2 <= 49, {x, 
y, z}];
Volume[RegionUnion[RegionUnion[R1, R2], R3]]

1709.77

in a long time.

PS. The command

Volume[CSGRegion[ "Union", {CSGRegion[
 "Union", {Ball[{0, 0, 0}, 3], Ball[{1/2, -1/3, 8}, 7]}], 
Ball[{-1/2, 1/3, 8}, 7]}]] // Timing

crashes my comp in few dozen minutes.

$\endgroup$
10
  • 1
    $\begingroup$ Nowadays, one would do something like Volume[CSGRegion["Union", {Ball[{0, 0, 0}, 3], Ball[{1/2, -1/3, 8}, 7]}]] or Volume[CSGRegion["Intersection", {Ball[{0, 0, 0}, 3], Ball[{1/2, -1/3, 8}, 7]}]]. Unfortunately, this does not seem to work yet for more than two spheres. The MO thread you link to is a little related to why I asked this previous question. $\endgroup$ May 30, 2022 at 8:12
  • $\begingroup$ @J. M.: Thank you so much. Please, convert your comment to an answer. $\endgroup$
    – user64494
    May 30, 2022 at 8:19
  • $\begingroup$ For older version without CSGRegion use Region* commands; see image. But what if an implicit region didn't fit a named primitive? $\endgroup$
    – Syed
    May 30, 2022 at 8:28
  • 1
    $\begingroup$ Have you checked using geometric primitives instead? Like {Volume@RegionUnion[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]], Volume@RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]]} $\endgroup$ May 30, 2022 at 9:13
  • 2
    $\begingroup$ Dear colleagues: I prefer answers over valuable comments. $\endgroup$
    – user64494
    May 30, 2022 at 9:21

2 Answers 2

3
$\begingroup$

Get the anylytical volume of intersection with variable transformation.

Transform center of R2 ball into z direction

R1[{x_, y_, z_}] = x^2 + y^2 + z^2 <= 9;
R2[{x_, y_, z_}] = (-1/2 + x)^2 + (y + 1/3)^2 + (z - 8)^2 <= 49;

rot = RotationTransform[{{1/2, -1/3, 8}, {0, 0, 1}}]

(*   TransformationFunction[{{(4*(2317 + 108*Sqrt[2317]))/30121, 
    -((6*(-2317 + 48*Sqrt[2317]))/30121), -(3/Sqrt[2317]), 0}, 
   {-((6*(-2317 + 48*Sqrt[2317]))/30121), 
    (3*(6951 + 64*Sqrt[2317]))/30121, 2/Sqrt[2317], 0}, 
   {3/Sqrt[2317], -(2/Sqrt[2317]), 48/Sqrt[2317], 0}, {0, 0, 0, 1}}]   *)

Due to symmetrie R1 coordinates are transformed, but ball equation remains the same.

rot[{x, y, z}]

(*   {(4 (2317 + 108 Sqrt[2317]) x)/30121 - (6 (-2317 + 48 Sqrt[2317]) y)/
  30121 - (3 z)/Sqrt[2317], -((6 (-2317 + 48 Sqrt[2317]) x)/30121) + (
  3 (6951 + 64 Sqrt[2317]) y)/30121 + (2 z)/Sqrt[2317], (3 x)/Sqrt[
  2317] - (2 y)/Sqrt[2317] + (48 z)/Sqrt[2317]}   *)

Total@(rot[{x, y, z}]^2) // Simplify

(*   x^2 + y^2 + z^2   *)

Center of R2 gets stretched a little bit in z direction. Use this new center, radius remains the same.

{rr = rot[{1/2, -1/3, 8}] // Simplify, rr // N}

(*   {{0, 0, Sqrt[2317]/6}, {0., 0., 8.02254}}   *)

red = Reduce[
  x^2 + y^2 + z^2 <= 9 && x^2 + y^2 + (z - Sqrt[2317]/6)^2 <= 49]

{int = Integrate[Boole[LogicalExpand@red], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}]//Simplify, 
  int // N} // AbsoluteTiming

{4.0469006, {(740/3 - 29879687/(2592 Sqrt[2317])) \[Pi], 22.5624}}
$\endgroup$
7
  • $\begingroup$ +1. How about three balls? In the article cited in that thread the statement "Theoretically, we establish a formula for the volume of a restriction, based on Gauss' divergence theorem. The proof being constructive, we develop the associated algorithm" is claimed. $\endgroup$
    – user64494
    May 30, 2022 at 13:50
  • $\begingroup$ "Three balls?" I don't know. Could be a distorted, non-rectangular coordinate system. But I think Integrate would get stuck then. $\endgroup$
    – Akku14
    May 30, 2022 at 14:33
  • $\begingroup$ I'd like to quote Joseph O'Rourke from that thread in my question "This impressive work was tested on ~60,000 models in the Protein Data Bank, computing the volume of the union of sometimes more than 50,000 atoms represented as balls of different radii. Implemented in C++". $\endgroup$
    – user64494
    May 30, 2022 at 14:59
  • $\begingroup$ @user: none of the built-ins are yet capable of handling more than two balls, as I alluded previously. If this is important enough for you, you should bug Support so they'd have an idea about the demand. $\endgroup$ May 30, 2022 at 16:02
  • $\begingroup$ Accepted. I can reproduce it in version 13 on Windows 10. $\endgroup$
    – user64494
    May 30, 2022 at 17:23
4
$\begingroup$

It seems like WL has some built-in answers for symbolic regions:

Volume @ RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]]

([Pi] (r1 + r2 - Sqrt[z^2])^2 (-3 r1^2 + 6 r1 r2 - 3 r2^2 + z^2 + 2 r1 Sqrt[z^2] + 2 r2 Sqrt[z^2]))/(12 Sqrt[z^2])

However, be careful because it seems like these formulas make some assumptions about the region existing in the first place. For example:

Volume @ RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]] /. {r1 -> 1, r2 -> 1, z -> 100}

(249704 [Pi])/3

This should obviously be zero (and it is, if you do the substitution before calculating the volume):

Volume[
 RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]] /. {r1 -> 1, r2 -> 1, z -> 100}
 ]

0

Edit

As halmir pointed out in the comments, you can use GenerateConditions to get a more reliable result:

Volume[
 RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]], 
 GenerateConditions -> True
]
$\endgroup$
3
  • $\begingroup$ +1. Thank you so much for your answer and analysis of the result. $\endgroup$
    – user64494
    May 30, 2022 at 15:12
  • 3
    $\begingroup$ to see the condition, Volume[RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]], GenerateConditions -> True] $\endgroup$
    – halmir
    May 31, 2022 at 15:55
  • $\begingroup$ Even Volume[RegionIntersection[Ball[{x1, y1, z1}, r1], Ball[{x2, y2, z2}, r2]], GenerateConditions -> True] and Volume[RegionUnion[Ball[{x1,y1,z1},r1],Ball[{x2,y2,z2},r2]],GenerateConditions->True] work very quickly. I think these formulas are implemented in Mathematica as table values. $\endgroup$
    – user64494
    May 31, 2022 at 16:34

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.