Updated 10/23/18 Barry Dayton

Classical vs. Numerical Perspective

The split between exact and numerical mathematics dates back 3,500 years. As urbanization began in Egypt and Mesopotamia, there was a need for mass production of essential commodities like bread and beer. Proportional reasoning was needed to get the recipes right. For more in-depth information, see George Gheverghese Joseph’s book The Crest of the Peacock (1991). The Egyptians solved this problem using very awkward unit fractions, which persisted in Europe through the Early Middle Ages. In Mesopotamia, a base-60 number system was easily modified by adding place values for (1/60),(1/3600),... etc. They were aware that a quotient, such as 1÷7, did not have a finite representation, but the Mesopotamians used their numbers for practical things like manufacturing bread and beer and building houses and roads; measurement was not exact anyway, so a sexagesimal place or two was sufficient.

It was Pythagoras, ironically a shaman, who saw in numbers a way to demonstrate nonphysical truths. He took the well-known Mesopotamian right triangle theorem (a^2)+(b^2)= (c^2) used in engineering and made it into a Diophantine problem by searching for integer solutions. For an accurate account of the life and work of Pythagoras, read Walter Burkert’s book Lore and Science in Ancient Pythagoreanism (1972). Pythagoras, who was no doubt aware of Mesopotamian mathematical proficiency, realized that it was the Egyptian numbers that were true.

Plato and Aristotle, who had little good to say about Pythagoras, liked this idea that mathematics could be repurposed to illustrate truth by logic and that Mesopotamian numbers were not acceptable for this purpose. Classical mathematics thus began with Euclid, a follower of Plato who eliminated all numbers by using only geometrical concepts and Aristotelian logic.

To be fair, while Mesopotamian mathematics produced many good algorithms, iterative methods for solving quadratic and cubic equations that had noninteger solutions, it is the classical method that allowed for proof and the development of abstract mathematics. The two views of mathematics have persisted to this day; often we call them pure and applied, but that is not an accurate distinction. It is more a matter of whether we need an actual answer that might not be a rational number or we just want an existence or structural theorem. In the journals of the applied mathematics world, like SIAM or ACM, you will rarely come across a decimal number.

The mathematicians who have been the most influential are often those who are able to work in both worlds, such as Archimedes, Newton, and Gauss. You will see the latter two names often in this book. Keep in mind as you read forward that often numerical methods are available only because we have the classical insight to go on.

As you will see, it is the numerical approach that leads to many simple but effective algorithms. In this book, our goal is to write simple code of a few lines each, and the numerical approach will allow us to do that. However, if your goal is to develop black box algorithms, algorithms that answer difficult situations with no human input, the classical method is preferable, but expect thousands of lines of code rather than a few. The longest code in this book is the non-numerical black box code for diamond diagrams (dDiagram), which, with constants and subroutines, takes about four notebook pages, and even then, desired options are omitted.

The reader who is knowledgeable about numerical methods and the Wolfram Language may skim the rest of this chapter, but there may be some things that are slightly non-standard that should be noted. In particular, I explain what is meant by the solution of f [x,y]=0 in the numerical case.

2.1 Numbers in the Wolfram Language

One reason I like the Wolfram Language is that it treats numeric data in a sophisticated way but takes care of the details, so for the most part, the user does not need to pay attention to the details of how it is done. The documentation is quite complete, so I will just give an overview.

Numbers can be exact, symbolic, or numeric. Exact numbers are integers or rational numbers, and arithmetic, addition, subtraction, multiplication, division, and exponentiation with an integer exponent all preserve exactness.

34+2/3

(3/17-7/19)^7

Symbolic numbers are given by a predefined symbol or letters e or E, ℼ or pi or an expression using exact or symbolic numbers, such as Sqrt[2] or Sqrt[e^4].

a=Log[E^7]

Sqrt[E^4]

sol1=x/.Solve[x^2-Pi x^2+E x-7==0,x][[1]]

The output sol1 is now a symbolic number.

