# Finding length of intersection of two surfaces

I would like to know how we find the length of the intersection of two surfaces. For instance, in the following example,a surface intersects with a plane:

How do we find the length of intersection which is highlighted in the region?

Is there a general method for computing the length in a given region for a surface and a plane?

Here is the code for the above example:

h = Exp[-x^3 - y] - 1 - t ;
g = y - t ;

ContourPlot3D[{h == 0, g == 0}, {x, -4, 4}, {y, -4, 4}, {t, -4, 4},
MeshFunctions -> { Function [{x, y, t, f}, h - g]},
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}},
ContourStyle ->
Directive[Orange, Opacity[0.3], Specularity[White, 30]]]


It is the intersection of the function $e^{-x^3-y}-1=z$ and $y=z$

• @RunnyKine True enough, and thanks for the notice. I have to say given how slow and how imprecise the Discretize* functions still are, they're not very good options for numerics except for very rough approximations. – Michael E2 Sep 19 '14 at 22:33
• @MichaelE2. I didn't get a notice of your comment. You're right about the slowness, plus they appear to be still buggy even in 10.0.1 – RunnyKine Sep 21 '14 at 20:22
• @RunnyKine My feeling is that they are pretty good up to codimension 1 - i.e. one dimension less than the ambient space. So boundaries of solids in space, boundaries of areas in the plane. Not so good with lines in space, though. And you usually get Plot-type precision, when discretized. It's nice to able to set up code in these terms though, say for graphical illustrations. – Michael E2 Sep 21 '14 at 20:27
• @MichaelE2, the bug about MaxCellMeasure not working wasn't fixed in 10.0.1, that would have helped a bit with precision. – RunnyKine Sep 21 '14 at 20:32

Fixed (see below)

Here's an approach:

r1 = Exp[-x^3 - y] - 1 == z;
r2 = y == z;


We create ImplicitRegions:

reg1 = ImplicitRegion[r1, {x, y, z}];
reg2 = ImplicitRegion[r2, {x, y, z}];


The intersection of these regions is the line you seek:

reg = RegionIntersection[reg1, reg2];


And here is the length (note the inclusion of the range of values in DiscretizeRegion)

RegionMeasure @ DiscretizeRegion[reg, {{-4, 4}, {-4, 4}, {-4, 4}}]


OR

ArcLength @ DiscretizeRegion[reg, {{-4, 4}, {-4, 4}, {-4, 4}}]


10.9488106

Note: Previously the approach shown gave a wrong answer as a result of a bug previously thought to be from DiscretizeRegion but I'm not so sure about that now. Anyways the new approach shown here has fixed that issue.

• How come the RegionPlot of a space curve is a 2D plot? – user484 Sep 3 '14 at 13:58
• @RahulNarain since one of the surfaces is a plane, the result is two dimensional on that plane. – Sparr Sep 3 '14 at 13:59
• The numerical answer here differs from @Pickett's numerical answer by a small but notable amount. – Sparr Sep 3 '14 at 13:59
• I wouldn't say "slightly affected" since the result differs from the symbolic one by 1 - 9.46689315/(10.9513) i.e. more than $10$%. – Artes Sep 3 '14 at 15:58
• The update to V10.0.1 has improved the result to 3-4 digits of precision, which is still not as good an approximation as the length of the line produced by ContourPlot3D. (It takes a long time to compute, though. Not much bang for the buck, all in all.) – Michael E2 Sep 19 '14 at 22:37

Since another answer provides misleading results I'm motivated to explain what's wrong there and suggest a correct solution. One can guess that the error there has been produced by DiscretizeRegion or by an improper parametrization of the curve with RegionIntersection (more likely the former reason).

Let's start from the begining defining the following functions:

h[x_, y_, t_] := Exp[-x^3 - y] - 1 - t
g[x_, y_, t_] := y - t


I'd like a neat method involving ImplicitRegion, RegionIntersection however it's rather surprising that RegionMeasure[reg] (defined in another answer) yields this warning:

 RegionMeasure::nmet: Unable to compute the measure of region
