Solve for the intersection point of 3 surfaces

Question

Basically, I have three 2-torus (tilted) in 3-dimensional Euclidean space. They are expressed by parametric equations:

torus1[a_, b_] = {5/6 (7 + 5.5 Cos[a]) Cos[b] - (2 Sin[a])/3 - (7 + 5.5 Cos[a]) Sin[b]/6,
                  -(7 + 5.5 Cos[a]) Cos[b]/6 - (2 Sin[a])/3 + 5/6 (7 + 5.5 Cos[a]) Sin[b],
                  -(7 + 5.5 Cos[a]) Cos[b]/6 + (10 Sin[a])/3 - (7 + 5.5 Cos[a]) Sin[b]/6};
torus2[a_, b_] = {5/6 (11 + 2. Cos[a]) Cos[b] - Sin[a] - 1/6 (11 + 2. Cos[a]) Sin[b],
                  -1/6 (11 + 2. Cos[a]) Cos[b] - Sin[a] + 5/6 (11 + 2. Cos[a]) Sin[b],
                  -1/6 (11 + 2. Cos[a]) Cos[b] - Sin[a] - 1/6 (11 + 2. Cos[a]) Sin[b]};
torus3[a_, b_] = {5/6 (7 + (4 + 2.5 Cos[a]) Cos[b]) - 1.25 Sin[a] -
                  (4 + 2.5 Cos[a]) Sin[b]/6,
                  (-7 - (4 + 2.5 Cos[a]) Cos[b])/6 - 1.25 Sin[a] -
                  (4 + 2.5 Cos[a]) Sin[b]/6,
                  (-7 - (4 + 2.5 Cos[a]) Cos[b])/6 - 1.25 Sin[a] +
                  5/6 (4 + 2.5 Cos[a]) Sin[b]};

where a and b are some angular parameters run from $0$ to $2\pi$.

I want to solve for the coordinates of the intersection point of these three tori. However, since they are not described by the canonical form $f(x,y,z)=0$, I don't know how to deal with them. I have tried to transform the parametric equation to the canonical form; however, it is too complicated to do even using Mathematica (with GroebnerBasis).

How can I solve the equation torus1=torus2=torus3?

Did you happen to notice that torus1 and torus2 are the same thing? — J. M.'s missing motivation, May 08 '13 at 15:45
Is it a torus? When I try Plot3D[torus2[a, b], {a, 0, 2 Pi}, {b, 0, 2 Pi}] it looks like a bunch of intersecting sheets. — bill s, May 08 '13 at 15:46
If you want to turn them into the form $f_i(x,y,z)=0$ try using Eliminate. — Spawn1701D, May 08 '13 at 15:55

score 9 · Answer 1 · answered May 08 '13 at 17:05

I upvoted the response by @J.M. and was tempted to leave it at that. This is similar but automates the process a bit further by explicitly implicitizing (is that an oxymoron?) the tori. Somehow I think that step deserves mention since it can be a useful thing in its own right.

We start with code to take the trig parametrized tori and find algebraic implicit forms.