The Wolfram Language also has approximate decimal numbers of two types, machine numbers and arbitrary-precision numbers. In this book, we will be satisfied with using machine numbers, which are faster and accurate enough for our needs. However, the reader who wants to extend our results to curves above degree 6 may need to consider using arbitrary-precision numbers and rewriting much of the code.

Any number entered with a decimal point is considered, by default, to be a machine number. The precision of a machine number is the constant MachinePrecision, which may be dependent on your hardware.

$MachinePrecision

For technical reasons, a number may be 53 base-2 bits, which is the approximate number of decimal digits available, so what that means is that we can expect about 16-digit precision.

The user must be particularly aware of the difference between the display number and the actual machine number it represents. The display number typically uses only four or five significant digits (leading or trailing zeros are not considered significant), and it is rounded. Internally, Mathematica stores it as a full machine-precision number. If you put your cursor at the end of the numerical output and evaluate the cell, you will get the full precision output.

a=Sqrt[2.0]

One of the basic rules for using numerical numbers is never manually type in a display output that is being used in a calculation, including in other software. This is part of the more general rule never round in the middle of a calculation.

5.12*4*Pi

{N[1/7],N[sol1]}

A machine number should never be converted back to an exact number in the context of this book. In particular, never use the function Rationalize in our context.

Remember that equals is not a concept in numerical mathematics, but the relations < and > are allowed. In particular, machine numbers are never equal to zero, even if they display as such. In certain cases, an algorithm is designed to return zero or perhaps a zero coefficient or a zero imaginary part of a complex number. In this case, judicious use of Chop is acceptable.

InverseFourier[Fourier[{20,10,5,.1,0}]]

Chop[%]

On the other hand, the central point of this book is understanding the solution of a bivariate polynomial equation, even when the polynomial has numerical coefficients. Exact zero is not a machine number. What we mean by f [x,y]=0 is a central issue in this book and will be discussed in the next section.

Certain Wolfram Language built-in functions work fine with numerical data, such as CoefficientRules, Expand, NSolve, and FindRoot. If the bivariate polynomial is exact, like in the previous example, one could use Factor, but it can’t be used with numerical polynomials. In general, trying to factor numerical polynomials is tricky, but in principle, numerical polynomials do not factor, so neither the Wolfram Language nor I give a numerical factoring procedure.

Clear[a,x,y];
Factor[2 (x^3) y-2 (a^2) x y-3 (a^2) (x^2)+3 (a^4)]

Factor[-2+(x^2)+2 x y+(y^2)]

a=(x+y+N[Sqrt[2]])
b=(x+y-N[Sqrt[2]])

Expand[a b]

In Chapter 1, I do have a reliable procedure to decide if a numerical polynomial is reduced (no repeated factors) or to attempt to find greatest common divisors. Built-in functions such as Simplify or Reduce should not be used on numerical polynomials. In fact, most of the polynomial algebra functions in the Wolfram Language should not be used on numerical polynomials. Mathematica is perhaps the only computer algebra software where computations like GroebnerBasis using numerical polynomials are allowed and do, in fact, work reasonably well, but I avoid these methods in this book. Some numerical methods I like for working with polynomial systems are discussed in Appendix 1.

2.2 Iteration and Newton’s Method

We often need to decide if a machine number is a zero of a polynomial, but since zero is not a machine number, evaluating the polynomial at that number may display zero but will not give zero.

