5
$\begingroup$

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?

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

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

$\endgroup$
  • $\begingroup$ Thanks very much! I am satisfied with approximate solution. $\endgroup$ – Shenghan Jiang May 8 '13 at 16:25
9
$\begingroup$

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 *)
$\endgroup$
  • $\begingroup$ Oxymoronizing an explicit implicitation $\endgroup$ – Dr. belisarius May 8 '13 at 18:06
  • 1
    $\begingroup$ @belisarius Yeh, it's something I learned how to do in Alchemistry. $\endgroup$ – Daniel Lichtblau May 8 '13 at 18:07
4
$\begingroup$

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

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

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.