Experiments on intersections of "bisector surfaces" for pairs of lines

in R^3. First we consider the two skew lines V(x,y) and V(y-1,z):

> with(linalg): with(Groebner):

Warning, new definition for norm

Warning, new definition for trace

> C1:=[x1-1,y1]; C2:= [y2-1,z2];

The conditions for (x,y,z) to be on the bisector surface are that

there are points (x1,y1,z1) in C1 and (x2,y2,z2) in C2 such

that (x-xi,y-yi,z-zi) are perpendicular to the tangent line to Ci

and ||(x-x1,y-y1,z-z1)|| = ||(x-x2,y-y2,z-z2)||:

> N1:=crossprod(grad(C1[1],[x1,y1,z1]),grad(C1[2],[x1,y1,z1]));

> N2:=crossprod(grad(C2[1],[x2,y2,z2]),grad(C2[2],[x2,y2,z2]));

> eqns1:=[op(C1),op(C2),dotprod(N1,[x-x1,y-y1,z-z1],orthogonal),dotprod(N2,[x-x2,y-y2,z-z2],orthogonal),(x-x1)^2+(y-y1)^2+(z-z1)^2-(x-x2)^2-(y-y2)^2-(z-z2)^2];

Use a lex Groebner basis to eliminate xi, yi, zi:

> G1:=gbasis(eqns1,plex(x1,y1,z1,x2,y2,z2,x,y,z));

The equation of the bisector surface:

> y=solve(G1[1],y);

This is a hyperbolic paraboloid (saddle surface)

> D1:=[x1+1,y1]; D2:= [y2+1,z2];

> M1:=crossprod(grad(D1[1],[x1,y1,z1]),grad(D1[2],[x1,y1,z1]));

> M2:=crossprod(grad(D2[1],[x2,y2,z2]),grad(D2[2],[x2,y2,z2]));

> eqns2:=[op(D1),op(D2),dotprod(M1,[x-x1,y-y1,z-z1],orthogonal),dotprod(M2,[x-x2,y-y2,z-z2],orthogonal),(x-x1)^2+(y-y1)^2+(z-z1)^2-(x-x2)^2-(y-y2)^2-(z-z2)^2];

> G2:=gbasis(eqns2,plex(x1,y1,z1,x2,y2,z2,x,y,z));

> Intersection:=[G1[1],G2[1]];

> gbasis(Intersection,plex(x,y,z));

"generic" picture -- intersection is a smooth curve of degree 4, genus 1 in 3-space; the

projection to a plane has two double points

> D1:=[2*x1-y1+z1+1,-x1+3*y1]; D2:= [-x2+y2-z2+2,z2];

> M1:=crossprod(grad(D1[1],[x1,y1,z1]),grad(D1[2],[x1,y1,z1]));

> M2:=crossprod(grad(D2[1],[x2,y2,z2]),grad(D2[2],[x2,y2,z2]));

> eqns2:=[op(D1),op(D2),dotprod(M1,[x-x1,y-y1,z-z1],orthogonal),dotprod(M2,[x-x2,y-y2,z-z2],orthogonal),(x-x1)^2+(y-y1)^2+(z-z1)^2-(x-x2)^2-(y-y2)^2-(z-z2)^2];

> G2:=gbasis(eqns2,plex(x1,y1,z1,x2,y2,z2,x,y,z));

> Intersectiongen:=[G1[1],G2[1]];

> GI:=gbasis(Intersectiongen,plex(x,y,z));

> GI[1];

> with(algcurves);

> evalf(singularities(GI[1],y,z));

> genus(GI[1],y,z);

>

> with(plots):

> implicitplot(GI[1],y=-50..50,z=-50..50,grid=[100,100]);

So the intersection of the bisector surfaces is a smooth quartic of genus 1 in P^3

Want to find a singular intersection case where the lines are in more "generic" position.

> g:='g';

> E1:=[x1-1,y1]; E2:= [a*x2+b*y2+c*z2+d,e*x2+f*y2+g*z2+h];

>

> Norm1:=crossprod(grad(E1[1],[x1,y1,z1]),grad(E1[2],[x1,y1,z1]));

> Norm2:=crossprod(grad(E2[1],[x2,y2,z2]),grad(E2[2],[x2,y2,z2]));

> eqns3:=[op(E1),op(E2),dotprod(Norm1,[x-x1,y-y1,z-z1],orthogonal),dotprod(Norm2,[x-x2,y-y2,z-z2],orthogonal),(x-x1)^2+(y-y1)^2+(z-z1)^2-(x-x2)^2-(y-y2)^2-(z-z2)^2];

> G3:=gbasis(eqns3,plex(x1,y1,z1,x2,y2,z2,x,y,z)):

> collect(G3[1],[x,y,z]);

Collecting the coefficients:

> G31F:=(b*g-c*f)^2*x^2-2*(a*g-c*e)*(b*g-c*f)*x*y+2*(b*g-c*f)*(a*f-b*e)*x*z + (a*g-e*c)^2*y^2-2*(a*f-b*e)*(a*g-c*e)*y*z-((b*g-c*f)^2+(a*g-e*c)^2)*z^2-2*((a*f-b*e)^2+(a*g-c*e)^2 + (b*g-c*f)^2)*x-(a*h-e*d)^2-(b*f-d*h)^2-(c*h-d*g)^2+(a*f-b*e)^2+(b*g-c*f)^2+(a*g-c*e)^2+(b*f-d*h)^2-(b*h-d*f)^2-2*((a*g-c*e)*(d*g-c*h)+(a*f-b*e)*(d*f-b*h))*x-2*((a*f-b*e)*(a*h-d*e)+(c*f-b*g)*(c*h-d*g))*y-2*((a*g-c*e)*(a*h-d*e)+(c*f-b*g)*(d*f-b*h))*z;

