Gauss Decomposition

In this section, we discuss our method for decomposing curves inside a circular region. The idea behind this comes from C. F. Gauss’s fourth proof of the fundamental theorem of algebra (FTA) written in 1849. Let’s start with an exploration of some of the topographical terminology.

4.1 Plane Topology

To avoid too much abstraction, we will work in the real Euclidean plane (ℝ^2) viewed as a metric space. The norm of a point is ||p||=||{(Subscript[x, 0]),(Subscript[y, 0])}||=(Sqrt[(Subscript[x, 0]^2)+(Subscript[y, 0]^2)]), and the distance between two points p, q will be denoted by ||p-q||.

Norm[{3.2,-5.6}-{-1.3,3.8}]

EuclideanDistance[{3.2,-5.6},{-1.3,3.8}]

In this subsection, 𝒮 will denote a subset of (ℝ^2), often called a curve or a semialgebraic set. A semialgebraic set is a subset of the plane that is obtained by finite unions, complements and/or intersections of sets of the form {{x,y}|f(x,y)>0} or {{x,y}|f(x,y)≥0} where f[x,y] is a polynomial function in variables x, y—that is f(x,y)—is a curve. Mostly, we will just use the atomic semialgebraic set {{x,y}|f(x,y)>0}. Real Algebraic Geometry (Bochnak, 1998), one of the most comprehensive texts in the eponymous field, points out that semialgebraic sets are very useful in studying algebraic curves, a fact that Gauss noted in his proof written in 1849.

A subset 𝒰 of 𝒮 is called open if, for every point p in 𝒰, there is a number ϵ>0 so that every point q in 𝒮 with ||p-q||<ϵ is also in 𝒰.

A point p of (ℝ^2) is called a boundary point of 𝒱⊂ 𝒮 if, given ϵ>0, there are points u∉𝒱, v∈𝒱 with ||u-p||<ϵ and ||v-p||<∈. Here, p itself may or may not be in 𝒱. The set 𝒱 is called closed in 𝒮 if it contains all its boundary points in 𝒮. A standard elementary result in topology is that 𝒱 is closed in 𝒮 if and only if the complement of 𝒱 in 𝒮 is open in 𝒮. The set consisting of 𝒱 and all its boundary points in 𝒮 is called its closure in 𝒮. A point in 𝒮 that is not a boundary point is called an interior point.

The set 𝒰={{x,y}|(x^2) + (y^2) < 1} is open in (ℝ^2) while 𝒱={{x,y}|(x^2) + (y^2) ⩽ 1} is closed in (ℝ^2) and is the closure of 𝒰.

The set 𝒞={{x,y}|(x^2) + (y^2) = 1} is the boundary of both 𝒰 and 𝒱 and is closed in (ℝ^2).

While 𝒲={{x,y}|(x^2) + (y^2) =1, y>0} is not open in (ℝ^2), it is open in 𝒞 with closure {{x,y}|(x^2) + (y^2) = 1, y ⩾ 0 }.

Additional easy exercises in point set topology are to show an arbitrary union of open sets in 𝒮 is open in 𝒮 while finite unions of closed sets are closed. Using complements, finite intersections of open sets are open but infinite intersections of closed sets are closed. Both the empty set and 𝒮 itself are open and closed in 𝒮.

The key definition of this subsection is:

The set 𝒮 is connected if and only if 𝒮 is not the union of two disjoint open subsets of 𝒮.

Alternatively, 𝒮 is connected if and only if the only subsets of 𝒮, which are both open and closed in 𝒮, are the empty set and 𝒮 itself. If p∈𝒮, then the union of all connected subsets of 𝒮 containing p are called the connected component of p in 𝒮.

Let 𝒞={{x,y}| (x^2) + (y^2) = 1} be the circle curve. Then the complement of 𝒞 in (ℝ^2) is the union of two open sets {{x,y}|(x^2) + (y^2) < 1}⋃{{x,y}|(x^2)+ (y^2) > 1}. Each of these is a connected component of the complement of 𝒞.

Another famous example of a disconnected open set of the plane comes from an analysis of the complex Newton method basin of attraction. Consider the complex polynomial p[z]=(z^3)-1, and let 𝒰 be the set of points p in (ℝ^2)such that Newton’s method, starting from initial point p, converges. It can be shown that 𝒰 is open in (ℝ^2). Let 𝒱 be the complement of 𝒰, which is then closed.

JuliaSetPlot[z-(z^3-3)/(3z^2),z, ColorFunction->"RedBlueTones"]

𝒰 consists of the points in the "RedBlueTones" colors, and 𝒱 consists of the white points where Newton’s method does not converge. If you could blow this picture up sufficiently, you would see that the complement 𝒰 of 𝒱 is the disjoint union of infinitely many connected components! It is also interesting to note that the white boundary of 𝒱 is neither one- nor two-dimensional; it is what is known as a fractal.

Now let 𝒮, 𝒯 be plane subsets, perhaps one or both is the real line ℝ. A function ϕ: 𝒮⟶𝒯 is called continuous if for every open set 𝒱 ⊂ 𝒯, the inverse image set (ϕ^(-1))(𝒱)={p ∈𝒮|ϕ(p)∈𝒱} is open in 𝒮. In this case, both 𝒮, 𝒯 are real intervals, which is equivalent to the Calculus I definition of continuity. By simple set considerations, the definition is equivalent to the statement that the inverse image of every closed set is closed.

This last statement shows that if f(x,y) is viewed as the polynomial function f:(ℝ^2)⟶ℝ, then the point set V(f)={p|f[p]=0} is closed in (ℝ^2) since {0} is closed in ℝ. Further, its complement {p|f(p)>0 or f(p)<0} is the union of two open sets in (ℝ^2). Thus a main application of the topological considerations above can be summarized as follows:

A square-free curve, f=0, as a point set V(f) disconnects its complement into two (not necessarily connected) open sets {p|f(p)>0}⋃{p|f(p)<0}.

The reason f needs to be square-free is that otherwise, one of the two sets above could be empty. It is somewhat harder to prove that in the square-free case, the curve is the complete boundary of each of these semialgebraic sets in (ℝ^2). We leave it to Gauss to make the pronouncement that this is true in the next subsection. It actually follows easily from the result in Chapter 5 that a square-free curve has a finite number of singular points.

We have already seen examples where the semialgebraic set {p|f(p)>0} is not connected. Nor is it true that the set V(f) is connected; however, since V(f) is a closed set of (ℝ^2), each component of V(f) is too.

