# How do I find and plot the intersection of these three surfaces?

I know how to find the region of intersection of these cylinders using integration, but how would I plot that in mathematica? This is what I have

ContourPlot3D[{x^2 + y^2 == 1, x^2 + z^2 == 1,
y^2 + z^2 == 1}, {x, -3, 3}, {y, -3, 3}, {z, -3, 3},
AxesLabel -> Automatic, PlotLegends -> "Expressions"]


• The intersection of the three surfaces would be eight points. You could use Solve[] to find them – Michael E2 Jun 22 '19 at 23:07
• Possibly related: mathematica.stackexchange.com/q/100825/4999 – Michael E2 Jun 22 '19 at 23:10
• How would I do that? by defining an InfinitePlane through the points on the cylinder? – sabunketohar Jun 23 '19 at 0:02
• Use Solve[] with the three equations – Michael E2 Jun 23 '19 at 4:50

You can use region functionality for this. Here is the intersection of the boundary of the cylinders:

reg = RegionIntersection[
RegionBoundary @ Cylinder[{{0,0,-2},{0,0,2}}],
RegionBoundary @ Cylinder[{{0,-2,0},{0,2,0}}],
RegionBoundary @ Cylinder[{{-2,0,0},{2,0,0}}]
];


Here is a visualization of the points of intersection:

mesh = DiscretizeRegion[reg, MeshCellStyle -> {0 -> Directive[Red, PointSize[Large]]}]


And the values of the intersection:

MeshPrimitives[mesh, 0]


{Point[{0.707107, 0.707107, -0.707107}], Point[{-0.707107, -0.707107, 0.707107}], Point[{0.707107, 0.707107, 0.707107}], Point[{-0.707107, -0.707107, -0.707107}], Point[{0.707107, -0.707107, 0.707107}], Point[{-0.707107, 0.707107, -0.707107}], Point[{-0.707107, 0.707107, 0.707107}], Point[{0.707107, -0.707107, -0.707107}]}

Visualization of the intersection of the solid cylinders (blue) and the intersection of the boundary of the cylinders (red):

reg2 = RegionIntersection[
Cylinder[{{0,0,-1},{0,0,1}}],
Cylinder[{{0,-1,0},{0,1,0}}],
Cylinder[{{-1,0,0},{1,0,0}}]
];

Show[Region @ reg2, mesh]


You can show the lines corresponding to intersection of pairs of three surfaces using the option BoundaryStyle. You can get the points where all three surfaces intersect using NSolve (as Michael suggested in comments) and use them with Graphics3D and put the two graphics together using Show:

ClearAll[x, y, z]
pnts = {x, y, z} /. NSolve[{x^2 + y^2 == 1, x^2 + z^2 == 1, y^2 + z^2 == 1}, {x, y, z}];

Show[ContourPlot3D[{x^2 + y^2 == 1, x^2 + z^2 == 1, y^2 + z^2 == 1},
{x, -3, 3}, {y, -3, 3}, {z, -3, 3},
AxesLabel -> Automatic,
PlotLegends -> {Placed[SwatchLegend[ "Expressions"], Right],
Placed[LineLegend[{Orange, Purple, Cyan},
{"intersection(1, 2)", "intersection(1, 3)", "intersection(2, 3)"}], Right],
Placed[PointLegend[{Red}, {"intersection(1, 2, 3)"}], Right]},
Mesh -> None,
ContourStyle -> Opacity[.4],
BoundaryStyle -> {{1, 2} -> Directive[Orange, Thick],
{1, 3} -> Directive[Purple, Thick],
{2, 3} -> Directive[Cyan, Thick]}],
Graphics3D[{Red, Sphere[#, .2] & /@ pnts}]]


RegionPlot3D[
x^2 + y^2 < 1 &&
x^2 + z^2 < 1 &&
y^2 + z^2 < 1,
{x, -1.2, 1.2},
{y, -1.2, 1.2},
{z, -1.2, 1.2},
PlotPoints -> 100,
PlotStyle -> Opacity[0.5]]