##---------------------------------------------------------------------- ## Title: Intersection ## Created: March 2010 ## Authors: Luu Ba Thang ## Advisor: Laurent Buse ## ---------------------------------------------------------------------- ## Description: This part consists in some algorithms to compute the intersection points curve/curve in the plane, ## curve/curve in the space and curve/surface. ##------------------------------------------------------------------------- ## Version last modified 04/2010 ##------------------------------------------------------------------------- __ := NULL: addtohelp := proc() global __; __ := __, args: end: with(LinearAlgebra): with(MatrixPolynomialAlgebra): with(PolynomialTools): with(Groebner): #### Intersection of curve/curve and curve/surface in the space ######## #---------------------------------------------------------------------- addtohelp(CurveSurfaceHomogeneous): #---------------------------------------------------------------------- ## HELP CurveSurfaceHomogeneous ## CALLING SEQUENCE: CurveSurfaceHomogeneous(F,var1,G,var2) ## DESCRIPTION: ## - Given a parameterized homogeneous curve F:=(f0,f1,f2,f3) with variable var1 ## - Given a parameterized homogeneous surface G:=(g0,g1,g2,g3) with variable var2 ## - d: degree of parameterized G ## - L: List of implicit variables correspond to G ## RETURN: Locus points of curve and surface correspond to var1[1]:=1;. ## EXAMPLES: ## Curve C with parameterization F:=[a,0,0,b] and Surface S with homogeneous parameterization ## G:=[s^2+t^2+u^2,2*s*u,2*s*t, s^2-t^2-u^2], var1:=[a,b], var2:=[s,t,u]; L:=[x,y,z,w], ## CurveSurfaceHomogeneous(F,var1,G,var2) ## {[1, 0, 0, -1.], [1, 0, 0, 1.]} CurveSurfaceHomogeneous:=proc(F,var1,G,var2) local F1,L1,Mt,i,j,m,n,P,d, varimpl; m:=nops(F); d:=0; for i from 1 to m do d:= max(d,degree(G[i],var2)); end do; F1:=subs(var1[1]=1,F); varimpl:=[x,y,z,w]; Mt:=MatrixRepSurface(G,var2,varimpl,2*d-1); for i from 1 to m do Mt:=subs(varimpl[i]=F1[i],Mt); end do; L1:= evalf(ComputationEigenvalues(Mt,var1[2])); n:=nops(L1); P:=[]; for j from 1 to n do F1:=subs(var1[1]=1,F); for i from 1 to m do F1[i]:=subs(var1[2]=L1[j],F1[i]); end do; P:=[op(P),F1]; end do; P:={op(P)}; RETURN(P); end proc: #---------------------------------------------------------------------- addtohelp(SpaceCurveIntersectionHomogeneous): #---------------------------------------------------------------------- ## HELP SpaceCurveIntersectionHomogeneous ## CALLING SEQUENCE: SpaceCurveIntersectionHomogeneous(F,var1,G,var2) ## DESCRIPTION: ## - Given a homogeneous parameterized curve F:=[f0,f1,f2,f3] with variable var1 ## - Given a homogeneous parameterized curve G:=[g0,g1,g2,g3] with variable var2 ## RETURN: Locus points of curve and curve. ## EXAMPLES: ## Curve C1 with homogeneous parameterization F:= [a^3, a^2*b, a*b^2, b^3]; ## Curve C2 with homogeneous parameterization G:=[s^4, s^3*t, s^2*t^2, t^4] ## SpaceCurveIntersectionHomogeneous([a^3, a^2*b, a*b^2, b^3], [a,b], [s^4, s^3*t, s^2*t^2, t^4], [s, t]) ## {[1, 0., 0., 0.], [1, 1., 1., 1.]} SpaceCurveIntersectionHomogeneous:=proc(F,var1, G,var2) local F1,L1,Mt,i,j,m,n,P,varimpl; F1:=subs(var1[1]=1,F); varimpl:=[x[1],x[2],x[3],x[4]]; Mt:=MatrixSpaceCurveBase(G,var2,varimpl); m:=nops(F); for i from 1 to m do Mt:=subs(varimpl[i]=F1[i],Mt); end do; L1:= evalf(ComputationEigenvalues(Mt,var1[2])); n:=nops(L1); P:=[]; for j from 1 to n do F1:=subs(var1[1]=1,F); for i from 1 to m do F1[i]:=subs(var1[2]=L1[j],F1[i]); end do; P:=[op(P),F1]; end do; P:={op(P)}; RETURN(P); end proc: ## HELP SpaceCurveIntersection ## CALLING SEQUENCE: SpaceCurveIntersection(f,var1,g,var2) ## DESCRIPTION: ## - Given a parameterized curve F:=[f0,f1,f2,1] with a variable var1 ## - Given a parameterized curve G:=[g0,g1,g2,1] with a variable var2 ## RETURN: Locus points of curve and curve. ## EXAMPLES: ## Curve C1 with the parameterization f:= [a^3, a^2, a, 1]; ## Curve C2 with the parameterization G:=[s^4, s^3, s^2,1] ## SpaceCurveIntersection([a^3, a^2, a,1], a , [s^4, s^3, s^2, 1], s) ## {[0., 0., 0., 1], [1., 1., 1., 1]} #---------------------------------------------------------------------- addtohelp(SpaceCurveIntersection): #---------------------------------------------------------------------- SpaceCurveIntersection:=proc(f,var1, g, var2) local F1,L1,Mt,i,j,m,n,P,varimpl,varhomo,y,g1,g2; F1:=f; m:=nops(g); varhomo:=y; varimpl:= [seq(x[i], i = 1 .. m)]; g1:=g; g2:=Multiply(Matrix(g),Vector(varimpl)); g2[1]:= Homogenize(g2[1], varhomo); for j from 1 to m do g1[j]:=coeff(g2[1],varimpl[j]); end do; Mt:= MatrixSpaceCurveBase(g1,[varhomo,var2],varimpl); for i from 1 to m do Mt:=subs(varimpl[i]=F1[i],Mt); end do; L1:= evalf(ComputationEigenvalues(Mt,var1)); n:=nops(L1); P:=[]; for j from 1 to n do F1:=f; for i from 1 to m do F1[i]:=subs(var1=L1[j],F1[i]); end do; P:=[op(P),F1]; end do; P:={op(P)}; RETURN(P); end proc: ############### Curve/Curve in the plane ################ #---------------------------------------------------------------------- addtohelp(PlaneCurveIntersectionHomogeneous): #---------------------------------------------------------------------- ## HELP PlaneCurveIntersectionHomogeneous ## CALLING SEQUENCE: PlaneCurveIntersectionHomogeneous(F,var1,G,var2) ## DESCRIPTION: ## - Given a homogeneous parameterized curve F:=[f0,f1,f2] with variable var1 ## - Given a homogeneous parameterized curve G:=[g0,g1,g2] with variable var2 ## RETURN: Locus points of curve and curve. ## EXAMPLES: ## Curve C1 with homogeneous parameterization F:= [-4*s^3+(1/2)*t^3, s*t^2+s^2*t, t^3]; ## Curve C2 with homogeneous parameterization G:=[u^3, u^2*v, v^3] ## PlaneCurveIntersectionHomogeneous([-4*s^3+(1/2)*t^3, s*t^2+s^2*t, t^3], [s,t], [u^3, u^2*v, v^3], [u, v]) ## [-4., 0., 0.], [-3.197948630, 2.540924237, 1.604102740], [7.54776894, 10.95775934, 23.09553787], ## [-1.317686880-1.089747447*I, -1.938625677+1.683795576*I, 5.364626241-2.179494894*I], ## [-1.317686880+1.089747447*I, -1.938625677-1.683795576*I, 5.364626241+2.179494894*I], ## [2.142776730-2.395430436*I, -4.810716109-1.824961439*I, 12.28555346-4.790860873*I], ## [2.142776730+2.395430436*I, -4.810716109+1.824961439*I, 12.28555346+4.790860873*I] PlaneCurveIntersectionHomogeneous:=proc(F,var1, G,var2) local F1,L1,Mt,i,j,m,n,P,varimpl,d; F1:=subs(var1[1]=1,F); d:= Degree(Matrix(G),var2[2]); print(d); varimpl:=[x[1],x[2],x[3]]; Mt:=MatrixPlaneCurve(G,var2,varimpl,2*d-1); m:=nops(F); for i from 1 to m do Mt:=subs(varimpl[i]=F1[i],Mt); end do; L1:= evalf(ComputationEigenvalues(Mt,var1[2])); n:=nops(L1); P:=[]; for j from 1 to n do F1:=subs(var1[1]=1,F); for i from 1 to m do F1[i]:=subs(var1[2]=L1[j],F1[i]); end do; P:=[op(P),F1]; end do; P:={op(P)}; RETURN(P); end proc: #---------------------------------------------------------------------- addtohelp(PlaneCurveIntersection): #---------------------------------------------------------------------- ## HELP PlaneCurveIntersection ## CALLING SEQUENCE: PlaneCurveIntersection(f,var1,g,var2) ## DESCRIPTION: ## - Given a homogeneous parameterized curve F:=[f0,f1,f2] with a variable var1 ## - Given a homogeneous parameterized curve G:=[g0,g1,g2] with a variable var2 ## RETURN: Locus points of curve and curve. ## EXAMPLES: ## Curve C1 with a parameterization f := [-4+(1/2)*t^3, t^2+t, t^3]; ## Curve C2 with a parameterization g := [1, v, v^3] ## PlaneCurveIntersection(f, t, g, v) ## [-4., 0., 0.], [-3.197948630, 2.540924237, 1.604102740], [7.54776894, 10.95775934, 23.09553787], ## [-1.317686880-1.089747447*I, -1.938625677+1.683795576*I, 5.364626241-2.179494894*I], ## [-1.317686880+1.089747447*I, -1.938625677-1.683795576*I, 5.364626241+2.179494894*I], ## [2.142776730-2.395430436*I, -4.810716109-1.824961439*I, 12.28555346-4.790860873*I], ## [2.142776730+2.395430436*I, -4.810716109+1.824961439*I, 12.28555346+4.790860873*I] PlaneCurveIntersection:=proc(f,var1, g, var2) local F1,L1,Mt,i,j,m,n,P,varimpl,varhomo,y,g1,g2,deg; F1:=f; m:=nops(g); varhomo:=y; varimpl:= [seq(x[i], i = 1 .. m)]; g1:=g; deg:= Degree(Matrix(g),var2); g2:=Multiply(Matrix(g),Vector(varimpl)); g2[1]:= Homogenize(g2[1], varhomo); for j from 1 to m do g1[j]:=coeff(g2[1],varimpl[j]); end do; Mt:= MatrixPlaneCurve(g1,[varhomo,var2],varimpl,2*deg-1); for i from 1 to m do Mt:=subs(varimpl[i]=F1[i],Mt); end do; L1:= evalf(ComputationEigenvalues(Mt,var1)); n:=nops(L1); P:=[]; for j from 1 to n do F1:=f; for i from 1 to m do F1[i]:=subs(var1=L1[j],F1[i]); end do; P:=[op(P),F1]; end do; P:={op(P)}; RETURN(P); end proc: [__];