> expand(G31F-G3[1]);

Write coefficients of G3[1] in terms of the Plucker coordinates of the general line

> G31FD:=subs({a*g-e*c=Delta[13],b*g-c*f=Delta[23],a*f-b*e=Delta[12],a*h-e*d=Delta[14],b*h-d*f=Delta[24],c*h-d*g=Delta[34],d*g-c*h=-Delta[34],d*f-b*h=-Delta[24],-b*g+c*f=-Delta[23]},G31F);

> eqsFD:=[subs({x=1,y=1/2,z=0},G31FD),subs({x=1,y=1/2,z=0},diff(G31FD,x)),subs({x=1,y=1/2,z=0},diff(G31FD,z))];

> Plucker:=Delta[12]*Delta[34]-Delta[13]*Delta[24]+Delta[14]*Delta[23];

> expand(subs({Delta[13]=a*g-c*e,Delta[23]=b*g-c*f,Delta[12]=a*f-b*e,Delta[14]=a*h-d*e,Delta[24]=b*h-d*f,Delta[34]=c*h-d*g},Plucker));

> GDelta:=gbasis([op(eqsFD),Plucker],plex(Delta[34],Delta[24],Delta[14],Delta[23],Delta[13],Delta[12]));

First equation here says: over R we must have = 0 since

> subs(Delta[13]=0,GDelta);

Use original equations rather than the Groebner basis:

> subs(Delta[13]=0,[op(eqsFD),Plucker]);

From equations 2 and 3 here, or , and or

CASE A: If

> caseA:=subs({Delta[13]=0,Delta[24]=Delta[12]},[op(eqsFD),Plucker]);

> GA:=gbasis(caseA,plex(Delta[12],Delta[23],Delta[14],Delta[34]));

Factoring the first equation:

> expand((Delta[14]^2+Delta[34]^2)*(Delta[23]-Delta[34]));

This says that (or, since everything is real, = 0)

> subs(Delta[34]=Delta[23],GA);

<> 0 implies

> Dsubs:={Delta[13]=0,Delta[24]=Delta[12],Delta[34]=Delta[23],Delta[14]=-Delta[12]};

> G31FDs:=subs(Dsubs,G31FD);

> GInt:=[G1[1],G31FDs];

> GBInt:=gbasis(GInt,plex(x,y,z));

> expand(GBInt[1]/Delta[23]^2)/(-4);

Let . Then this has the form which factors

as: so the curve

is the union of two parabolas.

The other case:

> Dsubs2:={Delta[13]=0,Delta[14]=0,Delta[34]=0,Delta[24]=Delta[12]};

> subs(Dsubs2,[op(eqsFD),Plucker]);

> Int2:=[G1[1],subs(Dsubs2,G31FD)];

> GInt2:=gbasis(Int2,plex(x,y,z));

>

Same as previous.

CASE B: If and , the following shows we get no real solutions at all unless the lines are identical:

> subs({Delta[13]=0,Delta[12]=0,Delta[23]=0},[op(eqsFD),Plucker]);

CASE C: If and

> subs({Delta[13]=0,Delta[12]=0,Delta[24]=0},[op(eqsFD),Plucker]);

So either or . The second case takes us back to situation of B, and no real solutions. The first

case says either or . Either case says the lines are identical.

> subs({Delta[13]=0,Delta[12]=0,Delta[24]=0,Delta[34]=0,Delta[14]=0},[op(eqsFD),Plucker]);

CASE D: If and

> subs({Delta[13]=0,Delta[24]=Delta[12],Delta[23]=0},[op(eqsFD),Plucker]);

This is similar to the previous case C.

An example:

> teqs:=[subs({x=1,y=1/2,z=0},G3[1]),subs({x=1,y=1/2,z=0},diff(G3[1],x)),subs({x=1,y=1/2,z=0},diff(G3[1],z))];

> subvals:={b=-2,c=1,d=3,

> f=1,g=2,h=1};

> teqss:=subs(subvals,teqs);

> Gt:=gbasis(teqss,plex(a,e));

> fsolve(Gt[1],e,complex);

> subs(e=-2,Gt);

So a = -1, and as in CASE A (!)

> teqcheck:=subs({op(subvals),e=-2,a=-1},teqs);

> E2s:=subs({op(subvals),e=-2,a=-1},E2);

> Norm1:=crossprod(grad(E1[1],[x1,y1,z1]),grad(E1[2],[x1,y1,z1]));

> Norm2:=crossprod(grad(E2s[1],[x2,y2,z2]),grad(E2s[2],[x2,y2,z2]));

> eqns3:=[op(E1),op(E2s),dotprod(Norm1,[x-x1,y-y1,z-z1],orthogonal),dotprod(Norm2,[x-x2,y-y2,z-z2],orthogonal),(x-x1)^2+(y-y1)^2+(z-z1)^2-(x-x2)^2-(y-y2)^2-(z-z2)^2];

> G3:=gbasis(eqns3,plex(x1,y1,z1,x2,y2,z2,x,y,z));

> Intersection13:=[G1[1],G3[1]];

> GI:=gbasis(Intersection13,plex(x,y,z));

> with(algcurves):

> singularities(GI[1],y,z);

> genus(GI[1],y,z);

> implicitplot(GI[1],y=-10..10,z=-10..10,grid=[100,100]);

> with(plots):

> SP1:=implicitplot3d(G1[1],x=-5..5,y=-5..5,z=-5..5,grid=[30,30,30],style=PATCHNOGRID):

> SP2:=implicitplot3d(G3[1],x=-5..5,y=-5..5,z=-5..5,grid=[30,30,30]):

> display({SP1,SP2},axes=BOXED);

>