f=x^4+2x^3-x-.2;
x0=-0.7813887111276197`;
f/.{x->x0}

The value returned is called the residue. A large residue would indicate that x0 is not a solution, but there is no general rule as to when a residue is small enough to be a solution.

A better way to identify a zero of a polynomial equation is to produce that zero by an iterative solution method.

Specifically, we will illustrate using Newton’s method, which is both simple and effective at identifying solutions to numerical equations. Actual software may use somewhat different iterative methods. Historically, what we now call Newton’s method was actually discovered by Thomas Simpson; Newton himself used a method dating back to Mesopotamia. Simpson named it after his mathematical hero, Newton. The argument as to whether Newton or Raphson is the discoverer of the modern Newton method is moot since neither used this method. Prior to the use of electronic computers, Newton’s method was frowned upon partly because it was not exact and thus not suitable for mathematics but also because it works best with complex initial points, and mathematicians did not have the confidence that Newton’s method would work to grind through the hand calculations of the complex division in sufficient precision.

In the late 40s, Uspensky wrote about Graeffe’s method, which is another numerical but not iterative method using logarithms, which “is the only practical method for [the] computation of imaginary roots” (Uspensky, 1948). With computers, it is now easy to do a lot of experiments to convince ourselves that the complex Newton’s method almost always works, but a theoretical analysis is not good since, in theory, there are, for degree 3 and higher, an infinite number of initial points where Newton’s method does not work (see Section 4.1).

Newton’s method in one variable works as follows. Start with a polynomial f [z]=(Subscript[a, n])(z^n)+(Subscript[a, n-1])(z^(n-1))+...+ (Subscript[a, 1])z+(Subscript[a, 0]), using the variable z because we may want to work with complex numbers even though the (Subscript[a, j]) are real. Pick a random initial value (Subscript[p, 0]), preferably complex, and then iteratively define (Subscript[p, j+1])=(Subscript[p, j])-((f [(Subscript[p, j])])/(f' [(Subscript[p, j])])) where f'[p] is the derivative of f [z] at (Subscript[p, j]).

newton1a[f_,p0_,n_,z_]:=Module[{p},
p=N[p0];
Do[Print[p]; p=p-(f/D[f,z])/.{z->p},{i,n}];
p]

f=-1-z-0.2` (z^2)+0.008` (z^3)+0.000064` (z^4)

p0=RandomComplex[{-10-10I,10+10I}]

newton1a[f,p0,10,z]

As noted earlier, in the modern use of Newton’s method, an important point is that the stopping criteria is not that the residue f [p] is small but that the change in p is small.

newton1[f_,p0_,z_]:=Module[{p1,p2},
p1=N[p0];
Chop[Catch[Do[p2=p1-(f/D[f,z])/.{z->p1};
If[Abs[p1-p2]<1.*^-14,Throw[p2]];
p1=p2,{30}]]]
]

p0=RandomComplex[{-10-10I,10+10I}]

z1=newton1[f,p0,z]

NumberForm[z1,16]

The polynomial we used is numerical since at least one coefficient has a decimal point, and even though these decimals are equivalent to rational numbers, they are, by default, treated as machine numbers.

f/.{z->z1}

If there is a question whether a result of Newton’s method with a complex starting point is real, one can take the real part of the answer to one of the Newton programs and run it from that starting point. If the procedure comes close to obtaining the same result, then it is likely a real number since a Newton function with a real function and a real initial point stays real.

z1=newton1a[f,p0,8,z]

z2=Re[z1]

z3=newton1a[f,z2,5,z]

We can conclude this is a real number.

In practice, rather than using either of the Newton functions I created and used perviously, one would use the built-in function FindRoot[f,{z,z0}], which is an implementation of Newton’s method.

NSolve[f,z,Reals]

In two variables, which will mostly be our case, one needs two equations and an initial point that has two coordinates. The two-variable version of the Newton method is often called the Gauss–Newton method because sophisticated implementations use ideas of both Gauss and Newton. However, we will use the built-in function FindRoot.

f=-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;
g=x^2+y^2-64;

{p1,p2}=RandomComplex[{-10-10I,10+10I},2]

q={x,y}/.FindRoot[{f,g},{{x,p1},{y,p2}}]

q0=Re[q]

FindRoot[{f,g},{{x,q0[[1]]},{y,q0[[2]]}}]

NSolve[{f,g},{x,y},Reals]

I can now address the main issue: what do I mean by a machine number pair, being a solution to f=0?  A working criteria is that we will accept the output of  NSolve or FindRoot as in the previous examples as a solution so long as this point is not a singular point.  Further functions in this book which directly report such output will also give solutions to f=0.

