Real Plane Curves - a numerical approach via Mathematica

Barry H Dayton

Oakton Community College
Northeastern Illinois University

Introduction

In classical curve theory all coefficients and points are exact, often rational numbers.  In this case most points can not be explicitly described, thus methods avoid working computationally with points.  In this talk we work over the field of Real numbers and describe an alternate approach using explicit computations using Mathematica Machine Numbers as coefficients.

In this talk we will describe procedures to construct curves,  identify special points, transform curves, decompose curves and plot and sketch curves.  Our main goal is to understand the topology of the curve.

By a curve we mean the set of all real zeros of a polynomial in several variables, plus infinite points.  Although we often begin with a curve in integer coefficients so that work can be replicated, our software  often converts the curve to machine coefficients and most output is in that form.

There are  special techniques in working numerically rather than exactly.  Some of our algorithms are non-deterministic and probabistic, they may not work at every execution, however when they do work they generally produce a correct answer.  Because of ranomization they may not give exactly the same answer each run.  However the randomization allows us to concentrate on generic curves rather than spending  time on special cases.

Only a few of the procedures are implemented in this worksheet.  For a full current implementation see  my website MouseOver and download  RealPlaneCurves.zip.

Constructing Curves

Random Curves

(* Total degree of multivariate polynomial  *)
tdeg[f_,X_]:= With[{T= Thread[X->Table[t_,{i,Length[X]}]]}, 
Max[Exponent[MonomialList[f]/.T,t_]]
];
(*  monomials of order exactly n in variables  X  *)
tdv[X_,n_]:=tdv[X,n] =
 If[Length[X]==1,{X[[1]]^n},If[n == 1,X,Join[X[[1]]*tdv[X,n-1],tdv[Drop[X,1],n]]]];
(* Monomials of order n or less in variables X  *)
ctdv[X_,n_]:=Prepend[Flatten[Table[tdv[X,j],{j,n}]],1];
randomintegercurve[X_,d_,k_]:=Module[{f,M,L,kk,b},
f=0;I
While[tdeg[f,X]<d,
b=Binomial[d+2,2];
If[k<2,kk=b; L=ctdv[X,d], kk=Min[k,b];
M=Drop[ctdv[{x,y},d],1];
L=Append[RandomSample[M,kk-1],1]];
f=RandomInteger[{-5,5},kk].L];
f
]
randomrealcurve[X_,d_]:=RandomReal[{-1,1},Binomial[d+2,2]].ctdv[X,d];
line[p_,q_,X_]:=Module[{e1,e2,a,b,c,M},
e1=Switch[Length[p],2,Append[p,1],3,p];
e2=Switch[Length[q],2,Append[q,1],3,q];
M={e1,e2,RandomReal[{-1,1},3]};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).{x,y,1}
]
f=randomintegercurve[{x,y},4,6];
g=randomrealcurve[{x,y},4];
line[p_,q_,X_]:=Module[{e1,e2,a,b,c,M},
e1=Switch[Length[p],2,Append[p,1],3,p];
e2=Switch[Length[q],2,Append[q,1],3,q];
M={e1,e2,RandomReal[{-1,1},3]};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).{x,y,1}
];
NewtonHyperbola[n_]:=Module[{v,c,k,n2dig,s,maxn},
maxn=32767; (*set maxn=1023 for cubics 32767 for quartics, 2097151 for quintics *)
If[(n> maxn) ||(n<0), Print["n out of range"];Return[Fail]];
v={{0,0},{1,0},{0,1},{2,0},{1,1},{0,2},{3,0},{2,1},{1,2},{0,3},{4,0},{3,1},{2,2},{1,3},{0,4},{5,0},{4,1},{3,2},{2,3},{1,4},{0,5}};
c={1,1,1,.2,1,.2,.008,.2,.2,.008,.000064,
.008,.2,.008,.000064,.2^10,.2^6,.008,.008,.2^6,.2^10};
k=Which[ n<8,3,n<64,6,n<1024 ,10,n<32768,15,True,21];
n2dig=IntegerDigits[n,2,k];
s=Reverse[n2dig];
Total[Table[(-1)^s[[i]]c[[i]]x^v[[i,1]] y^v[[i,2]],{i,k}]]
];

f=randomintegercurve[{x,y},4,7]

The last argument allows one to select the number of non-zero terms.  A significant number of these curves are pathological.  The random real curves are very generic.

g=randomrealcurve[{x,y},4]

Our Newton Hyperbolas are motivated by the study of cubics by Newton  and higher degree curves by the contemporary mathematician Oleg Viro.  These are not well scaled and can be difficult, and so make good examples for testing the algorithms.

h=NewtonHyperbola[425]

A very general line making routine is useful.

