5
$\begingroup$

I would like to find the intersection points of two surfaces which are constructed with ListPlot3D. The first list is

x={{0.687522, 0.672048, 0.657758, 0.644584, 0.632436, 0.621219, 0.61084,
   0.601212, 0.592255, 0.583898, 0.57608}, {0.718239, 0.70296, 
  0.688594, 0.675148, 0.662596, 0.650892, 0.639979, 0.629794, 
  0.620276, 0.611365, 0.603006}, {0.744103, 0.729406, 0.71536, 
  0.70202, 0.689407, 0.677517, 0.666327, 0.655806, 0.645912, 0.636603,
   0.627836}, {0.76593, 0.752001, 0.738512, 0.725536, 0.713119, 
  0.701286, 0.690045, 0.679387, 0.669294, 0.65974, 
  0.650698}, {0.784468, 0.771368, 0.75855, 0.746088, 0.734038, 
  0.722443, 0.711325, 0.700697, 0.690559, 0.680901, 
  0.671708}, {0.800344, 0.788064, 0.775953, 0.764077, 0.752495, 
  0.741253, 0.730386, 0.719917, 0.709859, 0.700215, 
  0.690982}, {0.814062, 0.802558, 0.791142, 0.779875, 0.768808, 
  0.757989, 0.747455, 0.737236, 0.727353, 0.717818, 
  0.708637}, {0.826019, 0.815231, 0.804478, 0.793809, 0.783271, 
  0.772908, 0.762756, 0.752847, 0.743206, 0.733851, 
  0.724796}, {0.836524, 0.826392, 0.816257, 0.806161, 0.796144, 
  0.786244, 0.776496, 0.766933, 0.757579, 0.748456, 
  0.73958}, {0.845825, 0.836288, 0.826724, 0.817166, 0.807648, 
  0.798205, 0.788868, 0.779666, 0.770625, 0.761767, 
  0.75311}, {0.854114, 0.845118, 0.836077, 0.82702, 0.817975, 
  0.808973, 0.80004, 0.791203, 0.782487, 0.773913, 0.765501}}

The second list is :

y={{0.86018, 0.823767, 0.783404, 0.740538, 0.697531, 0.656804, 0.619925,
   0.58739, 0.558988, 0.534209, 0.512495}, {0.868675, 0.835506, 
  0.798953, 0.760038, 0.720471, 0.682171, 0.646621, 0.61455, 0.586049,
   0.560854, 0.538561}, {0.876074, 0.845573, 0.812132, 0.776497, 
  0.739939, 0.703976, 0.66993, 0.638615, 0.610325, 0.584988, 
  0.562348}, {0.882593, 0.854335, 0.82349, 0.79062, 0.756691, 0.72291,
   0.690421, 0.660043, 0.632184, 0.606922, 0.584127}, {0.888395, 
  0.862052, 0.83341, 0.802903, 0.771278, 0.739504, 0.708556, 0.679212,
   0.65194, 0.62692, 0.604125}, {0.893601, 0.868915, 0.842169, 
  0.813704, 0.78411, 0.754168, 0.724707, 0.69644, 0.669857, 0.645204, 
  0.622538}, {0.898304, 0.875068, 0.849974, 0.823293, 0.795497, 
  0.767223, 0.739173, 0.711991, 0.68616, 0.661967, 
  0.63953}, {0.902579, 0.880625, 0.856984, 0.831874, 0.805679, 
  0.778923, 0.752202, 0.726087, 0.701042, 0.677375, 
  0.655245}, {0.906485, 0.885673, 0.86332, 0.839607, 0.814846, 
  0.789472, 0.763993, 0.738914, 0.714668, 0.691571, 
  0.66981}, {0.91007, 0.890283, 0.869082, 0.846617, 0.823147, 
  0.799035, 0.774713, 0.750629, 0.727182, 0.704681, 
  0.683333}, {0.913376, 0.894513, 0.874348, 0.853006, 0.830703, 
  0.807745, 0.784502, 0.761367, 0.738705, 0.716816, 0.695914}}

I join them with Join[{x}, {y}] and the ListPlot3D looks like:

enter image description here

I would like to find the intersection points of those two surfaces.