Singular points will be defined in Chapter 3 and may require more analysis.  Especially if the polynomial f itself has machine number coefficients or if the singular point is not given by rational or symbolic numbers Newton’s method, and the implimentations NSolve and FindRoot may not give accurate results.  See Appendix 1 and 2 for examples and discussion.  Also my code in Chapter 6 for complexProjectiveSingularPoints may not give accurate enough estimates of infinite singular points for a given application.  Of course, if the polynomial and calculated singular point are given with rational or symbolic numbers and evaluation gives an exact zero then we will accept that singular point as a solution.

In Chapter 4, we will give Gauss’s alternate definition of a point that is on a possibly numerical curve f=0 where p is a topological boundary point of the region f>0. We will then show how to find enough points to approximately plot the curve.

2.3 Randomness in Numerical Algorithms

An important way in which numerical algorithms are simpler than corresponding exact algorithms is the ability to use randomness. In many cases, an algorithm will work most of the time, but there are special cases in which it will not work. In Newton’s method, for example, we noted that most initial points gave rise to a sequence of values converging to a zero, but we also mentioned that for univariate polynomials of degree 3 or more, there are many initial points where Newton’s method does not, in theory, converge. It is not, however, possible to give a useful description of these points to avoid, especially before finding the solutions, but in most cases, these points are surrounded by points that are good initial points. For this reason, we choose random points as initial values and, particularly if we choose a random complex number, the probability that we choose a bad point is essentially zero.

f= 3x^4-6x^3-x^2+4x-4

newton1a[f,1,6,x]

Here we see that Newton’s method will jump between 0 and 1 indefinitely, so both 0 and 1 are not good initial points.

sol=x/.NSolve[f]