l = line[{1, 2}, RandomReal[{-3, 3}, 2], {x, y}]

Affine Points

tangentrealpoints[f_,kk_,X_,Q_:True]:=Module[{J, sol,rp, p,q,k},
If[(kk/.Thread[X->RandomReal[{-1,1},2]])==0, k=RandomReal[{-1,1},2].X, k=kk];
If[Q,Print["k= ",k]]; 
J=Det[D[{f,k},{X}]];
sol=NSolve[{f,J}] ;
rp=Cases[X/.sol,{_Real,_Real}];
If[rp=={}&& Q,Print["trp: no point found"];Return[]];
rp
];
infiniterealpoints[f_,X_]:=Module[{F,pp,sol,RR,iptr,hr, hp, ml,rp,k, degf,ips},
ml=MonomialList[f];
degf=tdeg[f,X];
hp=Total[Select[ml,tdeg[#,X]== degf&]];
k=RandomReal[{-1,1},2].X-1;
sol=NSolve[{hp,k}];
rp=Cases[X/.sol,{_Real,_Real}];
If[rp=={},Print["no infinite points"];Return[{}],
ips=rp.{{1,0,0},{0,1,0}}];
Normalize/@DeleteDuplicates[ips]
];
randomrealpoint[f_,X_]:=Module[{rp,ip,q,k,sol,p},
(* First Try -- ovals*)
rp=tangentrealpoints[f,0,X,False];
If[Length[rp]>0, q=rp[[RandomInteger[{1,Length[rp]}]]],
(* second attempt no ovals*)
rp=infiniterealpoints[f,X];
If[Length[rp]>0, ip=rp[[RandomInteger[{1,Length[rp]}]]], Return[Fail]];
k=(ip+RandomReal[{-.2,.2},3]).Append[X,1];
sol=NSolve[{f,k}];
rp=Cases[X/.sol,{_Real,_Real}];
If[Length[rp]>0, q=rp[[RandomInteger[{1,Length[rp]}]]], Return[Fail]]];
If[Norm[D[f,{X}]/.Thread[X->q]]<1.*^-8, Print["Singular Point"]];
q
];

Finding real points may not be so easy, the real locus of a real curve may be small and far away from the origin.  The locus may even be finite.  Our main trick is to find points with  tangent lines  parallel to a linear form k = ax+by.  Here is a simplified version of the code which was motivated by conversations with Daniel Lichtblau.

trp[f_,k_,X_]:=Module[{J, sol},
J=Det[D[{f,k},{X}]];  (*Jacobian determinant *)
sol=NSolve[{f,J}] ;
Cases[X/.sol,{_Real,_Real}]]

trp[f,x+y,{x,y}]

With our general routine we can ask for a random line by setting k=0.

tangentrealpoints[f,0,{x,y}]

Singular points are generally tangent to any line, they appear with multiplicity greater than 1:

tangentrealpoints[x^3-x^2-y^2,0,{x,y}]

"k= "0.4885514062129057` x+0.7815200264731335` y

{{-3.9652568315790904`*^-17,-6.34309426721409`*^-17},{0.`,0.`}}

This point is actually an isolated point.

randomrealpoint  has additional tricks to find points in most cases tangentrealpoints fails.  If randomrealpoint fails repeatedly the curve has no real points.

randomrealpoint[f,{x,y}]

Infinite points and Asymptotes

An infinite point is essentially a slope, specifically the slope of an asymptote.  Although we are working in the affine plane we use notation from projective geometry, the infinite point {1,m,0} is the point at the end of a line of slope m.  Infinite points are homogeneous, that is {a,b,0}={ra,rb,0} for any real r ≠  0.

infiniterealpoints[h,{x,y}]

To find the asymptote through one of the infinite points we have a procedure  Asymptotes[f,p,{x,y}]not implemented in this worksheet.  For example

c=x^2-y^2-1;
infpts=infiniterealpoints[c,{x,y}]
Asymptotes[c,infpts[[1]],{x,y}]

{{0.7071067811865474`,-0.7071067811865478`,0.`},{-0.7071067811865476`,-0.7071067811865475`,0.`}}

{-0.7071067811865481` x-0.7071067811865476` y}

For some good examples of this note that rational functions y = ((p(x))/(q(x))) give algebraic curves by f = yq(x)-p(x).  (Caution:  p(x), q(x) should not have  common zeros.)  If p(x)= (ax^2)+bx+c, q(x)=dx+e then one infinite point is {a,d,0} by inspection.  But the asymptote to this infinite point depends on the other coefficients.

Asymptotes[y (3x-1)- (2x^2+x-2),{3,2,0},{x,y}]

{-0.419590679148345`-0.5035088149780139` x+0.7552632224670197` y}

Note that one infinite point may support several asymptotes, for example the infinite point for any vertical asymptote is {0,1,0} but a rational function may have several vertical asymptotes.

Asymptotes[y x(x^2-1)-1,{0,1,0},{x,y}]

{-0.7071067811865474`-0.7071067811865475` x,0.7071067811865479` -0.7071067811865471` x,0.9999999999999998` x}

Note that for the line function either point on the line may be an infinite point, that is you can specify point and slope.

line[{2,3},{1,-2,0},{x,y}]

Fractional Linear Transformations

spec2=Take[#,2]/#[[3]]& ;
flt[p_,A_]:=spec2[A.Append[p,1]];
flti[p_,A_]:=Module[{q, o},
If[Length[p]==2,
q=A.Append[p,1];
If[Abs[q[[3]]]>1.*^-10,q=spec2[q]],
q=A.p;
If[Abs[q[[3]]]>1.*^-10,q=spec2[q]]];
q];

In keeping with our affine-with-infinite-point orientation the most useful transformation group is the group of  fractional linear transformations, FLT, of Abhyankar.  These are equivalent to projective linear transformations of the projective plane.

Let A = {{a1, a2, a3}, {b1, b2, b3}, {d1, d2, d3}}.  The FLT U associated with A is  U ({x, y}) = { ((a1x+a2y+a3)/(d1x+d1y+d3)), ((b1x+b2y+b3)/(d1x+d2y+d3))}.   FLT are compatible with matrix multiplication, that is if U,V are FLT associated with A,B respectively then V(U{x,y}) is the FLT associated with the matrix BA.  In particular if A is invertible so is U and the inverse of function U is the FLT associated with Inverse[A].

We have procedures to implement the FLT on points, let

A={{0.707107`,0.`,0.`},{0.`,-0.707107,-0.707107},{0.`,0.707107`,-0.707107}};
flt[{.6,.8},A]

flti[{0,1},A]

Here flti acts on and returns infinite points if necessary, while flt works only with affine points.  For action on a curve we have 
	g = FLT[f,A,{x,y}]
 which is more complicated and not implemented on this worksheet.  A classical approach requires substitution of the inverse function and simplification.  However simplification of rational functions is not practical when floating point coefficients are used.  The author uses instead his duality method which is too lengthy for this worksheet.

Fractional Linear Transformations are useful for viewing infinite points, that is bringing them into the affine plane.  One can see they are then indistinguisible from finite points so there is an obvious topology to put on the union of a curve and its infinite points.

Application of FLT to Conics and Cubics

One can identify the type of a non-singular (not union of lines) conic by counting infinite points,  zero gives an ellipse, 1  a parabola and 2 gives a hyperbola.

The following claims can be shown constructively, there are algorithms that produce the desired result:

1) Every non-singular plane conic can be reduced to one of the forms  a^2 x^2 ±  b^2 y^2 = 1with a<b, or y=ax^2, with a>0  using Euclidean FLT only.  A Euclidean FLT is one with last row {0,0,1} and upper left 2x2 submatrix orthogonal.

2) Every non singular plane conic can be reduced to both  x^2+y^2-1 and y-x^2 by a FLT.

From the later  we get  parameterizations of the conic for x real.  in fact we have the following correspondance for a non-singular matrix A

f=FLT[y-x^2,A,{x,y}]    ⟺ f parameterized by x-> flt[{x,x^2},Inverse[A]], x Real

f=FLT[x^2+y^2-1,A,{x,y}] ⟺ f parameterized by t-> flt[{Cos[t],Sin[t]},Inverse[A]], 0⩽ t ⩽2ℼ

For cubics we have

 Every numerical irreducible cubic (no lines) can be reduced by FLT to the form  y^2=x^3+ax^2+bx+c.  See my web package.
 
 This technique played a key role in my recent paper on Quadratic Surface Intersection Points.   See my Northeastern website http://homepages.neiu.edu/~bhdayton

Decompositions

A real curve can be decomposed into 1) affine topological components (ignore infinite points), 2) projective components (include infinite points), 3) algebraic components or 4) ovals.  