$\endgroup$
  • 1
    $\begingroup$ What is it that you really want? If you are talking about the 121 points in each x and y, then there are no points common for both: ContainsAny[Catenate@x, Catenate@y] yields False. I guess you'd want an Interpolation of the surfaces and their (numerical) intersection. But please explain in what form do you want it: a list of 3D points, a function that gives such point given one of the coordinates, an InterpolatingFunction describing the curve in 3D, an analytical approximation? $\endgroup$ – corey979 Jan 8 '17 at 19:00
  • $\begingroup$ Ye there are no point in common. I want the numerical intersection of the interpolated surfaces. I just want to see the values of k and l for each intersection point (x^*). A list of 3D points can be usefull I guess. $\endgroup$ – ecco Jan 8 '17 at 19:20
  • $\begingroup$ There are already good answers, but here is a way to get approximate common values. In[383]:= fullx = Flatten[MapIndexed[Flatten[{#2, #1}] &, x, {2}], 1]; fully = Flatten[MapIndexed[Flatten[{#2, #1}] &, y, {2}], 1]; nfx = Nearest[fullx]; nbrs = Map[{#, First[nfx[#]]} &, fully]; close = Select[nbrs, Norm[#[[1]] - #[[2]]] < 1/100 &][[All, 1]] Out[387]= {{1, 7, 0.619925}, {2, 7, 0.646621}, {3, 7, 0.66993}, {4, 7, 0.690421}, {5, 7, 0.708556}, {6, 7, 0.724707}, {7, 6, 0.767223}, {7, 7, 0.739173}, {8, 6, 0.778923}, {9, 6, 0.789472}, {10, 6, 0.799035}, {11, 6, 0.807745}} $\endgroup$ – Daniel Lichtblau Jan 8 '17 at 22:02
  • $\begingroup$ ListPlot3D[Join[{x}, {y}], Axes -> True] gives the range {1, 11}; your plot is inconsistent with your description. $\endgroup$ – corey979 Jan 9 '17 at 8:58
6
$\begingroup$

Using an approach first presented in this answer:

xf = ListInterpolation[x]; yf = ListInterpolation[y];
g1 = Plot3D[{xf[x, y], yf[x, y]}, {x, 1, 11}, {y, 1, 11}, 
            AxesLabel -> {"a", "b", "z"}, Mesh -> None, PlotStyle -> Opacity[0.5]];
g2 = Plot3D[{xf[x, y], yf[x, y]}, {x, 1, 11}, {y, 1, 11}, 
            BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {Directive[Thick, Red]}},
            Mesh -> None, PlotStyle -> None];
Show[g1, g2]

intersection of two interpolating functions

$\endgroup$
4
$\begingroup$

Let

data = {x}~Join~{y};

Create integer coordinates:

coords = CoordinateBoundsArray[{{1, Length@x}, {1, Length@y}}];

Next we join the data with the coordinates for the interpolation:

plane1 = Interpolation[Join @@ (Transpose /@ Transpose[{coords, x}])]
plane2 = Interpolation[Join @@ (Transpose /@ Transpose[{coords, y}])]

Their plot:

plot = Plot3D[{plane1[a, b], plane2[a, b]}, {a, 1, Length@x}, {b, 1, Length@y},
  AxesLabel -> {"a", "b", "z"}, Mesh -> None, PlotStyle -> Opacity[0.5]];

Some of the absolute options will be useful:

pr = AbsoluteOptions[plot, PlotRange][[1, 2]]
box = AbsoluteOptions[plot, BoxRatios][[1, 2]]

Then, a function to find the b coordinate for each a coordinate (starting value based on plot):

root[a0_] := With[{a = a0}, b /. FindRoot[plane1[a, b] == plane2[a, b], {b, 6.}]]

Their plot:

Plot[root[a0], {a0, 1, 11}, AxesLabel -> {"a", "b"}]

enter image description here

The 3D intersection is a parametric line:

line = ParametricPlot3D[{a0, root[a0], plane1[a0, root[a0]]}, {a0, 1, 11},
  AxesLabel -> {"a", "b", "z"}, PlotRange -> pr, BoxRatios -> box, PlotStyle -> {Red, Thick}]

enter image description here

lying on the intersection of the planes:

Show[plot, line]

enter image description here

$\endgroup$
2
$\begingroup$

I am going to work in the $k$-$l$-$x^*$ space you seem to prefer.

Here are the two sets of $x^*$ points, which I am presuming are evaluated on the grid

CoordinateBoundingBoxArray[{{0., 0.}, {1., 1.}}, .1]

That is, the grid with $k$ runs for 0 to 1 in steps of 0.1 and $l$ does the same.

xStar1 =
  {{0.687522, 0.672048, 0.657758, 0.644584, 0.632436, 0.621219, 
    0.61084, 0.601212, 0.592255, 0.583898, 0.57608}, 
   {0.718239, 0.70296, 0.688594, 0.675148, 0.662596, 0.650892, 
    0.639979, 0.629794, 0.620276, 0.611365, 0.603006}, 
   {0.744103, 0.729406, 0.71536, 0.70202, 0.689407, 0.677517, 
    0.666327, 0.655806, 0.645912, 0.636603, 0.627836}, 
   {0.76593, 0.752001, 0.738512, 0.725536, 0.713119, 0.701286, 
    0.690045, 0.679387, 0.669294, 0.65974, 0.650698}, 
   {0.784468, 0.771368, 0.75855, 0.746088, 0.734038, 0.722443, 
    0.711325, 0.700697, 0.690559, 0.680901, 0.671708}, 
   {0.800344, 0.788064, 0.775953, 0.764077, 0.752495, 0.741253, 
    0.730386, 0.719917, 0.709859, 0.700215, 0.690982}, 
   {0.814062, 0.802558, 0.791142, 0.779875, 0.768808, 0.757989, 
    0.747455, 0.737236, 0.727353, 0.717818, 0.708637}, 
   {0.826019, 0.815231, 0.804478, 0.793809, 0.783271, 0.772908, 
    0.762756, 0.752847, 0.743206, 0.733851, 0.724796}, 
   {0.836524, 0.826392, 0.816257, 0.806161, 0.796144, 0.786244, 
    0.776496, 0.766933, 0.757579, 0.748456, 0.73958}, 
   {0.845825, 0.836288, 0.826724, 0.817166, 0.807648, 0.798205, 
    0.788868, 0.779666, 0.770625, 0.761767, 0.75311}, 
   {0.854114, 0.845118, 0.836077, 0.82702, 0.817975, 0.808973, 
    0.80004, 0.791203, 0.782487, 0.773913, 0.765501}};

xStar2 =
  {{0.86018, 0.823767, 0.783404, 0.740538, 0.697531, 0.656804,
    0.619925, 0.58739, 0.558988, 0.534209, 0.512495},
   {0.868675, 0.835506, 0.798953, 0.760038, 0.720471, 0.682171,
    0.646621, 0.61455, 0.586049, 0.560854, 0.538561},
   {0.876074, 0.845573, 0.812132, 0.776497, 0.739939, 0.703976,
    0.66993, 0.638615, 0.610325, 0.584988, 0.562348}, 
   {0.882593, 0.854335, 0.82349, 0.79062, 0.756691, 0.72291,
    0.690421, 0.660043, 0.632184, 0.606922, 0.584127},
   {0.888395, 0.862052, 0.83341, 0.802903, 0.771278, 0.739504,
    0.708556, 0.679212, 0.65194, 0.62692, 0.604125},
   {0.893601, 0.868915, 0.842169, 0.813704, 0.78411, 0.754168,
    0.724707, 0.69644, 0.669857, 0.645204, 0.622538},
   {0.898304, 0.875068, 0.849974, 0.823293, 0.795497, 0.767223,
    0.739173, 0.711991, 0.68616, 0.661967, 0.63953},
   {0.902579, 0.880625, 0.856984, 0.831874, 0.805679, 0.778923,
    0.752202, 0.726087, 0.701042, 0.677375, 0.655245},
   {0.906485, 0.885673, 0.86332, 0.839607, 0.814846, 0.789472,
    0.763993, 0.738914, 0.714668, 0.691571, 0.66981}, 
   {0.91007, 0.890283, 0.869082, 0.846617, 0.823147, 0.799035,
    0.774713, 0.750629, 0.727182, 0.704681, 0.683333},
   {0.913376, 0.894513, 0.874348, 0.853006, 0.830703, 0.807745, 
    0.784502, 0.761367, 0.738705, 0.716816, 0.695914}};

Now I am going to make interpolating functions for these $x^*$ values over the grid.

pts[x_] :=
  With[{grid = CoordinateBoundingBoxArray[{{0., 0.}, {1., 1.}}, .1]},
    Catenate[Transpose /@ Transpose[{grid, x}]]]

xStar1Pts = pts[xStar1]; xStar2Pts = pts[xStar2];

I now have data in a form Interpolation likes.

xStar1PtsF = Interpolation[xStar1Pts]; xStar2PtsF = Interpolation[xStar2Pts];

Now I will find the points along the intersection of xStar1PtsF and xStar2PtsF, which are the roots of xStar1PtsF - xStar2PtsF at the 11 $k$-coordinates of the grid and near $l$ = .5.

kPts = Subdivide[1., 10];
lPts = 
  (FindRoot[xStar1PtsF[#, l] - xStar2PtsF[#, l], {l, .5}] & /@ kPts)[[All, 1, 2]]
{0.637846, 0.629023, 0.616597, 0.60183, 0.585783, 0.569174, 0.552539, 
  0.536255, 0.520555, 0.505568, 0.491337}

Now I put this all together to get the 3D coordinates of the intersection points.

 klPts = Transpose[{kPts, lPts, MapThread[xStar1PtsF, {kPts, lPts}]}]
{{0., 0.637846, 0.607112}, {0.1, 0.629023, 0.636951}, {0.2, 0.616597,  0.664536}, 
 {0.3, 0.60183, 0.689845}, {0.4, 0.585783, 0.712876}, {0.5, 0.569174, 0.733694}, 
 {0.6, 0.552539, 0.752417}, {0.7, 0.536255, 0.769201}, {0.8, 0.520555, 0.784227}, 
 {0.9, 0.505568, 0.797682}, {1., 0.491337, 0.809751}}

And, finally, here is the plot of everything.

Show[
  Plot3D[{xStar1PtsF[k, l], xStar2PtsF[k, l]}, {k, 0, 1}, {l, 0, 1}, 
    AxesLabel -> {"k", "l", "x"^"*"}],
  Graphics3D[{Red, Sphere[#, Scaled[.015]] & /@ klPts}]]

plot

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.