Finding all roots of a polynomial:

A combination of exercises from Section 2.6, Tea Time, p.94:

Andy Long, Spring, 2024

Use Horner’s method twice to produce a polynomial’s function and derivative values at a given value t (exercise #2, p. 93):

t -- the value at which one wishes to evaluate the polynomial
coefs -- the coefficients from highest power to constant term

The basic routine:

basicHorner[t_, coefs_]:=
Module[{n=Length[coefs],result=coefs[[1]],newcoefs={}},
For[i=2,i<=n,i++,
newcoefs=Append[newcoefs,result];
result=coefs[[i]]+result*t;
];
{result,newcoefs}
]
basicHorner[-3,{-4,2,3,-6}]

Use the basic routine twice, to get (and return) function value and derivative:

horner[t_,coefs_]:=
Module[{p=basicHorner[t,coefs],ft,qx}, 
ft=First[p]; (* p(t)*)
qx=First[Rest[p]];(* q(x)'s coefficients: *)
{ft,First[basicHorner[t,qx]]} (* This second call to basicHorner produces q'(t), which is fprime of t *)
]
horner[-3,{-4,2,3,-6}]

Using Newton with Horner on Polynomials:
	newtonhorner (exercise #4, p. 93):

(* help newtonit just does one step in Newton's method: *)
newtonit[x_,fx_,fpx_]:=x-fx/fpx
(* Now on to the function: *)
newtonhorner[a_,coefs_,tol_,maxits_,eps_,verbose_]:=
Module[{i=0,x0=a,x1,y1,pair=horner[a,coefs], diff},
x1=N[newtonit[x0,pair[[1]],pair[[2]]]];
diff=Abs[x1-x0];
While[diff>tol&&i<maxits && Abs[pair[[1]]]> eps,
x0=x1;
pair=horner[x0,coefs];
x1=newtonit[x0,pair[[1]],pair[[2]]];
If[verbose,Print[{x1,pair[[1]]}]];
diff=Abs[x1-x0];
i=i+1;
];
{x1,pair[[1]]}
]
newtonhorner[2,{1,0,-1},10^-6,100,10^-6,True]
newtonhorner[2,{1,0,-1},10^-6,100,10^-6,False]

deflate (exercise #4, p. 93):

deflate[c_,r_,verbose_]:=
Module[{n=Length[c],result,next=First[c]},
result={next};
For[i=2,i<=n,i++,
If[verbose,Print[c[[i]]," ",r*next]];
next=c[[i]]+r*next;
AppendTo[result,next]
];
If[verbose,Print["Error: ",Last[result]]];
Most[result]
]
deflate[{1,0,-1},1,False]
deflate[{-4,2,3,-6},-1.18985,True]

allPolynomialRoots:
	uses newtonhorner, deflate, and quadraticRoots

You need to start with a complex argument, to make sure that you find a root -- otherwise, the routine will surely produce an incorrect result, because it will give an (incorrect) “real root”.

We need a helper: Exercise #1, p. 93:

quadraticRoots[a_,b_,c_]:=
Module[{tmp=b+If[b==0,1,Sign[b]]Sqrt[b^2-4a*c]},
-{2*c/tmp,tmp/(2*a)}
]
quadraticRoots[1,0,-1]
quadraticRoots[1,0,1]

Now let’s put it all together into a recursive algorithm, allPolynomialRoots:

allPolynomialRoots[coefs_,x0_,tol_,maxits_,eps_,verbose_]:=
Module[{n=Length[coefs],root,roots},
If[n>3,
root=newtonhorner[x0,coefs,tol,maxits,eps,verbose][[1]];
roots=Append[allPolynomialRoots[deflate[coefs,root,verbose],root,tol,maxits,eps,verbose],root];
,
roots=quadraticRoots[coefs[[1]],coefs[[2]],coefs[[3]]];
];
Reverse[roots]
]
coefs=1.0{1,0,-1,-3,5};
z0=1+I; (* Complex values are frequently represented by "z" *)
allPolynomialRoots[coefs,z0,10^-6,100,10^-6,False]
(* Let's compare with what Mathematica finds: *)
Roots[
coefs[[1]]*x^4+coefs[[2]]*x^3+coefs[[3]]*x^2+coefs[[4]]*x^1+coefs[[5]]*x^0==0,
x]

#21, p. 94:

coefs=1.0{1,11,-34,-130,-275,819};
z0=1+I; (* There must be a real root since it's a quintic; but if I start with a real, I'll possibly get in trouble on the deflated quartic *)
allPolynomialRoots[coefs,z0,10^-6,100,10^-6,False]
coefs=1.0{5,3,-46,-102,365,1287};
z0=1+I; (* Complex values are frequently represented by "z" *)
allPolynomialRoots[coefs,z0,10^-6,100,10^-6,False]

#18, p. 94:
	Fibonacci finds a root (is he Babylonian?)

coefs=1.0{1,2,10,-20};
allPolynomialRoots[coefs,1,10^-6,100,10^-6,False]
fibonacciSolution=1+
22/60 +
7/60^2 +
42/60^3 +
33/60^4 +
4/60^5 +
40/60^6;
N[fibonacciSolution] (* pert' darn good! *)

#17, p. 94:

coefs=1.0{1,-1.255,-.9838,1.2712};
z0=1.13; (* I can use a real starting value, because the quadratic formula will find the other two roots *)
allPolynomialRoots[coefs,z0,10^-6,100,10^-6,True]

It finds another root, which is actually a little closer: 1.135; then the quadratic formula takes over and finds the desired root.

Plot[x^3-1.255x^2-.9838x+1.2712,{x,1.1,1.14}]