ParametricPlot3D Self Intersection

Question

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}}}]}]
 ]

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, -π, π}]

Answer 1

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]

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]]

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]

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]

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]

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]

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.

Answer 2

@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}]
]

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?