For later use, the following theorem of point set topology is easy to prove in our context of open sets:

The continuous image of a connected set is connected.

In particular, the image of a parametrically defined curve ϕ(t) from a connected interval of ℝ to (ℝ^2) is connected. It is the main goal of this section to show that not only is each connected component of the boundary of each component of {p|f(p)>0} the image of such a continuous parametric curve, but to actually find a numerically accurate approximation of these boundary curves. Thus curves and semialgebraic sets are much better behaved than arbitrary sets (such as the Newton basin of attraction).

4.2 The Gauss Principle

D’Alembert gave what the French considered to be the first acceptable proof of the fundamental theorem of algebra in 1748, but Gauss received credit for said proof (in the general world) in 1799. At the time of his fourth proof, Gauss was in his 70s, one of many examples that creativity in math does not end at 40. Overall, this last proof was similar to the first but with one new idea clarifying his earlier argument. I will call this the Gauss Principle. I state this using modern mathematical terminology:

Let 𝒞 be the point set given by the solution set of the square-free equation f(x,y)=0, and let 𝒟={{x,y}|(x^2)+(y^2)≤(r^2)} be the disk of radius r. Then the complete boundary of the semialgebraic set 𝒮={p∈𝒟|f(p)>0} consists of 𝒞 ⋂ 𝒟 together with portions of the boundary (x^2)+(y^2)=(r^2) of 𝒟. 𝒮 has finitely many connected components. Further, if 𝒦 is any connected component of 𝒮, then the connected components of the boundary of 𝒦 are images of continuously parameterized curves, which can be parameterized in the direction so that as the parameter increases, the interior of 𝒦 lies to the right, that is the regions are parameterized clockwise.

For low-degree or well-behaved nonsingular curves, it is easy to view this using RegionPlot. Note that {p ∈ 𝒟|f(p)>0} is the darker orange-blue region. For convenience, in the rest of this chapter, we define the utility function gaussPlot.

gaussPlot[f_,r_,pts_]:=Module[{plta,pltb,pltc},
plta:=RegionPlot[x^2+y^2<r^2,{x,-1.1r,1.1r},{y,-1.1r,1.1r}, PlotStyle->LightOrange, BoundaryStyle->Orange];
pltb:=RegionPlot[{f>0&& x^2+y^2<r^2 },{x,-r,r},{y,-r,r},BoundaryStyle->Blue];
pltc:=ListPlot[Table[Labeled[pts[[n]],n],{n,Length[pts]}],PlotStyle->Red, PlotMarkers->Automatic];
Show[plta,pltb,pltc, ImageSize->Small]];

{gaussPlot[y^2-x^3+x,3,{}],gaussPlot[-1-3 (x^2)+(x^3)+3 (y^2)-3 x (y^2),6,{}],gaussPlot[-1+x-0.2*x^2+0.008*x^3+y+x*y-0.2*x^2*y-0.2*y^2-0.2*x*y^2+0.008*y^3,15,{}]}

Note that there are one, three and four components of {p ∈ 𝒟|f[p]>0}, respectively.

f=(y-.75x-4)(y+.75x+4)(y-2x)(y+2x);
P0:=RegionPlot[x^2+y^2<100,{x,-10,10},{y,-10,10},PlotStyle->LightOrange, BoundaryStyle->None];
P1:=RegionPlot[y>.75x+4 && y>2x && y> -2x && x^2+y^2<100,{x,-10,10},{y,-10,10},PlotStyle->LightGreen,BoundaryStyle->Green];
P2:=RegionPlot[y<.75x+4 && y<2x && y>-2x && y>-.75x-4&&x^2+y^2<100, {x,-10,10},{y,-10,10},PlotStyle->LightBlue,BoundaryStyle->Blue];
P3:=RegionPlot[y<.75x+4&&y<-2x &&y>-.75x-4 &&y>2x && x^2+y^2<100,{x,-10,10},{y,-10,10}, PlotStyle->LightGray,BoundaryStyle->Gray];
P4:=RegionPlot[y<-.75x-4 && y>.75x+4 &&x^2+y^2<100,{x,-10,10},{y,-10,10},PlotStyle->LightPurple,BoundaryStyle->Purple];
P5:=RegionPlot[y<2x && y<-.75x-4 &&y<-2x && x^2+y^2<100,{x,-10,10},{y,-10,10},PlotStyle->LightMagenta,BoundaryStyle->Magenta];
P6=ContourPlot[x^2+y^2==100,{x,-10,10},{y,-10,10},ContourStyle->Orange];
{ContourPlot[f==0,{x,-10,10},{y,-10,10}],Show[P0,P1,P2,P3,P4,P5,P6]}

Note that the light orange region is f<0 while the other colored regions are the five connected components of the region f>0.

An important consequence of Gauss’s Principle is that the complement of a square-free curve f can be two-colored as it is the union of the two disjoint nonempty sets {f>0} and {f<0}.

Note that Abhyankar’s illustration (also in Chapter 1) violates this in several ways. More generally, there cannot be an odd number of curve segments meeting at a point.

Following the practice from Euclid until the end of the nineteenth century, Gauss most likely did not attempt to prove this topological statement as it was meant to be self evident. I am unaware of the original source for this statement; all modern accounts go back to the discussion in Appendix 1 of Uspensky (1948). I will not attempt a proof in this book but will give many examples and several algorithms illustrating this principle.

Here are some examples using Gauss curves, which Gauss had already introduced in his first proof. We use them to illustrate Gauss’s fourth proof.

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

f=gaussCurve[z^3-3z +1,x,y,z]
g=gaussCurve[z^3-1,x,y,z]
h= gaussCurve[z^3,x,y,z]
{gaussPlot[f,4,{}],gaussPlot[g,4,{}],gaussPlot[h,4,{}]}

Gauss noticed that for complex polynomial p(z) with leading term (z^n), the intersections of different Gauss curves with the circle of large radius r were nearly independent of the lower-degree terms and the plots near the circle were essentially the same.

