11
$\begingroup$

I would like to know where my ParametricPlot3D intersects itself and to extract all intersections (points, lines, surfaces...) to animate them by changing P:

ParametricPlot3D[{fX[u,v,P], fY[u,v,P], fZ[u,v,P]}, (* P is a vector of parameters *)
 {u,0,1}, {v,0,1}]

Example: The vector P is not important, so here is an example from the documentation for ParametricPlot3D. The problem is to find the set of points where the surface intersects itself. In this simple case, the self-intersection seems to consist of two line segments, which I was able simply to guess by looking. What I seek is a more general method.

Show[
 ParametricPlot3D[{Cos[u], Sin[u] + Cos[v], Sin[v]},
  {u, 0, 2 π}, {v, -π, π}, PlotStyle -> Opacity[0.5], 
  ColorFunction -> (ColorData["Rainbow"][#5] &)],
 Graphics3D[{Thick, Red, 
   Line[{{{1, 0, 1}, {-1, 0, -1}}, {{1, 0, -1}, {-1, 0, 1}}}]}]
 ]

Mathematica graphics

A small perturbation might serve as a working example:

ParametricPlot3D[
 {Cos[u + 0.2 Sin[v]], Sin[u] + Cos[v + 0.2 Sin[u]], Sin[v + 0.2 Cos[u]]},
 {u, 0, 2 π}, {v, -π, π}]
$\endgroup$
  • $\begingroup$ This is a duplicate of many duplicates, e.g. (126847) and (127161). $\endgroup$ – corey979 Nov 5 '16 at 16:57
  • $\begingroup$ Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 Nov 5 '16 at 17:05
  • 1
    $\begingroup$ That seem a duplicate due to the lack of information from my self: - fX, fY, fZ are from Interpolation Each time I use FindRoot, Reduce, NSolve Mathematica are not happy. "Not valid variable, Reduce was unable to solve the system with inexact coefficients..." Thanks PS: All answer I saw are on 2D that's why I post it again. $\endgroup$ – chkone Nov 5 '16 at 20:29
  • 1
    $\begingroup$ Having a simple working example will make it easier and more likely to get help. $\endgroup$ – Michael E2 Nov 6 '16 at 11:14
  • 1
    $\begingroup$ @corey979 I think the OP is right that obtaining the curves where a parameterized surface intersects itself is significantly distinct from getting the *points where a parameterized curve intersects itself (as in the current duplicate). Would you agree? $\endgroup$ – Michael E2 Nov 6 '16 at 11:35
6
$\begingroup$

This 3D problem can be solved based on two 2D solutions, 76737 and 33947, by Michael E2. For the second parametric system in the question,

plt = ParametricPlot3D[{Cos[u + 0.2 Sin[v]], Sin[u] + Cos[v + 0.2 Sin[u]], 
    Sin[v + 0.2 Cos[u]]}, {u, 0, 2 π}, {v, -π, π}, PlotStyle -> Opacity[0.5], 
    ColorFunction -> (ColorData["Rainbow"][#5] &), ImageSize -> Large]

enter image description here

generate slices for various values of z, for instance,

With[{z0 = .5}, RegionPlot3D[DiscretizeGraphics[plt], MeshFunctions -> {#3 &}, 
    Mesh -> {{z0}}, MeshStyle -> {Directive[Thick, Blue]}, PlotStyle -> None] /. 
    Graphics3D[g_, opts___] :> Graphics[g /. {x_Real, y_Real, Real} :> {x, y}, 
    FilterRules[{opts}, Graphics], Frame -> True]]

enter image description here

The points comprising the two curves then can be extracted from the plot, converted into InterpolatingFunctions, and FindRoot used to find the intersections. To provide initial guesses, examine the plots above, provide an expression for the guesses, and then use Nearest to determine where those nearest points are on the two curves. Finally, loop over slices in z to find the two lines of intersection, as follows:

pltdg = DiscretizeGraphics[plt]; tbl = {};
Do[guess = {sgn, 0}; Do[rp = RegionPlot3D[pltdg, MeshFunctions -> {#3 &}, 
    Mesh -> {{i}}, MeshStyle -> Red, PlotStyle -> None] /. Graphics3D[g_, opts___] :> 
    Graphics[g /. {x_Real, y_Real, z_Real} :> {x, y}, 
    FilterRules[{opts}, Graphics], Frame -> True];
  pts = Cases[((rp // InputForm)[[1, 1]]) // Normal, Line[z_] -> z, Infinity];
  i1x = Interpolation[pts[[1, All, 1]]];
  i1y = Interpolation[pts[[1, All, 2]]];
  i2x = Interpolation[pts[[2, All, 1]]];
  i2y = Interpolation[pts[[2, All, 2]]];
  s = Quiet@FindRoot[{i1x[t1] == i2x[t2], i1y[t1] == i2y[t2]}, 
    {t1, First@Nearest[First@pts -> Automatic, guess]}, 
    {t2, First@Nearest[Last@pts -> Automatic, guess]}];
  AppendTo[tbl, {i1x[t1], i1y[t1], i} /. s]; 
  guess = {sgn Abs[i1x[t1]], i1y[t1]} /. s, 
  {i, -.96, .96, .01}], 
{sgn, -1, 1, 2}]
Show[plt, ListPointPlot3D[tbl, BoxRatios -> {1, 2, 1}, 
  PlotRange -> {{-1, 1}, {-2, 2}, {-1, 1}}, PlotStyle -> Red], ImageSize -> Large]

enter image description here

Addendum: Largely Symbolic Solution

Because Sin[v + Cos[u]/5] == z0 can be inverted symbolically to obtain v as a function of u for given z0, a faster and more accurate solution exists than the general solution given above. First, obtain expressions for x and y as functions of u for fixed z0.

Clear[v]
Reduce[Sin[v + Cos[u]/5] == z0, v, Reals];
Flatten@Solve[Simplify[%, C[1] ∈ Integers && -1 <= z0 <= 1], v];
slice = FullSimplify[{Cos[u + Sin[v]/5], Sin[u] + Cos[v + Sin[u]/5]} /. # & /@ %, 
    C[1] ∈ Integers]
(* {{Cos[u + 1/5 Sin[ArcSin[z0] + Cos[u]/5]], 
     -Cos[1/5 (5 ArcSin[z0] + Cos[u] - Sin[u])] + Sin[u]}, 
    {Cos[u + 1/5 Cos[ArcCos[z0] + Cos[u]/5]], 
      Cos[1/5 (5 ArcSin[z0] - Cos[u] + Sin[u])] + Sin[u]}} *)

This expression gives the 2D slice at z0 through the first plot above. For instance, for z0 == 0.5,

Show[ParametricPlot[# /. z0 -> 0.5, {u, 0, 2 π}, PlotRange -> {{-1, 1}, {-2, 2}}] &
    /@ slice, PlotRangeClipping -> False]

enter image description here

This plot is equivalent to the second plot above but has higher resolution. A quick investigation shows that the two curves do not intersect in the vicinity of z0 == 0.2, as shown in this blow-up plot.

Show[ParametricPlot[# /. z0 -> 0.2, {u, 0, 2 π}, 
    PlotRange -> {{-.1, .1}, {-.01, .01}}] & /@ slice, AspectRatio -> 1]

enter image description here

The limits of the range of z0 for which no intersections exist are given by

zp = z0 /. FindRoot[Thread[slice[[1]] == {0, 0}], {z0, .23}, {u, Pi/2}]
(* 0.236925 *)
zm = z0 /. FindRoot[Thread[slice[[1]] == {0, 0}], {z0, .16}, {u, Pi/2}]
(* 0.159385 *)

With this information the 3D intersection curve can be obtained.

eq = Equal @@@ Transpose@MapThread[#1 /. u -> #2 &, {slice, {u, um}}]
(* { Cos[u + 1/5 Sin[ArcSin[z0] + Cos[u]/5]] == 
       Cos[um + 1/5 Cos[ArcCos[z0] + Cos[um]/5]], 
    -Cos[1/5 (5 ArcSin[z0] + Cos[u] - Sin[u])] + Sin[u] == 
       Cos[1/5 (5 ArcSin[z0] - Cos[um] + Sin[um])] + Sin[um]} *)

tbl = {}; guess = .2;
Do[If[i < zm || i > zp, s = Quiet@FindRoot[eq /. z0 -> i, {u, guess}, {um, guess}];
    AppendTo[tbl, Join[slice[[1]] /. z0 -> i /. s, {i}]]; 
    guess = Abs[slice[[1, 1]] /. z0 -> i /. s]], {i, -99/100, 99/100, 1/1000}]
Join[tbl, tbl /. {z1_, z2_, z3_} -> {-z1, -z2, z3}];
Show[plt, ListPointPlot3D[%, BoxRatios -> {1, 2, 1}, 
    PlotRange -> {{-1, 1}, {-2, 2}, {-1, 1}}, PlotStyle -> Red], ImageSize -> Large]

enter image description here

It agrees well with the third plot above except near z0 == 0.2, where the newly computed curve is the more accurate. Computing it also is about two orders of magnitude faster on a per-point basis.

$\endgroup$
  • $\begingroup$ This is a nice approach but it does not really find the intersection curve, only a finite number of points on the curve with no connectivity information between them. Also, if I change the surface to {Cos[u + 0.5 Sin[v]], Sin[u] + Cos[v + 0.5 Sin[u]], Sin[v + 0.5 Cos[u]]} it gives some strange results. $\endgroup$ – Rahul Nov 11 '16 at 23:20
  • $\begingroup$ @Rahul The points accurately represent the curves, and they can be connected in the plot, if desired. With respect to the surface in your comment, it is necessary to tune guess for each surface. $\endgroup$ – bbgodfrey Nov 11 '16 at 23:27
  • $\begingroup$ I see, it requires some manual intervention and some foreknowledge of where the intersections lie. $\endgroup$ – Rahul Nov 12 '16 at 22:46
  • $\begingroup$ @Rahul The foreknowledge is obtained from the second figure, which shows roughly where the intersections occur. $\endgroup$ – bbgodfrey Nov 12 '16 at 23:12
  • $\begingroup$ i think as a general approach you need to find the points then have a scheme to order/connect them. (the connection is perhaps worthy of a new question, but do a search first). $\endgroup$ – george2079 Nov 14 '16 at 13:10
3
+50
$\begingroup$

@bbgodfrey's answer got me thinking that this could be thought of as tracking roots as a function of a parameter (here u), so I might be able to use a function I put together for that purpose here. A self-intersection will be where {x,y,z}=={x2,y2,z2} for unique values of {u,v} and {u2,v2} (modulo 2π in this case). That gives us three equations and four unknowns {u,v,u2,v2}, so treat one of the unknowns as a parameter (I'll use u), then solve for {v,u2,v2} as a function of u.

First, define TrackRootPAL from here.

Next, define two versions of the surface:

x := Cos[u + 0.2 Sin[v]];
y := Sin[u] + Cos[v + 0.2 Sin[u]];
z := Sin[v + 0.2 Cos[u]];

x2 := Cos[u2 + 0.2 Sin[v2]];
y2 := Sin[u2] + Cos[v2 + 0.2 Sin[u2]];
z2 := Sin[v2 + 0.2 Cos[u2]];

Find a starting point somewhere on the self-intersection using FindRoot:

u = 0.5;
init = FindRoot[{x == x2, y == y2, z == z2}, {v, 0.1}, {u2, π}, {v2, π}]
(* {v -> 8.24744, u2 -> 5.44587, v2 -> -5.41526} *)

Then track it with TrackRootPAL, which returns a set of InterpolatingFunctions:

Clear[u];
tr = TrackRootPAL[{x - x2, y - y2, z - z2}, {v, u2, v2}, {u, 0, 2 π}, 0.5, {v, u2, v2} /. init]

Thanks to @bbgodfrey's analysis, we know there is another set of self-intersections. Repeat the above with a different initial point to find it:

u = 4;
init = FindRoot[{x == x2, y == y2, z == z2}, {v, 0.1}, {u2, π}, {v2, -π}];
Clear[u];
tr2 = TrackRootPAL[{x - x2, y - y2, z - z2}, {v, u2, v2}, {u, 0, 2 π}, 4, {v, u2, v2} /. init, NDSolveOpts -> {AccuracyGoal -> 6}]

Finally plot it all together:

Show[
 ParametricPlot3D[{x, y, z}, {u, 0, 2 π}, {v, -π, π}, PlotStyle -> Opacity[0.5]],
 ParametricPlot3D[{x2, y2, z2} /. {u2 -> u2[u], v2 -> v2[u]} /. tr, {u, 0, 2 π}, PlotStyle -> {Red, Thick}],
 ParametricPlot3D[{x2, y2, z2} /. {u2 -> u2[u], v2 -> v2[u]} /. tr2, {u, 0, 2 π}, PlotStyle -> {Red, Thick}]
]

Mathematica graphics

Of course this makes it look easier than it actually was. The trickiest parts were 1) knowing how many disjunct sets of self-intersections there are and 2) finding good initial guesses to locate them. @bbgodfrey's answer took care of 1) for me; 2) required some trial and error. Maybe this could be automated somehow?

$\endgroup$
  • $\begingroup$ Elegant solution (+1). Does your TrackRootPAL exist as a package? $\endgroup$ – bbgodfrey Nov 14 '16 at 14:00
  • $\begingroup$ @bbgodfrey No, I'm not confident that it's robust and well-formulated enough to package up. It'd be great if more experienced folks gave it a careful vetting! $\endgroup$ – Chris K Nov 14 '16 at 14:16

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.