Affine topological decomposition is the main computational step for 2) and 4).  The obvious tool would be curve tracing.  But curve tracing, the main engine of homotopy continuation, works well in the complex domain but poorly in the real domain.  So we use an other method.

For plane curves only, an algebraic component is a subset which itself is a complete solution set of a polynomial equation.  In fact a polynomial describing the component will be a factor of the original curves.  For curves with integer coefficients Mathematica’s factor is enough, for numerical coefficients we have a one point factorization procedure  CFactor which finds the algebraic component containing a given point.  This is mostly useful in degrees 2 and 3.

An oval is a subset of an algebraic curve which is a nice, almost 1-1, topological image of a circle.  Gauss’ oval theorem states that every real algebraic curve is a union of ovals.  We describe a data type of numerical data which describes an oval decomposition of a real curve.  We have an oval decomposition routine which can often decompose non-singular  curves of degree 4 or less into ovals.  Additionally there is a suite of procedures for working with and sketching this oval data.  With these tools one can, manually construct oval data for many singular curves also

Some Tools available
	OvalDecomp[f,a,{x,y}]  	Decomposes non-singular curve into ovals (no infinite tangents).
		OvalQ[f,p,OD,{x,y}]    	Decides which oval contains point p and where.
		OvalRefine[f,OD,{x,y}] 	Approximately doubles number of points in decomposition.
		flto[f,OD,A,{x,y}]		Transforms decompositon to a decomposition of image curve under FLT A
		
	See notebooks  OvalTools.nb and OvalExamples.nb

