# Solve for the intersection point of 3 surfaces

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 discontentment 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
• 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

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 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= {{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= {7.390000, Null} *)

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

(* Out= {0.900000, Null} *)

realptsexact // Length

(* Out= 4 *)

• Oxymoronizing an explicit implicitation – Dr. belisarius May 8 '13 at 18:06
• @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[]  -->  {10.5868, -2.73435, -2.34231}
torus2[x2 , y2] /. res[]  -->  {10.5868, -2.73435, -2.34231}
torus3[x3 , y3] /. res[]  -->  {10.5868, -2.73435, -2.34231} *)


The 3 torus and the point : • 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 – andre314 May 8 '13 at 17:04