success=0;
Do[p=newton1[f,RandomInteger[{0,4}],x];
If[Length[Select[sol,Abs[#-p]<.001&]]>0,success++],{100}]
success

success=0;
Do[p=newton1[f,RandomInteger[{-9,9}],x];
If[Length[Select[sol,Abs[#-p]<.001&]]>0,success++],{100}]
success

success=0;
Do[p=newton1[f,RandomComplex[{-9-9I,9+9I}],x];
If[Length[Select[sol,Abs[#-p]<.001&]]>0,success++],{100}]
success

For this exceptionally tricky polynomial, newton1 works, but in this experiment, it works only 66% of the time for a small integer range of initial points, 91% of the time for a larger range of integer initial points, and 94% of the time for a random complex number. A user-chosen integer point will fail even more often.

g=randomRealPolynomial[4,x,y]/.{y->0}

solc=x/.NSolve[g]

success=0;
Do[p=newton1[g,RandomComplex[{-9-9I,9+9I}],x];
If[Length[Select[solc,Abs[#-p]<.001&]]>0,success++],{1000}]
success

Another experiment the reader may try in Chapter 3, or better in Chapter 6 where we find complex and infinite singularities, is to compare randomintegerpolynomial[4,4,x,y] and randomrealpolynomial[4,x,y] to see which produces more singular polynomials. When I do this, I find that about 20% of the time randomintegerpolynomial[4,4,x,y] produces a singular (including non-square free) polynomial while randomrealpolynomial[4,x,y] has never given me a singular or non-square free polynomial. (See the next section for the code for random polynomials.) In fact, randomintegerpolynomial[d,d+2,x,y] will be singular most of the time for d≥5 although the singularities will often be infinite.

In Chapter 1, we used random real numbers in the system for finding the equation of a line so that our system would have a high probability of being of full rank.

In fact, the routine for finding singular points, including infinite and complex, in Chapter 6 uses randomization in an essential way. In Chapter 9, we will also consider an algorithm that starts by making a random transformation on the curve, and the answer depends somewhat on this random choice. This brings up a new question. Is it better to have an algorithm that returns a different random answer each time or a less good algorithm that always returns the same answer to a problem? One can certainly argue the latter, especially when trying to debug other algorithms using this result. One solution to this question is to choose the random item once and hope that no application just happens to conflict with this choice. In the case of the built-in function NSolve, the Wolfram Language makes this second choice when it encounters a system with infinite solutions and wants to return just one. I prefer the compromise of making this choice once per Wolfram Language session.

2.4 Constructing and Working with Bivariate Polynomials

Since an algebraic curve is the “solution set,” as defined previously, of a real bivariate polynomial, we discuss in this section some aspects of real bivariate polynomials.

In many Wolfram Language procedures handling multivariate functions, one of the inputs is a list of function variables. However, these procedures often expect each variable to be an atomic quantity, that is, a variable could be x, t, a, or even α or Ω but not X[[2]] where X is a list. So in coding user-defined procedures, as we do often in this book, the variables are treated as independent characters.

Considering the previous code to find real solutions of a pair of bivariate polynomials, we would write a shortcut function where F generally refers to the list {f,g}:

rsolve[F_,x_,y_]:=NSolve[F,{x,y},Reals];

rsolve[F_,X_]:= NSolve[F,X,Reals];

The best method in the Wolfram Language to extract coefficients is to use CoefficientRules. Since the coefficients are expected to be Wolfram Language machine numbers in this book, specifying variables is unnecessary unless some of your variables are user-undefined symbolic numbers.

Variables[3.2x y^2+ 5Pi x^2 y-Sin[.2334] x + c y]

CoefficientRules[0.2500319097983996` -0.8807521249099297` x+0.5727398772490768` (x^2)+0.6298021363965578` (x^3)-0.06969551948469332` (x^4)+0.12009965854459015` y+0.3557137203380485` x y+0.003982230515285234` (x^2) y-0.03985635320567571` (x^3) y-0.41567653449909425` (y^2)-0.5547000285700272` x (y^2)+0.6951203004292736` (x^2) (y^2)+0.909106063953359` (y^3)-0.7657856873127402` x (y^3)-0.9546094753841348` (y^4)]

tDeg[f_]:=Max[Table[Total[t],{t,Keys[Association[CoefficientRules[f]]]}]]

Do not ever use Round or Chop on coefficients in a numerical polynomial as that could affect the total degree, which is important information encoded in a polynomial!

N[.2^15]

It will be useful to construct random curves of degree d. To this end, we need a list of ordered pairs with non-negative components totaling d or less.

pExps[d_]:=Module[{pex},
pex[0]={{0,0}};
pex[n_]:=pex[n]=Join[pex[n-1],Table[{n-j,j},{j,0,n}]];
pex[d]]

pExps[3]

randomRealPolynomial[d_,x_,y_]:=FromCoefficientRules[Table[e->RandomReal[{-1,1}],{e,pExps[d]}],{x,y}]

g=randomRealPolynomial[4,x,y]

Since zero is not a numerical concept, a numerical polynomial should have a coefficient for each possible monomial. On the other hand, we may want to generate random exact polynomials.

randomIntegerPolynomial[d_,kk_,x_,y_]:=Module[{k,S,b,j,m,t},
b=Binomial[d+2,2];
If[kk<1 || kk>b, k=b,k=kk];
S=RandomSample[pExps[d],k];
m=Max[Table[Total[t],{t,S}]];
If[m<d,j=RandomInteger[{0,d}]; S=Append[Take[S,k-1],{j,d-j}]];FromCoefficientRules[Table[e->RandomInteger[{-9,9}],{e,S}],{x,y}]]

h=randomIntegerPolynomial[5,6,x,y]

At the start of Chapter 1, we mentioned that our, perhaps strange, algorithm for finding the equation of a line could be generalized for finding the equation of a curve through a number of points. As in the line algorithm, we must determine all the coefficients. A counting exercise based on the algorithm pExps shows that there are binomial coefficient  (({{d+2}, {2}}))=(((d+2)(d+1))/2) coefficients of a bivariate polynomial of degree d. We will use the last equation for getting a random solution out of the otherwise homogeneous system. We only need (({{d+2}, {2}}))-1 conditions. These could be points, tangent directions at points, specification of specific coefficients, etc. We will look at two examples.

ex=pExps[2]

We need a 6×6 coefficient matrix with a random last row leaving five conditions. We want the axes to be parallel to the coordinate axes, which means the monomial involving x y should be zero. x y corresponds to {1, 1} in ex, which is the fifth coefficient, so one row will have a 1 in column 5 and a zero otherwise. We want this conic to be symmetric about the y axes, so the coefficient of x should be zero. This gives a row with 1 in column 2 and zero elsewhere. Finally, we want the conic to go through points {5,0}, {2,3}, and {2.5,-1}. For the first point, the exponents where y has a positive degree have zero coefficients, and the x ones have the 5 to the power, hence the row r3={1,5,0,25,0,0}.

r4= Table[2^e[[1]]*3^e[[2]], {e,ex}]

r5= Table[2.5^e[[1]]*(-1)^e[[2]],{e,ex}]

M1={{0,0,0,0,1,0}, {0,1,0,0,0,0},{1,5,0,25,0,0},r4,r5,RandomReal[{-1,1},6]};
M1//MatrixForm
B1= {0,0,0,0,0,1};

X=LinearSolve[M1,B1]

XY2=Table[x^e[[1]]*y^e[[2]],{e,ex}]

f1=X.XY2

f1=Chop[Expand[f1/Coefficient[f1,y^2]]]

ContourPlot[f1==0,{x,-6,6},{y,-2,4},Epilog->{PointSize[Large],Red, Point[{{5,0},{2,3},{2.5,-1}}]}]

aCurve[pts_,x_,y_]:=Module[{d,P,M,B,n,c,pow},
pow[a_,n_]:=If[n==0,1,a^n];
d=Switch[Length[pts],2,1,5,2,9,3,14,4,20,5,27,6,_,Return["number of points must be 2,5,9,14,20,27"]];
P=pExps[d];
n=Length[P];
M=Table[If[Length[p]==2,pow[p[[1]],e[[1]]]*pow[p[[2]],e[[2]]],pow[p[[1]],e[[1]]]*pow[p[[2]],e[[2]]]*pow[p[[3]],d-e[[1]]-e[[2]]]],{p,pts},{e,P}];
AppendTo[M,RandomReal[{-1,1},n]];
B=Append[Table[0,{n-1}],1];
c=LinearSolve[M,B];
FromCoefficientRules[Table[P[[i]]->c[[i]],{i,n}],{x,y}]
];

Fbow = {(-1+t^2)/t^3, (-1+t^2)/t^4}

We see that t=0 is not in the domain; moreover, t near zero points are far away from the interesting part of the curve, so we restrict our parameter interval to -4≤t≤-.8 and .8≤t≤4. Since we expect a curve of degree 4, we need (({{6}, {2}}))-1=14 points, and by the remarks in the previous section, we should pick points using random t in the interval above.

Pts=Table[Fbow/.{t->r}, {r,Join[RandomReal[{-4,-.8},7],RandomReal[{.8,4},7]]}]

f2=aCurve[Pts,x,y]

f3=Expand[Chop[f2,1.*^-9]/Coefficient[f2,x^4]]

In case we picked the wrong degree, we would have a curve that went through these 14 points of the bow curve but not the actual bow curve.

Expand[x^4-x^2y+y^3/.Thread[{x,y}->Fbow]]

This says, at least, that the bow curve is contained in the curve defined by (x^4)-(x^2)+(y^3)=0. If we did not have an exact system, we could have randomly picked other random values of s to see if Fbow/.{t→s} still approximately satisfy f 3.

Tech Notes

newton1b[f_,p0_,n_,z_]:=Module[{p},
p=p0;
Do[Print[p]; p=p-(f/D[f,x])/.{z->p},{i,n}];
p]

f=x^2+3x-5;
newton1b[f,2/3,5,x]

f2=0.9999999996287839` (x^4)-0.999999999431914` (x^2) y+1.` (y^3);
Assoc=Association[CoefficientRules[f2,{x,y}]];
FromCoefficientRules[Table[k->Round[Assoc[k]],{k,Keys[Assoc]}],{x,y}]