Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
Did you happen to notice that torus1 and torus2 are the same thing? –  J. M. May 8 '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 8 '13 at 15:46
    
@J.M. I am sorry, I have already revised. –  Shenghan Jiang May 8 '13 at 15:47
    
@bills use ParametricPlot3D –  Spawn1701D May 8 '13 at 15:50
2  
If you want to turn them into the form $f_i(x,y,z)=0$ try using Eliminate. –  Spawn1701D May 8 '13 at 15:55

3 Answers 3

up vote 8 down vote accepted

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

share|improve this answer
    
Thanks very much! I am satisfied with approximate solution. –  Shenghan Jiang May 8 '13 at 16:25

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 *)
share|improve this answer
    
Oxymoronizing an explicit implicitation –  belisarius May 8 '13 at 18:06
1  
@belisarius Yeh, it's something I learned how to do in Alchemistry. –  Daniel Lichtblau May 8 '13 at 18:07

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

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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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