Graphics

realinflectionpoints[f_,X_]:=Module[{H,fh,d,tt,rp,L,X1,sol,p},
d=tdeg[f,X];
fh=Expand[tt^d*(f/.Thread[X->1/tt*X])];
X1=Append[X,tt];
H=Det[D[fh,{X1,2}]];

rp=Cases[X/.NSolve[{fh,H,tt-1}],{_Real,_Real}];
If[Length[rp]>0,L= rp,Print["no point found"];Return[]];
rp=Reap[Do[
If[Norm[{D[f,X[[1]]],D[f,X[[2]]]}/.Thread[X->L[[i]]]]> .01, Sow[L[[i]]]], {i,Length[L]}]][[2]];
If[Length[rp]>0,rp[[1]],{}]
];
autocontourplot[f_,M_,X_]:=Module[{ips,bps,tps,cp1,cp2,cp3,cp4,cp5,sol2,stps,clines, hps,lpts,cpts,k,xmin,xmax,ymin,ymax},
k=RandomReal[{-1,1},2].X;
tps=tangentrealpoints[f,k,X,False];
Quiet[hps=realinflectionpoints[f,X]];
stps=Sort[Table[k/.Thread[X->tps[[i]]], {i,Length[tps]}],Greater];
clines=Append[Prepend[Table[(stps[[i]]+stps[[i+1]])/2,{i,Length[stps]-1}],  stps[[1]]+.25*Abs[stps[[1]]]],stps[[-1]]-.25*Abs[stps[[-1]]]];
lpts={};
Quiet[Do[lpts=Join[Cases[X/.NSolve[{f,k-clines[[i]]}],{_Real,_Real}],lpts],{i,Length[clines]}]];
If[Length[hps]>0,cpts=Join[lpts,hps],cpts=lpts];
If[Length[tps]>0,cpts=Join[cpts,tps]];
xmin=Max[-M,Floor[Min[Flatten[Take[cpts,All,1]]]-1]];
xmax=Min[M,Ceiling[Max[Flatten[Take[cpts,All,1]]]+1]];
ymin=Max[-M,Floor[Min[Flatten[Take[cpts,All,-1]]]-1]];
ymax=Min[M,Ceiling[Max[Flatten[Take[cpts,All,-1]]]+1]];
Print[{X[[1]],xmin,xmax},{X[[2]],ymin,ymax}];
Print[ContourPlot[{f==0}, {x,xmin,xmax}, {y,ymin,ymax}]];
];

For a quick view of a curve use autocontourplot

f=NewtonHyperbola[425];
autocontourplot[f,15,{x,y}]

The  procedure ATCsketch finds the affine topological components.  Its engine is the basis for the other decompositions.  For example on the right is the ATCsketch of x^3-x-y^2 and on the left of NewtonHyperbola[425]

Oval Sketches

A sketch is a symbolic representation rather than literal view of a curve.  We have three types of oval sketches.  OvalFSketch shows the curve on a finite representation of the projective plane, here antipodal points on the red circle are equal infinite points.  The left hand sketch shows the graph of the rational function  y=1/x^2(x-1), the middle shows the 4 ovals of NewtonHyperbola[1876] and  the right shows the 3 topological components of NewtonHyperbola[613]

The OvalBoxSketch is an ordinary representation similar to a ContourPlot but shows different ovals in different colors.  The first shows the 4 ovals of NewtonHyperbola[1876], the second is a detail of the singularity of the rational curve above while the right hand is an ordinary contour plot of this singularity.

The OvalMSketch (still somewhat experimental) shows a representation of a single oval on a Moebius band model of the projective plane.  The outer boundary represents the origin while the middle circle represents the infinite line.  The left hand picture is of the parabola y = (x^2) + 1 while the right shows the oval NewtonHyperbola[613] which crosses the infinite line at three points.