1. [Kincaid and Cheney Problem 3.2#4] Steffensen’s method is given by the iteration

x_{n+1} = x_{n} − \frac{f(x_{n})}{g(x_{n})}  where  g(x) = \frac{f(x + f(x)) − f(x)}{f(x)}.

Show this method is quadratically convergent under suitable hypothesis.

Let r be a solution such that f(r) = 0 and f′(r) ≠ 0. Moreover assume |f′′| is bounded in a neighborhood of r. We shall show that Steffensen’s method is quadratically convergent provided that x_{0} is sufficiently close to r.

Consider the function

κ(x, ξ, η) = \frac{f′′(ξ)(f′(x) − \frac{1}{2}f′′(η)(x − r)) + f′′(η)}{2f′(x) + f′′(ξ)f(x)}.

By assumption the limit supremum of |κ(x, ξ, η)| exists as x → r, ξ → r and η → r. Let C denote this limit and choose δ_{1} > 0 so small that |κ(x, ξ, η)| ≤ 2C for every x, ξ and η such that |x−r| < δ_{1}, |ξ −r| < δ_{1} and |η −r| < δ_{1}. Choose δ_{2} > 0 so small that |x−r| < δ_{2} implies 2|f(x)| ≤ δ_{1}. Now let

δ = min \{\frac{δ_{1}}{2} , δ_{2}, \frac{1}{2C} \}.

Claim that if x_{0} ∈ (r − δ, r + δ) then x_{n} defined according to Steffensen’s method satisfies x_{n} ∈ (r − δ, r + δ) for all n ∈ N.

Let x_{0} ∈ (r − δ, r + δ) and define x_{n} according to Steffensen’s method. By Taylor’s theorem there exists ξ_{n} between x_{n} and x_{n} + f(x_{n}) such that

f(x_{n} + f(x_{n})) = f(x_{n}) + f′(x_{n})f(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n})^{2}.

Therefore

g(x_{n}) = f′(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n}).

It follows that

e_{n+1} = x_{n+1} − r = e_{n} − \frac{f(x_{n})}{f′(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n})}

where e_{n} = x_{n} − r.

Using Taylor’s series again yields

0 = f(r) = f(x_{n} − e_{n}) = f(x_{n}) − f′(x_{n})e_{n} + \frac{1}{2}f′′(η_{n})e_{n}^{2}

for some choice of η_{n} between r and x_{n}. Therefore

e_{n+1} = e_{n} − e_{n} \frac{f′(x_{n}) + \frac{1}{2}f′′(η_{n})e_{n}}{f′(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n})}

= e_{n} − e_{n} \frac{f′(x_{n}) − \frac{1}{2}f′′(ξ_{n})f(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n}) + \frac{1}{2}f′′(η_{n})e_{n}}{f′(x_{n}) + \frac{1}{2}f′′(ξ_{n})f(x_{n})}

= −e_{n} \frac{f′′(ξ_{n})f(x_{n}) + f′′(η_{n})e_{n}}{2f′(x_{n}) + f′′(ξ_{n})f(x_{n})}

= −e_{n} \frac{f′′(ξ_{n})(f′(x_{n})e_{n} − \frac{1}{2}f′′(η_{n})e_{n}^{2}) + f′′(η_{n})e_{n}}{2f′(x_{n}) + f′′(ξ_{n})f(x_{n})}

= −e_{n}^{2} \frac{f′′(ξ_{n})(f′(x_{n}) − \frac{1}{2}f′′(η_{n})e_{n}) + f′′(η_{n})}{2f′(x_{n}) + f′′(ξ_{n})f(x_{n})} = −e_{n}^{2}κ(x_{n}, ξ_{n}, η_{n}).

Claim that x_{n} ∈ (r−δ, r+δ) for every n ∈ N. For induction suppose x_{n} ∈ (r−δ, r+δ). Since ξ_{n} is between x_{n} and x_{n} + f(x_{n}) then

|ξ_{n} − r| ≤ |ξ_{n} − x_{n}| + |x_{n} − r| ≤ |f(x_{n})| + |x_{n} − r| < δ_{2} + δ ≤ δ_{1}.

Since η_{n} is between x and x_{n} then

|η_{n} − r| ≤ δ < δ_{1}.

Thus x_{n}, ξ_{n} and η_{n} satisfy |x_{n} − r| < δ_{1}, |ξ_{n} − r| < δ_{1} and |η_{n} − r| < δ_{1} so consequently |κ(x_{n}, ξ_{n}, η_{n})| < 2C. Therefore

|x_{n+1} − r| = |e_{n+1}| ≤ 2Ce_{n}^{2} = 2C|x_{n} − r|^{2} ≤ 2Cδ^{2} < δ.

This shows x_{n+1} ∈ (r − δ, r + δ) and completes the induction. Moreover, we have shown that x_{n} is quadratically convergent.

