# Volume of the intersection and union of two spheres

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)/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)/27804 < x < (2631 + 2 Sqrt)/ 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)/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.

• 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. May 30, 2022 at 8:12
• @J. M.: Thank you so much. Please, convert your comment to an answer. May 30, 2022 at 8:19
• For older version without CSGRegion use Region* commands; see image. But what if an implicit region didn't fit a named primitive?
– Syed
May 30, 2022 at 8:28
• 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]]} May 30, 2022 at 9:13
• Dear colleagues: I prefer answers over valuable comments. May 30, 2022 at 9:21

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))/30121,
-((6*(-2317 + 48*Sqrt))/30121), -(3/Sqrt), 0},
{-((6*(-2317 + 48*Sqrt))/30121),
(3*(6951 + 64*Sqrt))/30121, 2/Sqrt, 0},
{3/Sqrt, -(2/Sqrt), 48/Sqrt, 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) x)/30121 - (6 (-2317 + 48 Sqrt) y)/
30121 - (3 z)/Sqrt, -((6 (-2317 + 48 Sqrt) x)/30121) + (
3 (6951 + 64 Sqrt) y)/30121 + (2 z)/Sqrt, (3 x)/Sqrt[
2317] - (2 y)/Sqrt + (48 z)/Sqrt}   *)

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/6}, {0., 0., 8.02254}}   *)

red = Reduce[
x^2 + y^2 + z^2 <= 9 && x^2 + y^2 + (z - Sqrt/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)) \[Pi], 22.5624}}

• +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. May 30, 2022 at 13:50
• "Three balls?" I don't know. Could be a distorted, non-rectangular coordinate system. But I think Integrate would get stuck then. May 30, 2022 at 14:33
• 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++". May 30, 2022 at 14:59
• @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. May 30, 2022 at 16:02
• Accepted. I can reproduce it in version 13 on Windows 10. May 30, 2022 at 17:23

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
]

• +1. Thank you so much for your answer and analysis of the result. May 30, 2022 at 15:12
• to see the condition, Volume[RegionIntersection[Ball[{0, 0, 0}, r1], Ball[{0, 0, z}, r2]], GenerateConditions -> True] May 31, 2022 at 15:55
• 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. May 31, 2022 at 16:34