A Numerical Approach to
Real Curves
with the
Wolfram Language

Barry H Dayton
Northeastern Illinois University (Emeritus)

barryhdayton.space

Join the Conversation #WolframTechConf

This book explores how one can do real plane algebraic curve theory  using the Wolfram language.  Traditionally this is a subject that required a fair amount of mathematical prerequisites.  But replacing the abstract concepts with Wolfram Language this should now be accessible to those with a general mathematical background.  This book will be published by Wolfram Media.

(RowDefault[, ({{(RowDefault)}, {(RowDefault)}})])

Part 1, Plane Curves

A plane curve is the full solution set of a single equation with a  two variable polynomial on the left and zero on the right.  For example

x^3-x^2y+2x y^2+3x y+ 2y^2-x+3 ==0

An example is a line.  The equation of a real line can be given by solving for the coefficients of a x+b y+c-y=0.

rrv=RandomReal[{-1,1},3];
line[p_,q_,x_,y_,rr_]:=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,rr};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).{x,y,1}
];

Here p,q could be points {(Subscript[x, 1,])(Subscript[y, 1])}, {(Subscript[x, 2,])(Subscript[y, 2])} for (Subscript[x, 1,])(Subscript[y, 1,])(Subscript[x, 2,])(Subscript[y, 2])  reals or perhaps one could be a slope {1,m,0} or perhaps one could be a tangent direction {u,v,0}.  We will think of the triples as infinite points.

line[{1,2},{3,5},x,y,rrv]

First note that the coefficients are not the expected integers but Mathematica machine numbers.  Since most, or sometimes all, points on an algebraic curve are not integers I will be working numerically throughout this talk.

We will think of the triples with last coordinate zero as infinite points.  Note the code converts point {(Subscript[x, 1,])(Subscript[y, 1])} to {(Subscript[x, 1]),(Subscript[y, 1]),1}.  All these triples with at least one coordinate not zero are considered to be points in the projective plane.  Note that replacing {(Subscript[x, 1]),(Subscript[y, 1]),1} by r {(Subscript[x, 1]),(Subscript[y, 1]),1}={r (Subscript[x, 1]),r (Subscript[y, 1]),r}, r≠0, does not change the line.  Thus projective points are homogeneous  {(Subscript[x, 1]),(Subscript[y, 1]),1}={r (Subscript[x, 1]),r (Subscript[y, 1]),r}.

Curves (continued)

Curves could be quadratics

f=x^2-1/4 y^2 - 1

cubics

g=x^3-2 x^2+x-1- (y^2)

or higher degree.  Two families of examples I like are

Gaussian curves from Gauss’ first and fourth proof of the Fundamental Theorem of Algebra (1799-1849)

GaussCurve[p_,x_,y_,z_]:=ComplexExpand[Re[p/.{z->x+I*y}]]

where p is a complex polynomial in z.

g=GaussCurve[z^5-5z^4+9z^3-5z^2,x,y,z]
ContourPlot[g==0, {x,-5,5},{y,-5,5}]

Curves (continued)

Newton Hyperbolas are based on some of Newton’s explorations with cubics.  The code is  hidden here.  Enter an integer from 1 to (2^28).

f=NewtonHyperbola[1876,x,y]
ContourPlot[f==0,{x,-20,40},{y,-10,10}]

NewtonHyperbola[n_,x_,y_]:=Module[{v,c,k,n2dig,s},
If[(n>2^28) ||(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},{6,0},{5,1},{4,2},{3,3},{2,4},{1,5},{0,6}};
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,.2^15,.2^10,.2^6,.2^3,.2^6,.2^10,.2^15};
k=Which[ n<8,3,n<64,6,n<1024 ,10,n<32768,15,n<2^21,21,True,28];
n2dig=IntegerDigits[n,2,k];
s=Reverse[n2dig];
FromCoefficientRules[Table[v[[i]]->(-1)^s[[i]]*c[[i]], {i,k}],{x,y}]]

These curves in two variables are called affine curves, to get all the projective points we can homogenize.