2. [Kincaid and Cheney Problem 3.2#6] To compute reciprocals without division, we can solve x = 1/R by finding a zero of the function f(x) = x^{−1} − R. Write a short algorithm to find 1/R by Newton’s method applied to f. Do not use division or exponentiation in your algorithm. For positive R, what starting points are suitable?

Differentiating yields

f′(x) = −x^{−2}.

Therefore

ϕ(x) = x − \frac{f(x)}{f′(x)} = x + \frac{x^{−1} − R}{x^{−2}} = x + x − Rx^{2} = x(2 − Rx)

leads to the iteration formula x_{n+1} = x_{n}(2 − Rx_{n}). To determine what starting points x_{0} are suitable we look for an interval centered at 1/R over which |ϕ′| is strictly less than 1.

|ϕ′(x)| = |2 − 2Rx| < 1  implies  |\frac{1}{R} − x| < \frac{1}{2R}.

Therefore any starting point

x_{0} ∈ (\frac{1}{2R}, \frac{3}{2R})

is suitable. In fact, this interval can be extended to

x_{0} ∈ (0, \frac{3}{2R})

because if x ∈ (0, 1/2R] then

ϕ(x) − x = x(2 − Rx) − x = x(1 − Rx)

implies

\frac{x}{2} ≤ ϕ(x) − x < \frac{1}{2R}.

Consequently, if x_{0} ∈ (0, 1/2R], then x_{n} form an increasing sequence of points that are eventually contained in the interval where |ϕ′| < 1. Similarly

x_{0} ∈ [\frac{3}{2R}, \frac{2}{R})  implies  x_{1} ∈ (0, \frac{3}{4R}].

It follows that the method iteration will converge for any

x_{0} ∈ (0, \frac{2}{R}).

We now procede to find a choice for x_{0} that does not involve performing any divisions. Recall that R is stored on the computer using a binary expansion of the form

R = ∑_{i=1}^{m} b_{i}2^{n−i} = σ2^{n}  where  σ ∈ [1/2, 1).

Therefore taking

x_{0} = σ2^{−n} < \frac{1}{σ2^{n}} = \frac{1}{R}

provides a suitable initial condition by simply changing the sign of the exponent in the floating point representation while doing no division or exponentiation to the mantissa. The following C program implements this approach.

 1 /*
 2  newtdiv.c -- Find 1/R without hardware divide instruction
 3  Written October 17, 2012 by Eric Olson for Math 701
 4 */
 5 #include <stdio.h>
 6 #include <math.h>

7

 8 typedef union {
 9  struct {
10  long long unsigned sigma:52;
11  unsigned n:11;
12  unsigned s:1;
13  } bits;
14  double value;
15 } ieeebits;

16

17 double invert(double R){
18  double x,x0;
19  ieeebits guess;
20  guess.value=R;
21  guess.bits.n=2044-guess.bits.n;
22  x0=guess.value;
23  int i;
24  for(i=1;i<100;i++){
25  x=x0*(2-R*x0);
26  if(x==x0) return x;
27  x0=x;

28 }

29  return NAN;

30 }

31

32 int main(){
33  printf(
34  "newtdiv.c -- Find 1/R without hardware divide instruction\n"
35  "Written October 17, 2012 by Eric Olson for Math 701\n");

36

37  double R[]={ 1.0 , 2.0, 0.999, 214.32, 0.000128 };
38  int N=sizeof(R)/sizeof(double);
39  printf("\n%2s %22s %22s %22s\n","n","R","1/R","|1-R(1/R)|");
40  int n;
41  for(n=0;n<N;n++){
42  double x=R[n],y=invert(x);
43  printf("%2d %22.15e %22.15e %22.15e\n",n+1,x,y,fabs(1.0-x*y));

44 }

45  return 0;

46 }

The output is

newtdiv.c -- Find 1/R without hardware divide instruction
Written October 17, 2012 by Eric Olson for Math 701
n                      R                    1/R             |1-R(1/R)|
1  1.000000000000000e+00  1.000000000000000e+00  0.000000000000000e+00
2  2.000000000000000e+00  5.000000000000000e-01  0.000000000000000e+00
3  9.990000000000000e-01  1.001001001001001e+00  8.749511531958021e-17
4  2.143200000000000e+02  4.665920119447556e-03  9.063930161978817e-17
5  1.280000000000000e-04  7.812500000000000e+03  4.526544070126981e-17

3. [Kincaid and Cheney Problem 3.2#17] Which of the following sequences converge

quadratically.

a. 1/n^{2}, b. 1/2^{2^{n}}, c. 1/\sqrt{n}, d. 1/e^{n}, e. 1/n^{n}

A sequence x_{n} converging to a limit L, converges quadratically if there is a constant C and an integer N such that |x_{n+1} − L| ≤ C|x_{n} − L|^{2} for every n ≥ N. Since L = 0 for each of the sequences in question and each sequence is moreover positive, then this inequality may be rewritten as

\frac{x_{n+1}}{x_{n}^{2}} ≤ C  for every  n ≥ N.

This inequality will hold for some C and N large enough if and only if

lim sup_{n→∞} \frac{x_{n+1}}{x_{n}^{2}}

is bounded. We proceed as follows:

lim sup_{n→∞} \frac{1/(n + 1)^{2}}{(1/n^{2})^{2}} = lim_{n→∞} \frac{n^{4}}{n^{2} + 2n + 1} = ∞

lim sup_{n→∞} \frac{1/2^{2^{n+1}}}{(1/2^{2^{n}})^{2}} = lim_{n→∞} \frac{4^{2^{n}}}{4^{2^{n}}} = 1

lim sup_{n→∞} \frac{1/\sqrt{n + 1}}{(1/\sqrt{n})^{2}} = lim_{n→∞} \frac{n}{\sqrt{n + 1}} = ∞

lim sup_{n→∞} \frac{1/e^{n+1}}{(1/e^{n})^{2}} = lim_{n→∞} \frac{e^{2n}}{e^{n+1}} == lim_{n→∞} e^{n−1} = ∞

lim sup_{n→∞} \frac{1/(n + 1)^{n+1}}{(1/n^{n})^{2}} = lim_{n→∞} \frac{n^{2n}}{(n + 1)^{n+1}} ≥ lim_{n→∞} \frac{n^{2n}}{(2n)^{n+1}} = lim_{n→∞} \frac{1}{4}(\frac{n}{2})^{n} = ∞.

Therefore only 1/2^{2^{n}} converges quadratically while the others don’t.

4. [Kincaid and Cheney Computer Problem 3.2#1] Write a computer program to solve the equation x = tan x by Newton’s method. Find the roots nearest 4.5 and 7.7.

The program is

 1 /*
 2  newtontan.c -- Newton's method to solve tan(x) = x.
 3  Kincaid and Cheney Computer Exercise 3.2#1
 4  Written October 15, 2012 for Math 701 by Eric Olson
 5 */

6

 7 #include <stdio.h>
 8 #include <stdlib.h>
 9 #include <math.h>

10

11 #define N 100

12

13 double phi(double x){
14  double t=tan(x);
15  double f=t-x;
16  double df=t*t;
17  return x-f/df;

18 }

19 double newton(double x0){
20  double x[N];
21  x[0]=x0;
22  int n;
23  for(n=1;n<N;n++) {
24  x[n]=phi(x[n-1]);
25  if(x[n]==x[n/2]) return x[n];

26 }

27  printf("Newton failed to converge!\n");
28  exit(1);

29 }

30 double x0[]={ 4.5, 7.7 };
31 int K=sizeof(x0)/sizeof(double);
32 main(){
33 printf(
34  "Kincaid and Cheney Computer Exercise 3.2#1\n"
35  "Written October 15, 2012 for Math 701 by Eric Olson\n");

36

37  printf("\n%4s %24s %24s\n","n","x0","root");
38  int n;
39  for(n=0;n<K;n++)
40  printf("%4d %24.15e %24.15e\n",
41  n+1,x0[n],newton(x0[n]));
42  return 0;

43 }

and the output is

   Kincaid and Cheney Computer Exercise 3.2#1
   Written October 15, 2012 for Math 701 by Eric Olson
      n                       x0                     root
      1    4.500000000000000e+00    4.493409457909064e+00
      2    7.700000000000000e+00    7.725251836937707e+00

5. [Kincaid and Cheney Computer Problem 3.2#8] Program the Newton algorithm in complex arithmetic and test it on these functions with the given starting points. a. f_{1}(z) = z^{2} + 1, z_{0} = 3 + i b. f_{2}(z) = z + sin z − 3, z_{0} = 2 − i c. f_{3}(z) = z^{4} + z^{2} + 2 + 3i, z_{0} = 1

We shall use the Maple code generation library to create the functions ϕ(z) = z−f(z)/f′(z) needed for Newton’s iteration. The Maple script

1 restart;
2 #   makephi.mpl -- Make phi to solve F=x for Newton's Iteration
3 #   Written October 15, 2012 for Math 701 by Eric Olson
4 #
5 #   F is an expression in x
6 #   N is the number of the function
7 #
8 kernelopts(printbytes=false);
9 fname:=sprintf("phi%d.i",N);
10 fd:=open(fname,WRITE);
11 fprintf(fd,"/* %s -- Created by Maple */\n",fname);
12 close(fd);

13

14 with(codegen,C,optimize):
15 cmp:=[I=cI,cos=ccos,sin=csin,exp=cexp];
16 tmp:=[seq(t||i=T[i],i=1..64)];

17

18 f:=F;
19 df:=diff(f,z);
20 r1:=[phi=z-f/df];
21 r2:=[optimize(r1,'tryhard')];
22 r3:=subs(cmp,r2);
23 r4:=subs(tmp,r3);
24 C(r4,filename=fname);

takes an expression F and an integer N and creates the file phiN.i containing C code to compute ϕ in terms of z. The script can be used to generate code for the functions f_{n} given in the problem with the commands

      $ maple -q -DN=1 -DF="z^2+1" makephi.mpl
      $ maple -q -DN=2 -DF="z+sin(z)-3" makephi.mpl
      $ maple -q -DN=3 -DF="z^4+z^2+2+3*I" makephi.mpl

The resulting files are

1 /* phi1.i -- Created by Maple */
2  phi = z-(z*z+1.0)/z/2.0;
1 /* phi2.i -- Created by Maple */
2  phi = z-(z+csin(z)-3.0)/(1.0+ccos(z));
1 /* phi3.i -- Created by Maple */
2  T[0] = z*z;
3  phi = z-(T[0]*T[0]+T[0]+2.0+3.0*cI)/(4.0*T[0]+2.0)/z;

Notice that phi3.i contains an optimization where z^{2} is stored in a temporary variable. These files are included in our main program newtcomp.c which is given by

 1 /*
 2  newtcomp.c -- Newton's method to solve f(z) = 0.
 3  Kincaid and Cheney Computer Exercise 3.2#8
 4  Written October 15, 2012 for Math 701 by Eric Olson
 5 */

6

 7 #include <stdio.h>
 8 #include <stdlib.h>
 9 #include <math.h>
10 #include <complex.h>

11

12 #define cI I
13 typedef double complex Complex;
14 typedef Complex (*Function)(Complex z);

15

16 Complex phi1(Complex z){
17  Complex T[64],phi;
18 #include "phi1.i"
19  return phi;

20 }

21 Complex phi2(Complex z){
22  Complex T[64],phi;
23 #include "phi2.i"
24  return phi;

25 }

26 Complex phi3(Complex z){
27  Complex T[64],phi;
28 #include "phi3.i"
29  return phi;

30 }

31

32 #define N 100
33 Complex iterate(Function phi,Complex z0){
34  Complex z[N];
35  z[0]=z0;
36  int n;
37  for(n=1;n<N;n++) {
38  z[n]=phi(z[n-1]);
39  if(z[n]==z[n/2]) return z[n];

40 }

41  printf("Newton failed to converge!\n");
42  exit(0);

43 }

44 Complex z0[]={ 3+I, 2-I, 1 };
45 Function phi[]={ phi1, phi2, phi3 };
46 int K=sizeof(z0)/sizeof(Complex);

47

48 main(){
49  printf(
50  "Kincaid and Cheney Computer Exercise 3.2#8\n"
51  "Written October 15, 2012 for Math 701 by Eric Olson\n");

52

53  printf("\n%4s %16s %16s %16s %16s\n",
54  "n","re(z0)","im(z0)","re(z)","im(z)");
55  int n;
56  for(n=0;n<K;n++)
57  printf("%4d %16.8e %16.8e %16.8e %16.8e\n",
58  n+1,z0[n],iterate(phi[n],z0[n]));
59  return 0;

60 }

The output is

   Kincaid and Cheney Computer Exercise 3.2#8
   Written October 15, 2012 for Math 701 by Eric Olson
      n           re(z0)           im(z0)            re(z)            im(z)
      1   3.00000000e+00   1.00000000e+00   0.00000000e+00   1.00000000e+00
      2   2.00000000e+00  -1.00000000e+00   2.17975707e+00   0.00000000e+00
      3   1.00000000e+00   0.00000000e+00   1.02469479e+00  -7.88542500e-01

6. [Kincaid and Cheney Computer Problem 3.2#9] Use Steffensen’s method to find the first 10 positive roots of the equation x = tan x.

First graph the equation y = tan x and y = x to identify on what intervals the roots will be. The gnuplot script

    1 set terminal pstex
    2 set output "xtanx.tex"
    3 set key outside
    4 set key width 1
    5 set samples 1000
    6 set size 0.8,0.5
    7 plot [0:30] [0:30] tan(x) ti "$y=\\tan x$",x ti "$y=x$"

produces the output

y = x y = tan x 30 25 20 15 10 5 0 30 25 20 15 10 5 0

There is root between each vertical asymptote. The first 10 positive roots satisfy

x_{n} ∈ (\frac{(2n − 1)π}{2} , \frac{(2n + 1)π}{2} )  where  n = 1, 2, …, 10.

Steffensen’s method requires the initial guess be very close to the correct root to converge because of the aymptotes. To overcome this problem we check whether each Steffensen’s iterate remains in the interval where the n-th root is located. If it does not remain in the interval, then we take a bisection step and try again. The program is

    1 /*
    2  stefbisect.c -- Steffensen's method to solve tan(x) = x.
    3  Kincaid and Cheney Computer Exercise 3.2#9
    4  Written October 15, 2012 for Math 701 by Eric Olson
    5 */

6

    7 #include <stdio.h>
    8 #include <stdlib.h>
    9 #include <math.h>

10

   11 #define N 3000
   12 typedef long double Real;

13

   14 Real f(Real x){
   15  return tan(x)-x;

16 }

17 Real df(Real x){
18  Real t=tan(x);
19  return t*t;

20 }

21 Real g(Real x){
22  Real t=f(x);
23  return (f(x+t)-t)/t;

24 }

25 Real steffensen(Real a, Real b, Real x0){
26  Real x[N];
27  x[0]=x0;
28  int n;
29  for(n=1;n<N;n++){
30  Real r=x[n-1];
31  Real t1=f(r);
32  if(t1==0.0) return r;
33  Real t2=g(r);
34  if(t2==0.0) return r;
35  r=r-t1/t2;
36  if(r<a||r>b){
37  r=(a+b)/2;
38  if(f(r)<0) a=r;
39  else b=r;

40 }

41  x[n]=r;
42  if(x[n]==x[n/2]) return x[n];

43 }

44  printf("Steffensen failed to converge!\n");
45  exit(0);

46 }

47 int main () {
48  printf(
49  "Kincaid and Cheney Computer Exercise 3.2#9\n"
50  "Written October 15, 2012 for Math 701 by Eric Olson\n");
51  int n;
52  printf("\n%4s %30s\n","n","x[n]");
53  for(n=1;n<=10;n++){
54  Real r=steffensen((2*n-1)*M_PI/2,(2*n+1)*M_PI/2,3*n*M_PI/2);
55  printf("%4d %30.20Le\n",n,r);

56 }

57  return 0;

58 }

and the output is

   Kincaid and Cheney Computer Exercise 3.2#9
   Written October 15, 2012 for Math 701 by Eric Olson
      n                           x[n]
      1     4.49340945790906443839e+00
      2     7.72525183693770697282e+00
      3     1.09041216594289006193e+01
      4     1.40661939128314738452e+01
      5     1.72207552719307686595e+01
      6     2.03713029592875630482e+01

      7     2.35194524986890063350e+01
      8     2.66660542588126723597e+01
      9     2.98115987908929581697e+01
     10     3.29563890398224784671e+01

An alternative method is to invert the equations and work with arctan(y). This is implemented in the program

 1 /*
 2  stefatan.c -- Steffensen's method to solve tan(x) = x.
 3  Kincaid and Cheney Computer Exercise 3.2#9
 4  An alternative approach with arctan.
 5  Written October 15, 2012 for Math 701 by Eric Olson
 6 */

7

 8 #include <stdio.h>
 9 #include <math.h>

10

11 int k;
12 double f(double x){
13  return atan(x)-x+k*M_PI;

14 }

15 int main(){
16  printf(
17  "Kincaid and Cheney Computer Exercise 3.2#9\n"
18  "An alternative approach with arctan.\n"
19  "Written October 15, 2012 for Math 701 by Eric Olson\n");

20

21  printf("\n%4s %24s\n","n","x");
22  for(k=0;k<=10;k++){
23  double x=1.0;
24  int n;
25  for(n=1;n<=10;n++){
26  double y=f(x);
27  double d=f(x+y)-y;
28  if(d==0.0 || y==0.0) break;
29  x=x-y*y/d;

30 }

31  printf("%4d %24.15e\n",k,x);

32 }

33  return 0;

34 }

which produces the output

   Kincaid and Cheney Computer Exercise 3.2#9
   An alternative approach with arctan.
   Written October 15, 2012 for Math 701 by Eric Olson
      n                        x
      0    1.156162205308371e-02
      1    4.493409457909064e+00
      2    7.725251836937707e+00
      3    1.090412165942890e+01
      4    1.406619391283147e+01
      5    1.722075527193077e+01

        6    2.037130295928756e+01
        7    2.351945249868901e+01
        8    2.666605425881267e+01
        9    2.981159879089296e+01
       10    3.295638903982248e+01

It is worth noting that this second method involves transforming the original problem into something easier to solve and is significantly more efficient than the first.

7. [Kincaid and Cheney Problem 3.5#1] Use Horner’s algorithm to find p(4) where p(z) = 3z^{5} − 7z^{4} − 5z^{3} + z^{2} − 8z + 2.

We arrange the calculation as

| 3 −7 −5 1 −8  2
4 |  12 20 60 244 944
| 3 5 15 61 236 946

Therefore p(4) = 946.

8. [Kincaid and Cheney Problem 3.5#5] For the polynomial

p(z) = 3z^{5} − 7z^{4} − 5z^{3} + z^{2} − 8z + 2

find a disk centered at the origin that contains all the roots. By the localization theorem all roots lie inside an open disk centered at the origin of radius

\begin{aligned}ρ & = 1 + |a_{n}|^{−1} max \{ |a_{k}| : k = 0, 1, …, n − 1 \} \\ & = 1 + 3^{−1} max\{2, 8, 1, 5, 7\} = 1 + 8/3 = 11/3. \end{aligned}

9. [Kincaid and Cheney Problem 3.5#19] Prove that if p is a polynomial of degree n having real coefficients, then

(n − 1)(p′(x))^{2} ≥ np(x)p′′(x).

Assume that the roots are real.

We prove this statement by induction. Let P_{n} be the set of all polynomials of degree n having real coefficient and real roots.

Base case: If n = 1 then the statement reduces to the fact that 0 ≥ 0.

Induction step: Suppose (n − 1)(p′_{n}(x))^{2} ≥ np_{n}(x)p′′_{n}(x) for every p_{n} ∈ P_{n}. Let p_{n+1} ∈ P_{n+1}. Claim that n(p′_{n+1}(x))^{2} ≥ (n + 1)p_{n+1}(x)p′′_{n+1}(x).

Let r ∈ R be a root of p_{n+1}. Then there is p_{n} ∈ P_{n} such that

p_{n+1}(x) = (x − r)p_{n}(x)

and therefore

p′_{n+1}(x) = p_{n}(x) + (x − r)p′_{n}(x)
p′′_{n+1}(x) = 2p′_{n}(x) + (x − r)p′′_{n}(x).

For notational simplicity write p_{n} in place of p_{n}(x) and so forth. Estimate using the induction hypothesis and completing the square to obtain

(n + 1)p_{n+1}p′′_{n+1} = (n + 1)(x − r)p_{n}(2p′_{n} + (x − r)p′′_{n})

= 2(n + 1)(x − r)p_{n}p′_{n} + (n + 1)(x − r)^{2}p_{n}p′′_{n}

≤ 2(n + 1)(x − r)p_{n}p′_{n} + (n + 1)(x − r)^{2} \frac{n − 1}{n}(p′_{n})^{2}

= 2(n + 1)(x − r)p_{n}p′_{n} + \frac{n^{2} − 1}{n}(x − r)^{2}(p′_{n})^{2}

= n{[(x − r)p′_{n}]^{2} + 2p_{n}[(x − r)p′_{n}]}

− \frac{1}{n}\{[(x − r)p′_{n}]^{2} − 2np_{n}[(x − r)p′_{n}]\}

= n((x − r)p′_{n} + p_{n})^{2} − np_{n}^{2} − \frac{1}{n}((x − r)p′_{n} − np_{n})^{2} + np_{n}^{2}

= n((x − r)p′_{n} + p_{n})^{2} − \frac{1}{n}((x − r)p′_{n} − np_{n})^{2}

≤ n((x − r)p′_{n} + p_{n})^{2} = n(p′_{n+1})^{2}.

This completes the induction.

10. [Kincaid and Cheney Computer Problem 3.5#5] Write a computer routine that uses Newton’s method and deflation for finding all roots of a polynomial. Test the routine using Wilkinson’s polynomial

p_{20}(x) = ∏_{n=1}^{20}(x − n).

Using Newton’s method to find a root of an arbitrary polynomial requires care in choosing the initial guess x_{0}. In the code that follows we use Theorem 3 in Kincaid and Cheney page 110 on the localization of polynomial roots to guide choices for the initial guess and to determine whether Newton’s iteration has diverged.

Let ρ_{max} be ϕ times the bound obtained by Theorem 3 where ϕ = (1 + \sqrt{5})/2 and ρ_{n} = μϕ^{n} where μ is the unit roundoff error. Choose N to be the largest integer such that ρ_{N} ≤ ρ_{max}. The following heuristics are used to guess x_{0}. Since it is desired to find roots of smallest magnitude first, we always try x_{0} = 0 as the first guess. If this doesn’t work we guess randomly an x_{0} inside a circle centered ρ_{n} where n = 0, 1, …, N until a root is found.

It is worth noting our strategy for guessing x_{0} is simple, but may not lead to a root in all cases. If a certain maximum number of initial guesses K has been tried without finding a root the program exits with an error.

The program roots.c is given by

    1 /*
    2  roots.c -- Newton's method to find polynomial roots
    3  Kincaid and Cheney Computer Exercise 3.5#5
    4  Written September 25, 2012 for Math 701 by Eric Olson

5

    6  This code assumes all roots are simple roots.

7

    8  Input file format is

9

   10  n c0 c1 c2 c3 c4 ... cn

11

   12  where n is degree of polynomial and ck are real coefficients.
   13  End of data is indicated by n equal to -1.

14

   15 Version 1 to 2 changes:

16

   17  The input file format now allows complex coefficients.

18

   19  If n > 1000 then n - 1000 is the degree of the polynomial and
   20  the format for the coefficients is

21

   22  n c0.x c0.y c1.x c1.y ... cn.x cn.y

23

   24 */

25

   26 #include <stdio.h>
   27 #include <stdlib.h>
   28 #include <complex.h>

29 #include <math.h>
30 #include <fenv.h>
31 #define FE_BAD (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW)

32

33 /*
34  Compile time constants are

35

36  K maximum number of initial guesses
37  N maximum number of iterations per guess
38  L maximum degree of polynomial
39  P number of pollishing iterations
40 */

41

42 #define K 1000
43 #define N 1000
44 #define L 30
45 #define P 7

46

47 typedef long double complex Complex;
48 typedef long double Real;

49

50 typedef struct {
51  int d;
52  Complex c[L+1];

53

54 } polynomial;
55 typedef struct {
56  Complex y,dy;
57 } jet;

58

59 Complex guess(Real rho){
60  Real r=rho*random()/RAND_MAX;
61  Real theta=2*M_PI*random()/RAND_MAX;
62  return r*cexp(I*theta);

63 }

64

65 Real getrho(polynomial *p){
66  Real cmax=0.0;
67  int i;
68  for(i=0;i<p->d;i++){
69  Real r=cabsl(p->c[i]);
70  if(r>cmax) cmax=r;

71 }

72  return cmax/cabsl(p->c[p->d])+1.0;

73 }

74 jet polyjet(polynomial *p,Complex z){
75  int j;
76  jet r;
77  r.y=p->c[p->d];
78  r.dy=0;
79  for(j=1;j<=p->d;j++){
80  r.dy=r.dy*z+r.y;
81  r.y=r.y*z+p->c[p->d-j];

82 }

83  return r;

 85 polynomial deflate(polynomial *p,Complex z){
 86  polynomial q;
 87  int i,j;
 88  q.d=p->d-1;
 89  q.c[q.d]=p->c[p->d];
 90  for(j=1;j<=q.d;j++) q.c[q.d-j]=q.c[p->d-j]*z+p->c[p->d-j];
 91  return q;

92 }

 93 Complex findroot(polynomial *p){
 94  Complex z[N];
 95  Real phi=(1.0+sqrt(5.0))/2.0;
 96  Real rhomin=pow(2.0,-54.0);
 97  Real rhomax=phi*getrho(p);
 98  int i,j,k;
 99  z[0]=guess(0.0);
100  Real rho=rhomin;
101  for(k=0;k<K;k++){
102  for(i=1;i<N;i++){
103  jet r=polyjet(p,z[i-1]);
104  z[i]=z[i-1]-r.y/r.dy;
105  if(z[i]==z[i>>1]) return z[i];
106  if(fetestexcept(FE_BAD)){
107  feclearexcept(FE_BAD);
108  break;

109 }

110  if(cabsl(z[i])>rhomax) break;

111 }

112  z[0]=guess(rho);
113  rho*=phi;
114  if(rho>rhomax) rho=rhomin;

115 }

116  printf("Failed to find a root!\n");
117  exit(4);

118 }

119 Complex pollish(polynomial *p,Complex z0){
120  Complex z=z0;
121  int i;
122  for(i=0;i<P;i++){
123  jet r=polyjet(p,z);
124  z=z-r.y/r.dy;

125 }

126  if(fetestexcept(FE_BAD)){
127  feclearexcept(FE_BAD);
128  puts("Floating point exception while pollishing!");
129  exit(3);

130 }

131  return z;

132 }

133 int cmpmodulus(const void *z1,const void *z2){
134  Real r1=cabsl(*(Complex *)z1);
135  Real r2=cabsl(*(Complex *)z2);
136  if(r1>r2) return 1;
137  if(r1<r2) return -1;
138  return 0;

140

141 int main(){
142  printf(
143  "roots -- Newton's method to find polynomial roots Version 2\n"
144  "Written September 25, 2012 for Math 701 by Eric Olson\n\n");
145  int k;
146  for(k=1;;k++){
147  Complex z[L];
148  polynomial p;
149  int i,cflag;
150  if(scanf("%d",&i)!=1) exit(0);
151  if(i<0) exit(0);
152  if(i>1000) {
153  p.d=i-1000; cflag=1;
154  } else {
155  p.d=i; cflag=0;

156 }

157  if(p.d>L) {
158  puts("Memory exceeded!");
159  exit(1);

160 }

161  for(i=0;i<=p.d;i++){
162  Real x,y;
163  if(scanf("%Lg",&x)!=1){
164  puts("Incomplete input data!");
165  exit(2);

166 }

167  if(cflag) {
168  if(scanf("%Lg",&y)!=1){
169  puts("Incomplete complex number!");
170  exit(5);

171 }

172  } else y=0;
173  p.c[i]=x+I*y;

174 }

175  for(i=0;i<p.d;i++){
176  int j;
177  polynomial q=p;
178  for(j=0;j<i;j++) q=deflate(&q,z[j]);
179  z[i]=pollish(&p,findroot(&q));
180  qsort(z,i,sizeof(Complex),cmpmodulus);

181 }

182  printf("Coefficents of polynomial %d:\n",k);
183  printf("%4s %24s %24s\n","n","Re(Cn)","Im(Cn)");
184  for(i=0;i<=p.d;i++)
185  printf("%4d %24.15Le %24.15Le\n",
186  i,creall(p.c[i]),cimagl(p.c[i]));
187  printf("\nRoots of polynomial %d:\n",k);
188  printf("%4s %24s %24s\n","n","Re(Zn)","Im(Zn)");
189  for(i=0;i<p.d;i++)
190  printf("%4d %24.15Le %24.15Le\n",
191  i+1,creall(z[i]),cimagl(z[i]));
192  printf("\n");

193 }

and Wilkinson’s polynomial wilkpoly.dat is given by

 1 20
 2 2432902008176640000
 3 -8752948036761600000
 4 13803759753640704000
 5 -12870931245150988800
 6 8037811822645051776
 7 -3599979517947607200
 8 1206647803780373360
 9 -311333643161390640
10 63030812099294896
11 -10142299865511450
12 1307535010540395
13 -135585182899530
14 11310276995381
15 -756111184500
16 40171771630
17 -1672280820
18 53327946
19 -1256850
20 20615
21 -210

22 1

23

24 -1

Running the program with the command

   $ ./roots <wilkpoly.dat

produces the output

   roots -- Newton's method to find polynomial roots Version 2
   Written September 25, 2012 for Math 701 by Eric Olson
   Coefficents of polynomial 1:
      n                   Re(Cn)                   Im(Cn)
      0    2.432902008176640e+18    0.000000000000000e+00
      1   -8.752948036761600e+18    0.000000000000000e+00
      2    1.380375975364070e+19    0.000000000000000e+00
      3   -1.287093124515099e+19    0.000000000000000e+00
      4    8.037811822645052e+18    0.000000000000000e+00
      5   -3.599979517947607e+18    0.000000000000000e+00
      6    1.206647803780373e+18    0.000000000000000e+00
      7   -3.113336431613906e+17    0.000000000000000e+00
      8    6.303081209929490e+16    0.000000000000000e+00
      9   -1.014229986551145e+16    0.000000000000000e+00
     10    1.307535010540395e+15    0.000000000000000e+00
     11   -1.355851828995300e+14    0.000000000000000e+00
     12    1.131027699538100e+13    0.000000000000000e+00
     13   -7.561111845000000e+11    0.000000000000000e+00
     14    4.017177163000000e+10    0.000000000000000e+00
     15   -1.672280820000000e+09    0.000000000000000e+00
     16    5.332794600000000e+07    0.000000000000000e+00

17   -1.256850000000000e+06    0.000000000000000e+00
18    2.061500000000000e+04    0.000000000000000e+00
19   -2.100000000000000e+02    0.000000000000000e+00
20    1.000000000000000e+00    0.000000000000000e+00
Roots of polynomial 1:
 n                   Re(Zn)                   Im(Zn)
 1    1.000000000000000e+00    0.000000000000000e+00
 2    1.999999999999999e+00   -0.000000000000000e+00
 3    3.000000000000035e+00    0.000000000000000e+00
 4    3.999999999998008e+00   -0.000000000000000e+00
 5    5.000000000023786e+00    0.000000000000000e+00
 6    6.000000000154543e+00   -0.000000000000000e+00
 7    6.999999994901812e+00    0.000000000000000e+00
 8    7.999999968307939e+00   -0.000000000000000e+00
 9    8.999999979100598e+00    0.000000000000000e+00
10    9.999999912241915e+00   -0.000000000000000e+00
11    1.100000027036567e+01    0.000000000000000e+00
12    1.199999918601061e+01   -0.000000000000000e+00
13    1.300000450574287e+01    0.000000000000000e+00
14    1.400000206872498e+01    0.000000000000000e+00
15    1.500000473937224e+01    0.000000000000000e+00
16    1.600000224002417e+01   -0.000000000000000e+00
17    1.699999906543233e+01    0.000000000000000e+00
18    1.799999979181205e+01   -0.000000000000000e+00
19    1.900000006055330e+01    0.000000000000000e+00
20    1.999999998278379e+01   -0.000000000000000e+00

The roots are sorted and the polynomial is deflated in order of roots in increasing magnitude to reduce rounding error. To further reduce rounding error, each root obtained using Newton’s method on a reduced polynomial has been refined by applying Newton’s method to the original polynomial. However, note that although the long double complex data type used in this program retains over 18 digits of precision the resulting roots are accurate to only 7 significant digits. This is because the conditioning number corresponding to the problem of finding the roots of a polynomial is very large.

11. [Kincaid and Cheney Problem 3.6#1] Solve the system of equations

\{ \begin{matrix}x − 2y + y^{2} + y^{3} − 4 = 0 \\ −x − y + 2y^{2} − 1 = 0 \end{matrix}

by the homotopy method starting with the point (x_{0}, y_{0}) = (0, 0). Perform all calculations without recourse to numerical methods.

Following Kincaid and Cheney Example 2 on page 133 we write

\begin{aligned}h(t, x, y) & = tf(x, y) + (1 − t)(f(x, y) − f(x_{0}, y_{0})) \\ & = f(x, y) + (1 − t)f(0, 0) = [ \begin{matrix}x − 2y + y^{2} + y^{3} − 8 − 4t \\ −x − y + 2y^{2} − 2 − t \end{matrix} ] . \end{aligned}

Therefore

\frac{dh}{ds} = Dh\frac{d}{ds} [ \begin{matrix}t \\ x \\ y \end{matrix} ] = [ \begin{matrix}−4 & 1 & 3y^{2} + 2y − 2 \\ −1 & −1 & 4y − 1 \end{matrix} ] [ \begin{matrix}t′ \\ x′ \\ y′ \end{matrix} ] = 0.

We use Cramer’s rule via Lemma 1 on page 135 to choose a solution (t(s), x(s), y(s)) to this equation that is continuous in s. We obtain the system

t′ = det [ \begin{matrix}1 & 3y^{2} + 2y − 2 \\ −1 & 4y − 1 \end{matrix} ] = 3y^{2} + 6y − 3
x′ = − det [ \begin{matrix}−4 & 3y^{2} + 2y − 2 \\ −1 & 4y − 1 \end{matrix} ] = −3y^{2} + 14y − 2
y′ = det [ \begin{matrix}−4 & 1 \\ −1 & −1 \end{matrix} ] = 5

This system of ordinary differential equations is triangular and solved by integrating first for y, then for x and finally for t. The computation is

\begin{aligned}y(s) & = y_{0} + ∫_{0}^{s} 5 du = 5s \\ x(s) & = x_{0} + ∫_{0}^{s} ( − 3(5u)^{2} + 14(5u) − 2)du = −25s^{3} + 35s^{2} − 2s \\ t(s) & = ∫_{0}^{s} (3(5u)^{2} + 6(5u) − 3)du = 25s^{3} + 15s^{2} − 3s \end{aligned}

Now solve for s such that t(s) = 1. We obtain

25s^{3} + 15s^{2} − 3s − 1 = 0

which we first try to factor over Q. By the rational root theorem the possible roots are

\{ ± 1, ±\frac{1}{5}, ±\frac{1}{25}\}.

Testing these yields that s_{1} = −1/5 is a root so we obtain

(5s + 1)(5s^{2} + 2s − 1) = 0.

The remaining roots can be found by the quadratic formula as

s_{2} = \frac{−2 + \sqrt{4 + 20}}{10} = \frac{−1 + \sqrt{6}}{5}  and  s_{3} = \frac{−1 − \sqrt{6}}{5}.

The corresponding solutions are

\begin{aligned}(x(s_{1}), y(s_{1})) & = (2, −1) \\ (x(s_{2}), y(s_{2})) & = (14 − 5\sqrt{6}, −1 + \sqrt{6} ) \\ (x(s_{3}), y(s_{3})) & = (14 + 5\sqrt{6}, −1 − \sqrt{6} ). \end{aligned}

We now check the solutions using Maple. The Maple script is

1 restart;
2 f1:=(x,y)->x-2*y+y^2+y^3-4;
3 f2:=(x,y)->-x-y+2*y^2-1;

4

5 x1:=2; y2:=-1;
6 'f1(x1,y1)'=simplify(f1(x1,y1));
7 'f2(x1,y1)'=simplify(f2(x1,y1));

8

9 x2:=14-5*sqrt(6); y2:=-1+sqrt(6);
10 'f1(x2,y2)'=simplify(f1(x2,y2));
11 'f2(x2,y2)'=simplify(f2(x2,y2));

12

13 x3:=14+5*sqrt(6); y3:=-1-sqrt(6);
14 'f1(x3,y3)'=simplify(f1(x3,y3));
15 'f2(x3,y3)'=simplify(f2(x3,y3));

and the output is

                                                  2    3
                       f1 := (x, y) -> x - 2 y + y  + y  - 4
                                                     2
                         f2 := (x, y) -> -x - y + 2 y  - 1
                                      x1 := 2
                                      y2 := -1
                                                    2     3
                         f1(x1, y1) = -2 - 2 y1 + y1  + y1
                                                       2
                            f2(x1, y1) = -3 - y1 + 2 y1
                                               1/2
                                 x2 := 14 - 5 6

           1/2
y2 := -1 + 6
f1(x2, y2) = 0
f2(x2, y2) = 0
            1/2
x3 := 14 + 5 6
           1/2
y3 := -1 - 6
f1(x3, y3) = 0
f2(x3, y3) = 0

which shows we have indeed found three roots.

12. [Kincaid and Cheney Problem 3.6#2] Consider the homotopy

h(t, x) = tf(x) + (1 − t)g(x)

in which

f(x) = x^{2} − 5x + 6 and g(x) = x^{2} − 1.

Show that there is no real path connecting a root of g to a root of f. If there were such a path then it would have to pass through a point when t = 1/2 of the form (1/2, x) where h(1/2, x) = 0. However

h(1/2, x) = \frac{1}{2}(x^{2} − 5x + 6) + \frac{1}{2}(x^{2} − 1) = \frac{1}{2}(2x^{2} − 5x + 5) = 0

implies by the quadratic formula that

x = \frac{5 ± \sqrt{25 − 40}}{4} ∉ R.

Therefore there is no real path connecting a roots of g to a root of f.