RegionIntersection[
ImplicitRegion[-1+E^(Times[<<2>>]+Times[<<2>>])-t==0 &&
-4<=x<=4&&-4<=y<=4&&-4<=t<=4, {x,y,t}],
ImplicitRegion[ y-t==0 &&-4<=x<=4&&-4<=y<=4&&-4<=t<=4, {x,y,t}]]. >>


and does not provide a symbolic result. Then DiscretizeRegion produces after a lapse of a few seconds 9.46689315 which seems slightly too small.

sol = Solve[ h[x, y, t] == 0 && g[x, y, t] == 0 &&
-4 <= x <= 4 && -4 <= y <= 4 && -4 <= t <= 4, {x, y, t},
MaxExtraConditions -> All]

Solve::svars: Equations may not give solutions for all "solve" variables. >>

{{ y -> ConditionalExpression[ -1 + ProductLog[E^(1 - x^3)],
-(4 + Log[5])^(1/3) <= x <= 4],
t -> ConditionalExpression[ -1 + ProductLog[E^(1 - x^3)],
-(4 + Log[5])^(1/3) <= x <= 4]},
{x -> 0, y -> 0, t -> 0}}


Of course the second solution is a single point (it's an exceptional solution which we simply omit) and we need not worry about warning since that can be verifed by another methods like Reduce.

Now we define the curve:

curve = Simplify[{x, y, t} /. First @ sol, -(4 + Log[5])^(1/3) <= x <=


4]

Now since x is restricted to the range -(4 + Log[5])^(1/3) <= x <= 4 we can compute:

ArcLength[ curve, {x, -(4 + Log[5])^(1/3), 4}]


N @ %

10.9513


This reuslt can be verified by

NIntegrate[ Sqrt[ 1 + 2 (-((3 x^2 ProductLog[E^(1 - x^3)])/(
1 + ProductLog[E^(1 - x^3)])))^2], {x, -(4 + Log[5])^(1/3), 4}]


We have put

D[-1 + ProductLog[E^(1 - x^3)], x]


for y'[x] and t'[x] to

ArcLength[{x, y[x], t[x]}, {x, x0, x1}]

Integrate[ Sqrt[ 1 + y'[x]^2 + t'[x]^2], {x, x0, x1}, Assumptions -> x0 < x1]


Brute force approach:

     g[y0_] :=
x /. FindRoot[ Exp[-(x^3 + y)] - 1 == y /. y -> y0 , {x, -4, 4}]
line = Table[{g[y], y, y}, {y, -1, 4, .0001}];
Graphics3D[Line@line]
Total[Norm@(Subtract @@ #) & /@ Partition[line, 2, 1]]


10.9513

• On the contrary, to me your approach is quite elegant, concise with a flair for mathematics. And it works OK at least since MMA V5. – Sigis K Sep 26 '14 at 11:45

Thanks to Artes for helping me sort out bugs in this answer.

Solve[E^(-x^3 - y) - 1 - y == 0, y]


{{y -> -1 + ProductLog[E^(1 - x^3)]}}

To account for that fact that $-4 < y = z < 4$ we need to find the range of allowed $x$ values.

FindRoot[-1 + ProductLog[E^(1 - x^3)] == 4, {x, 1}]


{x -> -1.77681}

FindRoot[-1 + ProductLog[E^(1 - x^3)] == -4, {x, 1}]


{x -> 4.2856}

4.28 is larger than the upper limit for $x$, so in this case we keep 4. Finally we obtain:

ArcLength[{x, -1 + ProductLog[E^(1 - x^3)], -1 + ProductLog[E^(1 - x^3)]}, {x, -1.7768050578807908, 4}]


10.9513

• The numerical answer here differs from @RunnyKine's numerical answer by a small but notable amount. – Sparr Sep 3 '14 at 14:00
• @Sparr Bugs fixed and the new result is correct. – C. E. Sep 3 '14 at 21:08

A more geometrical approach based on CP3D surface-surface intersections boundary style...

inter = ContourPlot3D[{h == 0, g == 0}, {x, -4, 4}, {y, -4,
4}, {t, -4, 4}, Mesh -> None,
ContourStyle ->
Directive[Orange, Opacity[0.3], Specularity[White, 30]],
BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> Blue}]

reg = DiscretizeGraphics[inter, Line]

RegionMeasure[reg]

(* 10.9518 *)


Rough approximation via plotting

V9: One can estimate the length of the polygonal path in the ContourPlot from the graphics. It's a little easier to process if I adapt Daniel Lichtblau's (undocumented?) use of BoundaryStyle in his answer to Plotting implicitly-defined space curves. If the points in the plot lie exactly on the intersection, the length should be a lower bound on the true arc length. Since the points are not, we should expect some error.

contourplot =
ContourPlot3D[{h, g}, {x, -4, 4}, {y, -4, 4}, {t, -4, 4}, Contours -> {0},
ContourStyle -> Opacity[0], Mesh -> None,
BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {{Black}}},
Boxed -> False]


Norm /@ Differences@
First@Cases[Normal@contourplot, Line[p_] :> p, Infinity] // Total

(* 10.9512 *)


V10: The above does not work in V10. One can use {g == 0, h == 0} instead of Contours -> {0}, but the performance is about ten times slower (~11 sec. in V10) than in V9. Here's similar way that is faster:

Needs["GeneralUtilities"];

contourplot =
ContourPlot3D[{g == 0}, {x, -4, 4}, {y, -4, 4}, {t, -4, 4},
PlotPoints -> 35, MaxRecursion -> 4,
MeshFunctions -> {Function[{x, y, t}, Evaluate@h]}, Mesh -> {{0}},
ContourStyle -> None, BoundaryStyle -> None, Boxed -> False,
Axes -> None,
WorkingPrecision -> \$MachinePrecision]; // AccurateTiming

Norm /@ Differences @ First @ Cases[Normal@contourplot, Line[p_] :> p, Infinity] // Total
(*
0.17347  (* timing *)

10.9508  (* length *)
*)


By comparison DiscretizeRegion@RegionIntersection[ImplicitRegion[g == 0,..], ...] takes an amazing 20 seconds. Probably these function will be developed for general usefulness in the future.

Integration via NDSolve

With NDSolve we can integrate the speed of a parametrization over the segment of the curve in the plot region. One thing we need is a starting point on the curve. If Solve fails to give us one, we can get them from the contour plot in this way:

endpoints = First@Cases[Normal@contourplot, Line[p_] :> p[[{1, -1}]], Infinity]

(* {{4., -1.00395, -1.00125}, {-1.77835, 4., 4.}} *)


In this case, Solve works, and we'll use the first point where x == 4. (One could also use FindRoot.) Another we need is a starting velocity that points into the plot region. We can get a tangent vector from the cross product of the gradient vectors of the two surfaces. We will set up a differential equation for a unit speed parametrization, so we normalize the cross product. The direction is chosen by looking at the sign of the dot product of the tangent vector with the outward normal of the boundary of the plot region, that is, the vector {1, 0, 0} which is normal to x == 4. Finally we'll useWhenEvent to stop the integration when the parameterization leaves the plot region.

bdynormal = {1, 0, 0};                          (* start at x == 4 *)
startpoint = {x, y, t} /. First@Solve[{h == 0, g == 0, x == 4}, {x, y, t}];
startvelocity = -Sign[Dot[#, bdynormal]] Normalize[#] &@ (
Cross[D[h, {{x, y, t}}], D[g, {{x, y, t}}]] /. Thread[{x, y, t} -> startpoint]);

dae = {
h == 0 /. v : x | y | t :> v[u],             (* constrains solution to h == 0 *)
g == 0 /. v : x | y | t :> v[u],             (* constrains solution to g == 0 *)
Sqrt[#.# &@{x'[u], y'[u], t'[u]}] == 1};     (* unit speed *)
ics = {{x[0], y[0], t[0]}    == startpoint,
{x'[0], y'[0], t'[0]} == startvelocity}; (* used by NDSolve to pick which root
of the dae to follow *)
reg = Cuboid[{-4, -4, -4}, {4, 4, 4}];          (* the plot region *)

s = NDSolveValue[{
dae, ics,
arclength'[u] ==
Sqrt[#.# &@{x'[0], y'[0], t'[0]}], arclength[0] == 0,
WhenEvent[! RegionMember[reg, {x[u], y[u], t[u]}], "StopIntegration"],
arclength, {u, 0, Infinity}];
s[s["Domain"][[-1, -1]]]


Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

(* 10.9513 *)