Homogenize[f_,x_,y_,z_]=Module[{h,AA,PA,d},
AA=<|CoefficientRules[f,{x,y}]|>;
d=Max[Total[#]&/@Keys[AA]];
PA=Table[Append[k,d-Total[k]]->AA[k],{k,Keys[AA]}];
FromCoefficientRules[PA,{x,y,z}]];

Homogenize[x^2-2x y+y^2+x-2,x,y,z]

Points - Numerical Viewpoint

The concept of a point on a curve, that is, a numerical solution of an equation f=0 where f is a bivariate polynomial s is an oxymoron since zero is not a machine number.  Evaluating the left hand side gives instead a small residue.  For example

f=x^2-3 x y +y^4-2;
f/.Thread[{x,y}->{-2.659583803082112`,-0.6595838030821121`}]

Defining a point by putting some condition on the residue is not satisfactory.  Rather we will define a point as a limit of a quadratically convering iterative system solving algorithm, such as Newton’s Method, where convergence is not the size of the residue but that further iterations do not change the point.  In particular we will consider the output of Mathematica’s  NSolve or FindRoot satisfactory for  identifying points and also for identifying which points are real.

We can find points by intersecting a curve with a line or another curve.  A more important family of points are the critical points.

CriticalPoints[f_,x_,y_]:={x,y}/.NSolve[{f,y D[f,x]-x D[f,y]},{x,y},Reals]

CriticalPoints[y^2-x^3+3x-5,x,y]

A critical point which is a multiple solution, i.e. it appears more than once in the list, is called a singular point.

CriticalPoints[y^2-x^3-x^2,x,y]

Here {0, 0} is a singular point of the curve (y^2)=(x^3)+ (x^2)

Points - Gauss Viewpoint

Gauss studied real curves an applied this to his first and fourth proofs of the Fundamental Theorem of Algebra.   To Gauss the curve f=0 is the complete topological boundary of the set f>0.  Moreover this curve can be considered oriented so that when the curve travels in the positive direction the set f>0 lies on the right.

Technical Note:  It is assumed the curve f is square-free, but this is easy to check/fix.

Two applications of this

The curve can be plotted using ContourPlot or, better, RegionPlot.

One can path trace the oriented curve.  There are several methods that can be used including a simple Euler' s method and a more complicated version using Mathematica' s NDSolveValue.  Typically one can trace from one critical point to the next on the curve.  However one cannot trace out of a singular point, but tracing into a singular point works well, replace f by -f to change direction.  Extra points may be needed.

The code and  data for curve f=(x^2)y-(x^4)-(y^3) is hidden.  Region plot and tracing are shown.  The red dots are the critical and other endpoints in tracing.

f=-(x^4-x^2y +y^3);
ap=
{{-0.38089187219472304`,0.23798525400275947`},{0.38089187219472304`,0.23798525400275947`},{0,0},{0.33035355027080476`,-0.3753219042681578`},{-0.33035355027081426`,-0.37532190426816436`},{-0.10616104058422673`,0.1`},{0.10616104058422673`,0.1`}};
L=
{{-0.33035355027081426`,-0.37532190426816436`},{-0.28003098262573173`,-0.31313734855836123`},{-0.2289050287233655`,-0.25161301083080767`},{-0.17688963584897385`,-0.19084075078656212`},{-0.12387910132995686`,-0.13093741339821693`},{-0.06974058257555478`,-0.07205588341692795`},{0,0},{0,0},{-0.06286470080985908`,0.003967776827250588`},{-0.10222034463525448`,0.010561753223017011`},{-0.14092269658220552`,0.020279146286835897`},{-0.17871815069601932`,0.03307277192952807`},{-0.21534561911033015`,0.04889433316423835`},{-0.25052247770886366`,0.06770697302898489`},{-0.28391608781453126`,0.08950310684798661`},{-0.31508378293312966`,0.11433161827651463`},{-0.3433294543410152`,0.1423418103843027`},{-0.36727633820168987`,0.17383372643224745`},{-0.38310850279833253`,0.2087489402329529`},{-0.38089187219472304`,0.23798525400275947`},{-0.38089187219472304`,0.23798525400275947`},{-0.35797101154339395`,0.24983319754485345`},{-0.32098696028421325`,0.24396631169995922`},{-0.2844444765315698`,0.22864399176634476`},{-0.24963640521979402`,0.20916893821714277`},{-0.21620843184247313`,0.18729349908217702`},{-0.18386943468922065`,0.16379872508595847`},{-0.15241852316476037`,0.13911065202159864`},{-0.10616104058422673`,0.1`},{-0.10616104058422673`,0.1`},{-0.07639471324543984`,0.07328992709528848`},{-0.04715536897866882`,0.04600155184979335`},{0,0},{0,0},{0.06286470080985908`,0.003967776827250588`},{0.10222034463525448`,0.010561753223017011`},{0.14092269658220552`,0.020279146286835897`},{0.17871815069601932`,0.03307277192952807`},{0.21534561911033015`,0.04889433316423835`},{0.25052247770886366`,0.06770697302898489`},{0.28391608781453126`,0.08950310684798661`},{0.31508378293312966`,0.11433161827651463`},{0.3433294543410152`,0.1423418103843027`},{0.36727633820168987`,0.17383372643224745`},{0.38310850279833253`,0.2087489402329529`},{0.38089187219472304`,0.23798525400275947`},{0.38089187219472304`,0.23798525400275947`},{0.35797101154339395`,0.24983319754485345`},{0.3209869602842136`,0.24396631169995875`},{0.28444447653156996`,0.22864399176634478`},{0.2496364052197942`,0.20916893821714283`},{0.21620843184247326`,0.1872934990821771`},{0.18386943468922076`,0.16379872508595855`},{0.15241852316472104`,0.13911065202157086`},{0.10616104058422673`,0.1`},{0.10616104058422673`,0.1`},{0.07639471324543984`,0.07328992709528848`},{0.04715536897866882`,0.04600155184979335`},{0,0},{0,0},{0.06974058257554722`,-0.0720558834169199`},{0.12387910132994946`,-0.13093741339820872`},{0.17688963584896655`,-0.19084075078655377`},{0.22890502872335833`,-0.2516130108307992`},{0.2800309826257247`,-0.3131373485583526`},{0.33035355027080476`,-0.3753219042681578`}};
{RegionPlot[f>0,{x,-.5,.5},{y,-.5,.5}],Show[ContourPlot[x^2+y^2==.25,{x,-.5,.5},{y,-.5,.5},ContourStyle->LightGreen],ListLinePlot[L,PlotStyle->Black],ListPlot[ap,PlotMarkers->Automatic,PlotStyle->Red]]}

Fundamental Theorem

Every curve restricted to a bounded region can be traced as the previous example from critical points, perhaps augmented by a few more. We can produce a Graph, directed or undirected,  where the vertices are these points  and edges represent  the arcs traced.  For the tracing above the graph would be

Graph[{1->2,2->3,3->4,4->2,2->5,5->6,6->2,2->7},VertexLabels->"Name"]

The observation is that each vertex is even except the boundary ones.  If we let the curve be unbounded and add the infinite points then we get the Fundamental Theorem:

The Graph of a plane algebraic curve in the projective plane is an Euler Graph, that is, every vertex is even.

Gauss identified the infinite points of Gauss Curves.  Here are some projective graphs of Gauss Curves.

{({{ (z^3)-1}, {}, {}}), ({{z((z^2)-4)((z^2)+z+1)}, {}, {}}),({{5 (z^2)-5 (z^3)+(z^5)}, {}, {}}),({{-5 (z^4)+(z^5)}, {}, {}})}

Intersection Multiplicity, Bézout and Caley

Intersection Multiplicity: I use numerical linear algebra to find the intersection multiplicity an intersection point of 2 curves using a Wolfram Language implementation of Macaulay’s 1916 Algorithm

Bézout’s Theorem:  Curves of degree m,n  respectively intersect in the complex projective plane in exactly  m n  points counting multiplicity.  (note:  Not proven here but all intersection points can be found.)

An application is that an irreducible curve of degree d has at  most  (((d-1)(d-2))/2) singular points. We show this by implementing an algorithm in the Wolfram Language factoring a curve with  (((d-1)(d-2))/2)+1 or more singular points showing that it is not irreducible.

Caley’s Theorem (1865):  An irreducible curve with  (((d-1)(d-2))/2) double points can be parameterized by rational parameters.

Example:   The hypocycloid

hc=Expand[3 ((x^2)+(y^2))^2+8 x (-(x^2)+3 (y^2))+6 ((x^2)+ (y^2))-1]

can be parameterized by

X=-((0.3946283194077289` (-0.021739216895879716`+0.09692422977782154` t+0.7836858499940199` (t^2)+1.2150140700632728` (t^3)+1.` (t^4)))/(0.01843408157341412` +0.11686381789887891` t+0.4567603971204911` (t^2)+0.8607347125368988` (t^3)+1.` (t^4)));
Y=-((0.5577361971417719` (0.0046462193954929015` +0.09876727042260801` t+0.7203891230964178` (t^2)+1.9396274324562477` (t^3)+1.` (t^4)))/(0.01843408157341412` +0.11686381789887891` t+0.4567603971204911` (t^2)+0.8607347125368988` (t^3)+1.` (t^4)));

Show[ContourPlot[hc==0,{x,-2,2},{y,-2,2},ContourStyle->Directive[Thick,Orange]],ParametricPlot[{X,Y}, {t,-10,10},PlotStyle->Directive[Dashed,Black]]]

Representations -- Gauss Hemisphere, Möbius Band

For the Gauss Hemisphere we map onto the interior of a disk preserving arguments while the boundary of the disk consists of the infinite points not identified with their antipodal points.  Often an augmented ContourPlot works better, features far away are denoted with dashed curves.  These are shown below for NewtonHyperbola[1 050 752].

{                                                         }

The projective plane may be topologically modeled by a Möbius band where the entire boundary is pinched together to represent the affine origin {0,0}.  The center line of the band represents the infinite projective line.  This model can be implemented by twisting the strip and gluing the green sides and two halves of the orange top together.

Some plots in this model are a line, not through origin, GaussCurve[(z^3)+3] and rational function y=(4/(x((x^2)-4))).

{,,}

Descartes - Viro and Gauss-Viro Diagrams.

A Descartes-Viro diagram of a curve f=0 is obtained by coloring lattice points {i,j} in the first quadrant, including positive axes, green for positive coefficients of  (Subscript[c, i,j])(x^i)(y^j) and red for negative.  In other quadrants the color depends on sign of the term evaluated at a typical point in the interior of the quadrant.  A curve is drawn separating the colors.  Note, middle graphic, sometimes this is ambiguous.  These give possible curves that can be obtained from these coefficient signs.

{}

A Gauss-Viro Diagram is a careful modification of the Descartes-Viro diagram showing actual positive and negative regions of a specific curve.  I claim that every plane algebraic curve has such a diagram (generalized Harnack Theorem). The Descartes-Viro and Gauss-Viro diagrams of  curve NewtonHyperbola[1 050 752] of the previous slide are, letters representing infinite points,

Plane Curves Conclusion

The main goal of this approach is that with this information we should be able to identify a plane curve given possibly by an equation with machine number coefficients, identify this topologically (Graph) and geometrically (Gauss-Viro Diagram).  If necessary trace the curve to make it portable out of the plane curve setting.

Steps to Analyze Curve

Construct Descartes-Viro Diagram

Find Critical points

Find infinite points  (set homogenized equation to z=0, solve for x,y with auxilliary equation).

Draw ContourPlots where useful.

Path trace if necessary, especially into singular points, more points may be needed.

Draw Graph and Gauss-Viro Diagrams and augmented ContourPlot and/or ListLinePlot

Part 2, Space Curves

Naive Space Curves in (ℝ^3)

In multivariable calculus a space curve is given by two equations in three variables.  This situation is simple enough to borrow ideas from our study of plane curves thanks to the availability of Cross product.  We can easily compute a tangent vector

TangentVector3D[{f_,g_},p_,x_,y_,z_]:=Module[{n1,n2,bi},
n1={D[f,x],D[f,y],D[f,z]}/.Thread[{x,y,z}->p];
n2={D[g,x],D[g,y],D[g,z]}/.Thread[{x,y,z}->p];
bi=Cross[n1,n2];
If[Norm[bi]<.0001,Echo[p,"No tangent vector at "]; bi,Normalize[bi]]]

or a critical point

CriticalPoints3D[{f_,g_},x_,y_,z_]:=Module[{J},
J=D[{f,g,x^2+y^2+z^2},{{x,y,z}}];
{x,y,z}/.NSolve[{f,g,Det[J]},{x,y,z},Reals]]

We can then path trace to fill out the curve.  Note that the Fundamental Theorem holds for real space curves.

General Space Curves

More generally I will consider curves in (ℝ^n) for arbitrary n≥3 given by a system of  k ≥n-1 polynomials in n variables.  These will require heavy numerical linear algebra using Macaulay and Sylvester matrices and duality.

An important tool is a function which uses a reduced Macaulay matrix to determine if a point has a unique tangent vector and, if so, find it.

A polynomial system of equations defines a space curve if at least solution one point has a unique tangent vector and at most finitely many points  do not have tangent vectors.

The technique for analyzing space curves is to project the equation to a plane curve.  Analyze and trace the plane curve and then use fiber lifting to describe a space curve in (ℝ^3).

Example 1:  Cyclic-4 curve,  C4={w+x+y+z,w x+x y+w z+y z,w x y+w x z+w y z+x y z,-1+w x y z}⊂(ℝ^4)

C4 system projects to a degree 6 reducible plane curve, two hyperbolas with common asymptotes and two lines which are determined to be ghost lines as they are not part of the point projection.  It is determined lift to (ℝ^3)or (ℝ^4) lies in a plane so we are done.  However it is noticed that {w+y,x+z,-1+(w^2) (x^2)} is a simpler system describing this curve.

Example 2:  A sytem of 3 random quadratic hypersurfaces through origin of (ℝ^4).  This projects to a degree 8 plane curve with 7 singular points all of which occur because 2 points of the original project to the same plane point.  The original curve is non-singular.  This plane curve is analyzed and traced with some difficulty, partly because of the high degree.  This curve is then lifted to (ℝ^3) as the solution of a system of 7 equations of degrees 4 and 5 in three variables.  On the left is the projection to the plane, the right gives part of the curve in 3-space, note the green and blue components are disjoint from each other.