implicitize[tor_, vars_, newvars_] := Module[
  {c, s, cvars, svars, tvars, newtor, tpolys, allpolys},
  cvars = Map[c, vars];
  svars = Map[s, vars];
  tvars = Join[cvars, svars];
  newtor = TrigExpand[tor] /. {Cos[v_] :> c[v], Sin[v_] :> s[v]};
  tpolys = Map[c[#]^2 + s[#]^2 - 1 &, vars];
  allpolys = Join[Thread[newvars - newtor], tpolys];
  GroebnerBasis[allpolys, newvars, tvars]
  ]

We do this for the three tori.

implicittori = 
  Map[implicitize[Rationalize[#[a, b]], {a, b}, {x, y, z}] &,
   {torus1, torus2, torus3}];

Now solve for intersection points.

intersections = NSolve[implicittori == 0, {x, y, z}];
realpts = Select[intersections, FreeQ[#, Complex] &]

(* Out[26]= {{x -> 9.83321047045, y -> 0.115094126862, 
  z -> -2.79379942833}, {x -> 10.5867900819, y -> -2.73435297475, 
  z -> -2.34230565965}, {x -> 9.15662800421, y -> -0.48697702728, 
  z -> -0.858834384976}, {x -> 9.68512446602, y -> -3.01961259194, 
  z -> -0.256211546268}} *)

We can get exact solutions with modest further expenditure.

Timing[
 solns = Solve[implicittori == 0, {x, y, z}, Cubics -> False, 
    Quartics -> False];]

(* Out[28]= {7.390000, Null} *)

Timing[
 realptsexact = Select[solns, FreeQ[N[#, 25], Complex] &];]

(* Out[40]= {0.900000, Null} *)

realptsexact // Length

(* Out[41]= 4 *)

@belisarius Yeh, it's something I learned how to do in Alchemistry. — Daniel Lichtblau, May 08 '13 at 18:07

score 8 · Accepted Answer · answered May 08 '13 at 16:23

If you're happy with an approximate solution, you can use NSolve[]. As I mentioned in an answer to an earlier question of yours, GroebnerBasis[] can be used for parameter elimination. Let's do that for your three "tori":

t1 = First @ GroebnerBasis[Thread[{x, y, z} == Rationalize[torus1[a, b]]] ~Join~
                           {Cos[a]^2 + Sin[a]^2 == 1, Cos[b]^2 + Sin[b]^2 == 1},
                           {x, y, z}, {Cos[a], Sin[a], Cos[b], Sin[b]}];
t2 = First @ GroebnerBasis[Thread[{x, y, z} == Rationalize[torus2[a, b]]] ~Join~
                           {Cos[a]^2 + Sin[a]^2 == 1, Cos[b]^2 + Sin[b]^2 == 1},
                           {x, y, z}, {Cos[a], Sin[a], Cos[b], Sin[b]}];
t3 = First @ GroebnerBasis[Thread[{x, y, z} == Rationalize[torus3[a, b]]] ~Join~
                           {Cos[a]^2 + Sin[a]^2 == 1, Cos[b]^2 + Sin[b]^2 == 1},
                           {x, y, z}, {Cos[a], Sin[a], Cos[b], Sin[b]}];

(I used Rationalize[] to get rid of the inexact numbers, which don't play nice with GroebnerBasis[].)

Having obtained implicit Cartesian equations for your three surfaces, one can then do this:

NSolve[{t1 == 0, t2 == 0, t3 == 0}, {x, y, z}, Reals, WorkingPrecision -> 25]
   {{x -> 9.156628004210181380094443, y -> -0.4869770272804687489069949,
     z -> -0.8588343849755176760113928},
    {x -> 9.685124466015780195179393, y -> -3.019612591935055742262203,
     z -> -0.2562115462678606856207100},
    {x -> 9.833210470445689785623189, y -> 0.1150941268622902562908370,
     z -> -2.793799428332895281611490},
    {x -> 10.58679008185369737845052, y -> -2.734352974750579554184241,
     z -> -2.342305659647396259844311}}

One could use Solve[] instead to get exact solutions, but this is taking a long time on my box. (I'll edit this answer to include the results if and when it finishes.)

Thanks very much! I am satisfied with approximate solution. – Shenghan Jiang May 08 '13 at 16:25 — Shenghan Jiang, May 08 '13 at 16:25

andre314 · Answer 3 · 2013-05-08T17:15:30.760

4

You can minimize the distances between 3 points that each belong to a torus :

res = NMinimize[
        Norm[torus1[x1 , y1] - torus2[x2 , y2]] 
      + Norm[torus3[x3 , y3] - torus2[x2 , y2]],
      {x1, y1, x2, y2, x3, y3}]

(* torus1[x1 , y1] /. res[[2]]  -->  {10.5868, -2.73435, -2.34231}
   torus2[x2 , y2] /. res[[2]]  -->  {10.5868, -2.73435, -2.34231}
   torus3[x3 , y3] /. res[[2]]  -->  {10.5868, -2.73435, -2.34231} *)

The 3 torus and the point :

enter image description here

edited May 08 '13 at 17:15

answered May 08 '13 at 16:35

andre314

18,474
1
36
69

Smart method! However, in fact this equation has four solution. How to get another three? – Shenghan Jiang May 08 '13 at 16:51
not so simple, I think – andre314 May 08 '13 at 17:04

Solve for the intersection point of 3 surfaces

3 Answers3

Linked