WP:=RegionPlot[x^2+y^2<9,{x,-4,4},{y,-4,4},PlotStyle->White,BoundaryStyle->LightGray];
sol={x,y}/.NSolve[{x^2+y^2-15,h}];
bp=SortBy[sol,ArcTan@@#&];
{Show[gaussPlot[f,4,bp],WP],Show[gaussPlot[g,4,bp],WP],Show[gaussPlot[h,4,bp],WP]}

Focusing on the picture for p(z)=(z^n) (right plot), we see p((Subscript[x, j])+(Subscript[y, j])i) for j=1,3,5 has real part 0 and negative imaginary part while p((Subscript[x, j])+ (Subscript[y, j])i) for j=2,4,6 also has real part 0 (the definition of a Gauss curve) but imaginary part positive. If you extrapolate from the full pictures above to the numbered pictures in the previous output, we see that following the curve inwards from 1, 3 or 5 exits at one of the points 2, 4 or 6 (although which one depends on the lower-degree curve we see. In any case, Gauss observes that the imaginary part of p(x+yi) goes continuously from a negative value to a positive value along the trip. However, the image of the connected path to ℝ is connected, which means that it must pass through 0. At that point, p(x+yi)=0, so Gauss has found a zero of p(z) on the Gauss curves. In fact, since there are n-curves he has found n-zeros. This is his fourth proof of the FTA. I remind the reader that the argument in this paragraph only holds for Gauss curves of monic complex polynomials, but the Gauss Principle holds for all square-free curves.

There is a slight problem with the p (z)=(z^n) case because of the singularity at the origin, but the curve is removed, so there are still three components of {p|f(p)>0}. Therefore, this is essentially like the middle picture. More generally, using the Gauss Principle for singular curves, we can use perturbation to approximate the curve for RegionPlot.

In classical theory, a plane curve singularity must satisfy the three equations in two unknowns: f[x,y]=0,D [f,x]=0,D[f,y]=0, but an over-determined system tends to be very unstable, so any small change (perturbation) in the system causes a singular point to no longer be a point and/or singular. In our current semialgebraic point of view, adding a small positive constant to f expands the region f(p)>0 while subtracting a small positive constant contracts it. In either case, a pre-existing singular point is no longer on the perturbed curve.

p=z^5-5z^4+9z^3-5z^2;
f=-gaussCurve[p,x,y,z]
f1=f-.1;
f2 = f+.1;
sol={x,y}/.NSolve[{x^2+y^2-25,f}];
bp=Reverse[SortBy[sol,ArcTan@@#&]];
{gaussPlot[f,5,bp],gaussPlot[f1,5,bp], gaussPlot[f2,5,bp]}

Gauss’s example is on the left. The paths that Gauss is taking (that is, his decomposition of the left-hand singular system is approximated by the middle picture), so starting from 1, Gauss ends up at 10; from 3, he ends at 8; from 5, he ends at 4; from 7, he ends at 6; and from 9, he ends at 2. The picture on the right would give an alternate way to construct the paths; for example, from 1, he would end at 2. In fact, this is the way Gauss would decompose the curve -f=0. In either case, the original curve on the left would be the union of these paths.

The example shows that for singular curves, the Gauss Principle gives different decompositions for the curve f and -f. While this may cause a bit of confusion, the important thing is that the Gauss decomposition does allow for the decomposition of singular curves.

4.3 Strategy for Gauss Decomposition of Curves in a Disk

Based on the ideas of the last section, what we need to do to decompose an algebraic curve into continuously parameterized arcs (parameterized curve, that is, 1-1) is the following:

Find critical points on the curve that include the points where the curve intersects the boundary of the disk, singular points and at least two points on every topological connected component.

Connect these critical points with continuously parameterized arcs parameterized in the correct direction.

The critical points from Section 3.5 take care of most of the requirements of Step 1, but there may be critical points that are not in our disk because they are too far from the origin. Ideally, the radius r of the disk will be large enough to include all critical points, but in later sections, we will see examples of curves with critical points that have very large norms.

gPoints[f_,r_,x_,y_]:= Module[{bp,cp,sol},
sol=NSolve[{f,x^2+y^2-r^2},{x,y}, Reals];
If[Length[sol]>0,bp={x,y}/.sol,bp={}];
bp=Reverse[SortBy[bp,ArcTan@@#&]];
cp= DeleteDuplicates[ Select[{x,y}/.NSolve[{f,y D[f,x]-x D[f,y]},{x,y},Reals],Norm[#]<r&],Norm[#1-#2]<.1&];
Join[bp,cp]];

Note that for later convenience, we sort the boundary points clockwise around the circle and put them first on our list.

p=z^5-5z^4+9z^3-5z^2;
f=-gaussCurve[p,x,y,z]
cpf=gPoints[f,5,x,y];
f1=f-.2;
cpf1=gPoints[f1,5,x,y];
{gaussPlot[f,5,cpf], gaussPlot[f1,5,cpf1]}

g1=-(x^4+ 4x y +y^4 - 1/4 ) ; (* peanut curve *)
cpg1=gPoints[g1,3,x,y];
g2=-(1-6x^2+x^4-5y^2+5x^2y^2+4y^4);
cpg2=gPoints[g2,3, x,y];
{gaussPlot[g1,3,cpg1],gaussPlot[g2,3,cpg2]}

h1=-(1-8x^2+8x^4+3y-4y^3); (* Chebyshev quartic *)
cph1=gPoints[h1,2,x,y];
h2= -(x^4-x^2y+y^3) ; (*bow curve*)
cph2=gPoints[h2,.5,x,y];
{gaussPlot[h1,2,cph1],gaussPlot[h2,.5,cph2]}

We see on the Chebyshev quartic that points 5 and 6 and points 10 and 12 are very close together.

cph1a=Delete[cph1,{{6},{10}}];
cph2a=Join[cph2,{x,y}/.NSolve[{h2,x^2+y^2-.02},{x,y},Reals]];
{gaussPlot[h1,2,cph1a],gaussPlot[h2,.5,cph2a]}

Finally, we recall that choosing the positive semialgebraic set determines a direction for the curve—we go clockwise around each component of the positive region inside the disk. An important observation is that of the two unit tangent vectors at any regular point on a plane algebraic curve, we can select the one going in the correct clockwise direction with the following simple Mathematica procedure.

gTvec[f_,p_,x_,y_]:=Cross[Normalize[Chop[{D[f,x],D[f,y]}/.Thread[{x,y}->p],1.*^-6]]];

Note that this procedure returns {0,0} for a near-singular point; otherwise, it returns a vector of norm 1.

f=gaussCurve[z^4-2z^3-z+1,x,y,z]
cp=gPoints[f,5,x,y];
pltc=Graphics[Table[{Black,Arrow[{p,p+gTvec[f,p,x,y]}]},{p,cp}]];
Show[gaussPlot[f,5,cp],pltc, ImageSize->Medium]

Now that we have critical points and a direction of travel to complete Step 2, we need to be able to travel from one regular critical point to the next. We can treat singular points by removing them with a small perturbation or by finding nearby nonsingular points in the direction of travel at start (such as we did with the bow curve) and we can get these by using a straight line from the singular point to these nearby points.

4.4 Naive Curve Tracing

Curve (or path) tracking is a basic tool in numerical algebraic geometry (see Sommes and Wampler [2005] for examples). The simple idea, borrowed from numerical differential equation solving, is to use the predictor/corrector method we discuss here. Mathematica has good numerical differential equation solvers and we could use them for curve tracking (we will discuss this in a later section).

However, for the predictor/corrector method, given a point on the curve, we can make use of our known direction of travel. So at point p of the curve, we choose a step size s (which, for now, will be fixed) and predict that the next point will be near q=p+s gTvec[f,p,x,y]. We correct this by finding a nearby point to q on the curve. The most common method is to find the point of intersection of the line through q orthogonal to the tangent vector with the curve, typically using the two-variable Newton’s method. In practice, one or two iterations are a sufficient correction. In the Wolfram Language, we can use FindRoot.

Consider the curve f=(x^3)-(y^2)+3=0 and point (Subscript[p, 0])={1,2}.

f=x^3+3-y^2;
p0={1.,2.};
tv=gTvec[f,p0,x,y]

q=p0+1.*tv

nv=Cross[tv]

l=line[q, Append[nv,0],x,y]

lx={x,y}/.NSolve[{l,x==2}][[1]]
ly={x,y}/.NSolve[{l, y==2.9}][[1]]

p1={x,y}/.FindRoot[{f,l}, {x,q[[1]]},{y,q[[2]]}]

p2=closestPoint[f,q,x,y]

P1=RegionPlot[f>0,{x,.8,2},{y,1.7,2.9}];
P2=Graphics[{{Black,PointSize[Medium],Point[{p0,p1,p2}]},{Black,Text[Style["p0",16],p0+{.02,.08}]},{Black,Text[Style["p2",16],p2+{-.07,.02}]},{Black,Text[Style["p1",16],p1+{.04,.08}]}, {Black,Thick,Arrow[{p0,p0+{.8,.6}}]},{Black,Thick,Line[{lx,ly}]},{Black,Text[Style["tv",16],{1.5,2.3}]}}];
Show[P1,P2]

This method is implemented as pathFinderT in Appendix 2.

p2=closestPoint[f,q,x,y]

We can now write a procedure that, given a point p on the curve, will attempt to find a list of points on the segment to the next critical point q in the direction of the curve, assuming it is known by the user. The step size s is one of the arguments, and it is assumed that all the points are more than 2s from each other.

This procedure returns a list of points on the curve starting from the given initial point p and, hopefully, ending with q. Since we may miss it, we limit the number of iterations (this limit maxit can be adjusted in the code).

By default, this does 30 iterations and quits if it has not gotten to within distance 2s of the endpoint. An option maxit→m will change the maximum number of iterations. Note that even if this function does not achieve it’s objective, it will return a list of the the points it has found putting the endpoint q a the end of the list.

Options[pathFinder]={maxit->30};
pathFinder[f_,p_, q_,s_,x_,y_,OptionsPattern[]]:= Module[{L,p0,p1,u, mxt,k},
mxt=OptionValue[maxit];
p0=p;
L=Reap[Sow[p];
k=0;
While[Norm[q-p0]>2s && k<mxt,
u=p0+s*gTvec[f,p0,x,y];
p1 =closestPoint[f,u,x,y];
Sow[p1];
p0=p1;
k++];
Sow[q]][[2,1]];
If[k≥ mxt, Print[Style["Warning: iteration limit reached at",Orange],p1]];
L];

f=x^3+3-y^2;
L=pathFinder[f,p0,p2,.2,x,y]
Show[gaussPlot[f,4,{p0,p2}], ListLinePlot[L, PlotStyle->Black]]

In general, path finding is difficult, especially in real cases. Constructing a black box program, which will decompose any real plane curve into paths, is very difficult. Receiving guidance from a person who knows which points lie on a path can be helpful. If we want the computer to reliably and quickly figure out the curve by using path tracing to find the next point, the naive approach used previously will not work well. We will attempt an improvement in the next section.

f=gaussCurve[z^5-5z^4+9z^3-5z^2,x,y,z];

f1=f-.1

gp=gPoints[f1,5,x,y];
gaussPlot[f1,5,gp]

Timing[
L1=pathFinder[f1,gp[[2]],gp[[15]],.2,x,y];
L2=pathFinder[f1,gp[[15]],gp[[1]],.2,x,y];
L3=pathFinder[f1,gp[[4]],gp[[11]],.2,x,y];
L4=pathFinder[f1,gp[[11]],gp[[5]],.2,x,y];
L5=pathFinder[f1,gp[[6]],gp[[12]],.2,x,y];
L6=pathFinder[f1,gp[[12]],gp[[7]],.2,x,y];
L7=pathFinder[f1,gp[[8]],gp[[13]],.2,x,y];
L8=pathFinder[f1,gp[[13]],gp[[3]],.2,x,y];
L9=pathFinder[f1,gp[[10]],gp[[14]],.2,x,y];
L10=pathFinder[f1,gp[[14]],gp[[9]],.2,x,y];]

Show[RegionPlot[f>0 && x^2+y^2<25,{x,-5,5},{y,-5,5},BoundaryStyle->None], ListLinePlot[{Join[L1,L2],Join[L3,L4],Join[L5,L6],Join[L7,L8],Join[L9,L10]}, PlotStyle->{Black,Red,Brown,Blue,Purple}]]

We got a good plot and decomposition, but it took a long time.

For singular curves, it is important to note that pathFinder can trace into a singular point but not out of one. Thus we should augment our critical points with points a little upstream from the singular points and trace into the critical point by replacing the curve equation with its negative.

h2= -(x^4-x^2y+y^3) ;
cph2=Chop[gPoints[h2,.5,x,y]];
cph2a=Join[cph2,{x,y}/.NSolve[{h2,x^2+y^2-.02},{x,y},Reals]]
plta:=gaussPlot[h2,.5,cph2a];
pltc=Graphics[Table[{Black,Arrow[{p,p+.1gTvec[h2,p,x,y]}]},{p,Delete[cph2a,5]}]];
Show[plta,pltc]

H1=pathFinder[h2,cph2a[[2]],cph2a[[6]],.05,x,y];
H2=pathFinder[h2,cph2a[[6]],cph2a[[5]],.025,x,y];
H3=pathFinder[-h2,cph2a[[7]],cph2a[[5]],.025,x,y];
H4=pathFinder[h2,cph2a[[7]],cph2a[[3]],.025,x,y];
H5=pathFinder[h2,cph2a[[3]],cph2a[[8]],.025,x,y];
H6=pathFinder[h2,cph2a[[8]],cph2a[[5]],.025,x,y];
H7=pathFinder[-h2,cph2a[[9]],cph2a[[5]],.025,x,y];H8=pathFinder[h2,cph2a[[9]],cph2a[[4]],.02,x,y];
H9=pathFinder[h2,cph2a[[4]],cph2a[[11]],.025,x,y];
H10=pathFinder[h2,cph2a[[11]],cph2a[[5]],.025,x,y];
H11=pathFinder[-h2,cph2a[[10]],cph2a[[5]],.025,x,y];
H12=pathFinder[h2,cph2a[[10]],cph2a[[1]],.05,x,y];

{Show[gaussPlot[h2,.5,cph2a],ListLinePlot[{H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12}]],ListLinePlot[{H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12},PlotStyle->Blue, Axes->False, Frame->True]}

Length[Join[H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12]]

This takes some trial and error, including using different step sizes and some experimentation to generate these 113 points, but this gives the best plot of the bow curve.

4.5 The Affine Graph of a Curve

I have tried to use the word plot for a picture of a curve in this book. The word graph will be used in the same sense as it is used in the Graphs & Networks section of the Wolfram Language & System Documentation Center. Graphs provide a combinatorial way to view the topology of a curve in a simple way. We will build curve graphs using critical points as vertices and edges where there is a segment directly connecting one critical point to another. We will see that curve graphs have a special property that allows us to decompose curve graphs analogously to how we decompose curves.

CGdata={2->10, 10->3, 3->6, 6->7, 7->10, 10-> 4, 4->5, 5->9,9->6,6->8,8->5,5->1};
CG=Graph[CGdata, VertexLabels->"Name"]

CGlist=EdgeList[CG]

This last format is preferred because it distinguishes between directed edges and undirected edges.

e1=DirectedEdge[1,2]
u1=UndirectedEdge[1,2]
UG=UndirectedGraph[CG]
EdgeList[UG]

{VertexDegree[CG,2], VertexDegree[CG,4],VertexDegree[CG,10]}

Gaussian decomposition shows a very important fact about curve graphs. We call vertices not on the boundary of our disk interior vertices.

Every interior vertex has an even degree.

Like curves, graphs have connected components.

WeaklyConnectedComponents[CG]

c2=y^2-x^3+x;
gpc2=gPoints[c2,3,x,y];
{gaussPlot[c2,3,gpc2], GC2=Graph[{4->5,5->4,1->3,3->2},VertexLabels->"Name"]}

WeaklyConnectedComponents[GC2]

We will see the theoretical importance of a curve graph in Chapter 5. One immediate use is to feed information on a graph to our pathFinder algorithm to automatically find a piecewise linear approximation to a curve component. For the next procedure defined by plotPath, we will take a curve f=0 and a set of points (that is, vertices V), including enough boundary and critical points to identify segments for pathFinder.

plotPath[f_,G_,V_,x_,y_,s_]:=Module[{n,L},
n=Length[G];
L=Flatten[Table[pathFinder[f,V[[G[[i,1]]]],V[[G[[i,2]]]],s,x,y],{i,n}],1];
Append[L,V[[G[[n,2]]]]]]

V=gpc2

G1 = {1->3,3->2};
G2={4->5,5->4};

L1=plotPath[c2,G1,V,x,y,.1];
L2=plotPath[c2,G2, V,x,y,.1];

Show[RegionPlot[c2>0 && x^2+y^2<9,{x,-3,3},{y,-3,3},BoundaryStyle->{None,Pink}],ListLinePlot[{L1,L2}, PlotStyle->{Black,Orange}],ListPlot[Table[Labeled[V[[n]],n],{n,5}],PlotStyle->Red, PlotMarkers->Automatic]]

f=newtonHyperbola[376,x,y]
cp376=Join[criticalPoints[f,x,y],{x,y}/.NSolve[{f,-x+y-20},{x,y},Reals]];

cp376={{-5.28009203049686`,-17.206102954932472`},{11.659821180891406`,3.1684250826109506`},{5.167252947836531`,0.5278467102679588`},{1.4273504023589225`,2.0419148559610174`},{-0.39812614887988745`,-0.3862418222579013`},{-19.645226575920613`,0.3547734240793867`},{0.05173064096128144`,20.05173064096128`}}

gaussPlot[f,20,cp376]

G={7->4,4->3,3->5,5->6};
L=plotPath[f,G,cp376,x,y,.8];
ListLinePlot[L, Axes->False,Frame->True]

Note that this example will reappear in the next section.

4.6 Curve Graphs and Plots with a Differential Equation Solver

Ideally, we could develop a black box algorithm to find a graph of an affine algebraic curve and then plot it like we did in Section 4.4, but as was mentioned in Chapter 2, this is difficult when working numerically. Moreover, in Chapters 8 and 9, we will be more interested in plotting specific components rather than the whole curve, so we will go partway, letting several procedures do much of the work, but they will still require some direction from the user. Our main goal in this section is to give a randomized black box procedure, which will suggest a graph. The user can accept the result or try again. (The reader who is not interested in the details should jump to the summary at the end of this section.)

Our first step is to introduce a faster plotter than our naive plotter, which will allow smaller step sizes and give more chance of our software finding the next critical point. We have already seen gTvec and can now give our basic routine deTrace using the Wolfram Language’s NDSolveValue.

Here, f is the curve equation, p is the initial point, s is the maximum Euclidean distance from the initial point that can be traced and r is the diameter of a circle around the origin containing the bounded part of the curve to be traced. We also note that there is a randomization, so running this procedure several times with the same input will give different outputs.

deTrace[f_,p_,x_,y_,s_,r_]:=Module[{X,Y,F,u,t,tv,ttv,XS,YS},
F=f/.Thread[{x,y}->{X[t],Y[t]}];
tv=gTvec[f,p,x,y];
u=Normalize[RandomReal[{.8,1.2},2]];
ttv=u.tv;
NDSolveValue[{u.{X'[t],Y'[t]}==ttv,D[F,t]==0, {X[0],Y[0]}==p,
WhenEvent[X[t]^2+Y[t]^2>r^2||u[[1]]*X'[t]+u[[2]]*Y'[t]==0||(X[t]-p[[1]])^2+(Y[t]-p[[2]])^2>s^2,"StopIntegration"]},{X,Y},{t,0,20}]
]

The disadvantage to this, compared to pathFinder, is that it is hard to tell NDSolveValue exactly where to stop. As we will see, this method does not like to go around sharp curves and may stop on its own. Also, it will stop when its machine-defined parameter hits 20. On the other hand, deTrace is much faster than pathFinder and, therefore, can use a machine-defined smaller step size and a higher order of interpolation.

f=newtonHyperbola[376,x,y];
p=cp376[[7]]
{X,Y}=deTrace[f,p,x,y,20,20]

dom=deEnd[{X,Y}][[2]]

ParametricPlot[{X[t],Y[t]},{t,0,dom}, PlotRange->{{-20,6},{-5,20}},Axes->False,Frame->True]

Here we barely got started around the bend when it stopped, so it is useful to have a procedure to tell us where and when (with t as parameter) these charts end.

deEnd[{X_,Y_}]:=Module[{t1},
t1=X[[1,1,2]];
{{X[t1],Y[t1]},t1}]

deEnd[{X,Y}]

Norm[cp376[[6]]-cp376[[7]]]

deEnd[deTrace[f,p,x,y,30,20]]
deEnd[deTrace[f,p,x,y,30,20]]
deEnd[deTrace[f,p,x,y,30,20]]

We don’t get anywhere near our desired endpoint in any of these examples. In fact, we don’t get to negative x at all.

p=cp376[[7]]
{X1,Y1}=deTrace[f,p,x,y,30,20];
e1=deEnd[{X1,Y1}]
{X2,Y2}=deTrace[f,e1[[1]],x,y,30,20];
e2=deEnd[{X2,Y2}]
{X3,Y3}=deTrace[f,e2[[1]],x,y,30,20];
e3=deEnd[{X3,Y3}]
{X4,Y4}=deTrace[f,e3[[1]],x,y,30,20];
e4=deEnd[{X4,Y4}]
{X5,Y5}=deTrace[f,e4[[1]],x,y,30,20];
e5=deEnd[{X5,Y5}]
{X6,Y6}=deTrace[f,e5[[1]],x,y,30,20];
e6=deEnd[{X6,Y6}]
cp376[[6]]

My experiment took six applications to get to the desired endpoint. The reader with Mathematica can rerun the experiment, but often it will take even more applications. However, the user will see that this experiment is still much faster than running plotPath once!

Unfortunately, InterpolatingFunction does not combine nicely. Since the domains are different, they cannot be used in the same ParametricPlot statement. You can construct Piecewise functions out of them, but this is also messy.

ListLinePlot[ {p,e1[[1]],e2[[1]],e3[[1]],e4[[1]],e5[[1]],e6[[1]]},Axes->False,Frame->True,PlotRange->All]

pathFinderDE[f_,p_,q_,r_,x_,y_,eps_]:=Module[{i,e,X,Y,s},
Off[NDSolveValue::ndsz];
e=p;
s=Norm[p-q];
i=1;
Reap[Sow[p];
While[s>eps && Norm[e]<2*r&& i<60,
{X,Y}=deTrace[f,e,x,y,s*eps,r];
e=deEnd[{X,Y}][[1]];
Sow[e];
s=Norm[e-q];
i++];
Sow[q]][[2,1]]];

f=newtonHyperbola[376,x,y];
p=cp376[[7]]
q=cp376[[6]]
G={7->4,4->3,3->5,5->6}
{td,L1}=Timing[pathFinderDE[f,p,q,20,x,y,.05]];
{tp,L2}=Timing[plotPath[f,G,cp376,x,y,.5]];
{td,tp}
{Length[L1],Length[L2]}
ListLinePlot[{L1,L2},PlotStyle->{Red,Blue},Frame->True,Axes->False, PlotRange->All]

Note that plotPath took about 10 times longer with about 1.5 times the number of points, but the resulting plot is almost the same (see the difference marked in red).

plotPathDE[f_,G_,V_,x_,y_,s_]:=Module[{r,n,L,LL,PF},
n=Length[G];
r=Max[Norm/@V];
L=Flatten[Table[pathFinderDE[f,V[[G[[i,1]]]],V[[G[[i,2]]]],r,x,y,s],{i,n}],1];
Append[L,V[[G[[n,2]]]]]]

{tdp,L3}=Timing[plotPathDE[NH376,G,cp376,x,y,.25]];
tdp
Length[L3]
ListLinePlot[{L1,L3},PlotStyle->{Red,Blue},Frame->True,Axes->False, PlotRange->All]

The difference here, marked in red, is much smaller. The next step is to wean pathFinderDE from needing the endpoint. The preferred state would be pathFinderDE having the choice of endpoints—from a critical point set CP to being able to find the next endpoint in the Gauss direction of the curve, but for nonsingular curves, this is just a bookkeeping problem. We replace q by the critical point set CP and s by the minimum distance to a point in CP. We do have one technical problem: we first need to delete the starting point from the critical point set, but in Mathematica, using Drop and Delete on lists requires a position, so we will use a position instead. However, the Wolfram Language Position function does not work well with machine numbers, so the following procedure finds the number of the point in CP following the (n^th) point in CP.

endFinder[f_,n_,CP_,x_,y_,eps_]:=Module[{r,i,k,e,X,Y,s,CPD},
Off[NDSolveValue::ndsz];
r=Max[Norm/@CP];
e=CP[[n]];
CPD=CP;
CPD[[n]]=Nothing;
s=Norm[e-First[Nearest[CPD,e]]];
i=1;
e=Catch[
While[s>eps && Norm[e]<1.01*r&& i<60,
{X,Y}=deTrace[f,e,x,y,s,r];
e=deEnd[{X,Y}][[1]];
If[Norm[e]> r, Throw[e]];
s=Norm[e-First[Nearest[CPD,e]]];
i++];
Throw[e]];
If[Norm[e]>r,Return[0]];
k=1;
While[k≤Length[CP] && Norm[e-CP[[k]]]>s,k++];
If[k<=Length[CP],k,Fail]]

NH376=newtonHyperbola[376,x,y];
gaussPlot[NH376,21,cp376]

{endFinder[NH376,7,cp376,x,y,.2],
endFinder[NH376,4,cp376,x,y,.2],
endFinder[NH376,3,cp376,x,y,.2],
endFinder[NH376,5,cp376,x,y,.2]}

gaussGraphNS[f_,CP_,x_,y_,eps_]:=Module[{n,G},
n=Length[CP];
G=Select[Table[k->EndFinder[f,k,CP,x,y,eps],
{k,n}],#[[2]]>0&];
graph[G,VertexLabels->"Name",ImageSize->Small]];

gaussGraphNS[NH376,cp376,x,y,.2]

g=newtonHyperbola[425,x,y];
ggp=gPoints[g,10,x,y];
{gaussGraphNS[g,ggp,x,y,.2],
Show[ContourPlot[g==0,{x,-10,12},{y,-10,12}],ListPlot[Table[Labeled[ggp[[n]],n],{n,Length[ggp]}]]]}

n=26536;
f=newtonHyperbola[n,x,y]
CP1=gPoints[f,15,x,y]
gaussGraphNS[f,CP1,x,y,.4]
r=1.1*Max[Norm[#]&/@CP1];

P1=plotPathDE[f,{2->10,10->3},CP1, x,y,.1];
P2=plotPathDE[f,{4->9,9->5},CP1,x,y,.1];
P3=plotPathDE[f,{6->12,12->7},CP1,x,y,.1];
P4=plotPathDE[f,{8->11,11->1},CP1,x,y,.1];
{Show[RegionPlot[f>0, {x,-r,r},{y,-r,r}],ListPlot[Table[Labeled[CP1[[n]],n],{n,Length[CP1]}]]],
Show[RegionPlot[f>0, {x,-r,r},{y,-r,r},BoundaryStyle->None],
ListLinePlot[{P1,P2,P3,P4}, Axes->False,Frame->True, PlotRange->All],
ListPlot[Table[Labeled[CP1[[n]],n],{n,Length[CP1]}]]]}

In this example, the graph out-performs the plot!

If there is a singular (or nearly singular) point, then gaussGraphNS does a decent job of tracing into the point.

CP2={{0,0},{1,0},{0,1},{-1,0},{0,-1}};
gaussGraphNS[x y,CP2,x,y,.4]

To rectify this situation, we use a weak criterion for singularity, Norm[gTvec[f,p,x,y]]<eps/2 and do not allow this point as a starting point.

nearPoints[f_,r_,x_,y_,pt_,ep_]:=Module[{np,stnpts},
np={x,y}/.NSolve[{f,(x-pt[[1]])^2+(y-pt[[2]])^2-ep^2},{x,y},Reals];
np=Select[np,Norm[#-pt+ep/2*Normalize[gTvec[f,#,x,y]]]> ep&];
Select[np,Norm[#]<r&]];

For each near point, we add an edge to the graph from the singular point to the near point. This leads to our main (almost black box) routine to find the graph of a curve on a bounded disc. It assumes that a set of critical points (including boundary and singular points) has been calculated.

gaussGraph[f_,gp_,x_,y_,eps_]:=Module[{n,k,e,G,np,CP,r,S},
CP=gp;
r=Max[Norm[#]&/@gp];
n=Length[CP];
G={};
k=1;
S=Reap[While[k≤n, 
If[Norm[gTvec[f,CP[[k]],x,y]]<eps/2, 
np=nearPoints[f,r,x,y,CP[[k]],1.25*eps];
Do[Sow[k->n+j],{j,Length[np]}];
CP=Join[CP,np];
n=Length[CP];
e=0,
e=endFinder[f,k,CP,x,y,eps]];
If[e >0,Sow[k->e]];
k++]];
G=S[[2,1]];
Echo[{Show[ContourPlot[f==0, {x,-1.2r,1.2r},{y,-1.2r,1.2r}],ListPlot[Table[Labeled[CP[[j]],j],{j,Length[CP]}]]],
Graph[G,VertexLabels->"Name"]}];
{CP,G}];

f=y^2-x^3-x^2;
gp = gPoints[f,3,x,y]
gaussGraph[f,gp,x,y,.1]

ch4=1-8x^2+8x^4+3y-4y^3;
gp=gPoints[ch4,2,x,y];
gaussGraph[ch4,gp,x,y,.1]

Show[ContourPlot[ch4==0, {x,-2,2},{y,-2,2}],ListPlot[Table[Labeled[gp[[n]],n],{n,Length[gp]}]]]

gp1=Delete[Delete[gp,10],6]
{ch4pts,ch4G}=gaussGraph[ch4,gp1,x,y,.1]

Due to the complexity of this software, it is possible that the labels may change when this is run on different sessions, versions of Mathematica and machines.

Ch4A=\[LeftAssociation]1->{-1.2878893413505916`,1.5301441253802013`},2->{1.2878893413505876`,1.5301441253802006`},3->{-0.8672290665235644`,-0.938587657300242`},4->{0.`,1.0000000000004285`},5->{0.7071067811865476`,0.5000000000000014`},6->{0.`,-0.49999999999999856`},7->{0.8672290665235052`,-0.9385876573001991`},8->{-0.2504274405673725`,-0.18503072151060296`},9->{-0.707106781186551`,0.5000000000000037`},10->{0.2504274405674338`,-0.18503072151063518`},11->{0.7718939492398669`,0.6069000601292057`},12->{0.6413097819106776`,0.3937185110835575`},13->{0.08339503632735218`,-0.5931142734276572`},14->{-0.08043678659376786`,-0.4043186362844401`},15->{-0.6378263144263652`,0.604044302705585`},16->{-0.7686094142232938`,0.391177088214161`}\[RightAssociation];

G1={1->9,9->16,16->3,3->6,6->13,13->7,7->5,5->11,11->2};
G2={9->15,15->4,4->5,5->12,12->10,10->6,6->14,14->8,8->9};
{Labeled[Graph[G1,VertexLabels->"Name"],"G1"],Labeled[Graph[G2,VertexLabels->"Name"],"G2"]}

These subgraphs traverse the boundaries of our two connected components of the complement of the curve ch4=0. This is the decomposition of the curve proposed by Gauss.

In Chapters 5 and 8, we will learn how to extend these bounded graphs to global graphs. For now, we finish by extending our plotPath algorithms to this singular case.

plotPathS[f_,G_,V_,x_,y_,s_]:=Module[{n, L,r},
n=Length[G];
r=Max[Norm/@V];
L=Flatten[Table[If[Norm[gTvec[f,V[[G[[i,1]]]],x,y]]<s/2, {V[[G[[i,1]]]],V[[G[[i,2]]]]},
pathFinder[f,V[[G[[i,1]]]],V[[G[[i,2]]]],s,x,y]],
{i,n}],1];
AppendTo[L,V[[G[[n,2]]]]];
L]

plotPathDES[f_,G_,V_,x_,y_,eps_]:=Module[{r,n, L},
n=Length[G];
r=Max[Norm[#]&/@V];
L=Flatten[Table[If[Norm[gTvec[f,V[[G[[i,1]]]],x,y]]<eps/2, {V[[G[[i,1]]]],V[[G[[i,2]]]]},
pathFinderDE[f,V[[G[[i,1]]]],V[[G[[i,2]]]],r,x,y,eps]],
{i,n}],1];
AppendTo[L,V[[G[[n,2]]]]];
L]

f=x^4-4x y+y^4;
gp=Chop[gPoints[f,2.5,x,y]]

{CP,GG}=gaussGraph[f,gp,x,y,.05]

P=FindEulerianCycle[GG][[1]]

{ListLinePlot[plotPathS[f,P,CP,x,y,.1],Axes->False,Frame->True],
ListLinePlot[plotPathDES[f,P,CP,x,y,.2],Axes->False,Frame->True]}

plotPathS produces a better curve, but plotPathDES runs more quickly. Note that the linear approximation near the singular point is not the problem in the right image.

cp370={{-1.28789`,1.53014`},{1.28789`,1.53014`},{0.867229`,-0.938588`},{0.`,1.`},{-0.707107`,0.5`},{0.`,-0.5`},{0.250427`,-0.185031`},{-0.867229`,-0.938588`},{-0.250427`,-0.185031`},{0.707107`,0.5`},{-0.637826`,0.604044`},{-0.768609`,0.391177`},{0.083395`,-0.593114`},{-0.0804368`,-0.404319`},{0.771894`,0.6069`},{0.64131`,0.393719`}};

G1={1->9,9->16,16->3,3->6,6->13,13->7,7->5,5->11,11->2};
G2={9->15,15->4,4->5,5->12,12->10,10->6,6->14,14->8,8->9};
Show[RegionPlot[ch4>0,{x,-2,2},{y,-2,2}, BoundaryStyle->None, PlotStyle->LightOrange],ListLinePlot[plotPathS[ch4,G1,Values[Ch4A],x,y,.1],PlotStyle->Blue],ListLinePlot[plotPathS[ch4,G2,Values[Ch4A],x,y,.1],PlotStyle->Green],Axes->False, Frame->True]

4.7 Paradigm for Analyzing Curves in a Bounded Circle about the Origin

In the social sciences, which are not as precise as mathematics, the analog of an algorithm is a paradigm. Since analyzing a curve numerically is not a precise science either, we will summarize our approach by a paradigm.

First, we choose a radius r and find boundary and critical points using, for example, gPoints.

For now, the region should be a finite disk centered at the origin. Later in Chapter 8, we show how to use a more general, possibly infinite, region bounded by a conic. It is important that our points contain all the points on the boundary of the region and that all singular points are present. Also, each connected component of the curve, as well as each algebraic component, should have at least two points represented. gPoints will accomplish this but may give too many or too few points.

Next, we find the curve graph, in this case using gaussGraph. Separate into weak components.

If it is a simple curve, gaussGraph should easily find the graph, but note that this algorithm is sensitive to the choice of r. If incorrect results are returned, try a larger or smaller r. Deleting unnecessary points or adding additional points to the critical point set may help. As a last resort, a RegionPlot with critical points labeled can help you give a graph manually. If the curve has more than one algebraic component, manual graph finding may be the best option—better even than algebraically decomposing the curve. The topology of the curve is easiest found from the picture given by Mathematica, but unfortunately, the decomposition tools do not work as well as preserving the paths and directions do, in the picture.

Optionally, plot using plotPathS or plotPathDES.

It is important to have actual paths—that is, the directed graph should be given by edges with the initial point of one edge being the terminal point for the preceding edge. In particular, EdgeList and EdgeRules are not always in path order. In theory, each path should begin and end at border points or should be a closed path. In the first case, use FindPath to convert EdgeList or EdgeRules to a path starting and ending in a boundary point. In the second case, you can use FindEulerianCycle to get a path. Paths should always be in the Gauss direction—clockwise about the positive set of the equation. Again (as indicated in in the previous step), the picture, if correct, gives the decomposition into paths we want. For a relatively small critical point set, manually working from the picture is often the easiest strategy. However, if there are singular points, make sure that the paths go around connected components of the positive set f>0.

f=gaussCurve[z^5-5z^4+9z^3-5z^2,x,y,z]

gp=Chop[gPoints[f,4,x,y]]

{cp,G}=gaussGraph[f,gp,x,y,.2]

gOrigA=\[LeftAssociation]1->{-3.710503137713464`,1.4940436623534898`},2->{-1.694287437347172`,3.6234500244446064`},3->{1.0256836809227585`,3.8662608534200804`},4->{2.993193894204237`,2.653448758069135`},5->{3.885912999933733`,0.9485147109801363`},6->{3.885912999933767`,-0.9485147109801597`},7->{2.993193894204127`,-2.6534487580692407`},8->{1.0256836809227194`,-3.866260853420112`},9->{-1.6942874373471617`,-3.623450024444607`},10->{-3.7105031377133324`,-1.4940436623535953`},11->{2.015235448916427`,0.8415109240162972`},12->{2.015235448916427`,-0.8415109240162972`},13->{1.`,0},14->{0,0},15->{0.13687218865500828`,-0.209203259947802`},16->{-0.19697325992704043`,0.15395302814077477`}\[RightAssociation];

Show[RegionPlot[f>0,{x,-4,4},{y,-4,4}],ListPlot[Table[Labeled[gOrigA[i],i],{i,16}]]]

P1={4->11,11->5};
P2={6->12,12->7};
P3={8->13,13->3};
P4={2->14,14->16,16->1};
P5={10->14,14->15,15->9};

PP1=plotPathS[f,P1,cp,x,y,.2];
PP2=plotPathS[f,P2,cp,x,y,.2];
PP3=plotPathS[f,P3,cp,x,y,.2];
PP4=plotPathS[f,P4,cp,x,y,.2];
PP5=plotPathS[f,P5,cp,x,y,.2];

Labeled[Show[RegionPlot[{f>0 && x^2+y^2<16},{x,-4,4},{y,-4,4}, BoundaryStyle->None],ContourPlot[x^2+y^2==16,{x,-4,4},{y,-4,4},ContourStyle->Pink],
ListLinePlot[{PP1,PP2,PP3,PP4,PP5}],Axes->False,Frame->True], "Components of gaussCurve[(z^5)-4(z^4)+9(z^3)-5(z^2)]"]