Finding length of intersection of two surfaces

Question

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$

Answer 1

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.

Answer 2

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.

Let's start with parametrizing the curve

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]

Answer 3

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

Answer 4

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

Answer 5

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 *)

Answer 6

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 *)