1. [Kincaid and Cheney Problem 5.

2#1] Find the Schur factorization of





4 7
3 8
.
and A2 =
A1 =
1 2
2 3
This problem can be worked by hand. For accuracy we use Maple to perform each step in
the proof on Schurs theorem. The Maple script is
1
2

# Kincaid and Cheney Problem 5.2#1

# Schurs factorization for 2x2 matrix

3
4

# Written December 2012 by Eric Olson for Math 701

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

restart;
kernelopts(printbytes=false);
with(LinearAlgebra):
A:=Matrix(A0);
Lambda,X:=Eigenvectors(A);
lambda:=Lambda[1];
x:=X[1..2,1];
x:=x/Norm(x,2);
e1:=Vector([1,0]);
beta:=x[1]/abs(x[1]);
y:=beta*e1;
alpha:=simplify(rationalize(sqrt(2)/Norm(x-y,2)));
v:=simplify(alpha*(x-y));
U:=simplify(IdentityMatrix(2)-v.HermitianTranspose(v));
T:=simplify(U.A.HermitianTranspose(U));
U:=evalf(U);
T:=evalf(T);

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

matprint:=proc(outfile,M)
local fd,i,j;
fd:=open(outfile,WRITE);
for i from 1 to 2 do
for j from 1 to 2 do
if(type(M[i,j],realcons)) then
fprintf(fd,"%g",M[i,j])
elif(type(M[i,j],imaginary)) then
fprintf(fd,"%gI",Im(M[i,j]))
else
fprintf(fd,"%Zg",M[i,j])
fi;
if(j<2) then
fprintf(fd,"&")
else
fprintf(fd,"\\cr\n")
fi
od
od;
close(fd)
end proc;

45
46
47

matprint("U.tex",U);
matprint("T.tex",T);

This script is run with the commands

maple -DA0="[[3,8],[-2,3]]" schur.mpl >schurA1.out
maple -DA0="[[4,7],[1,12]]" schur.mpl >schurA2.out

The output indicates that A1 has a factorization where



0.894427
U1
0.447214i
and

3 + 4i
T1
0

and that A2 has a factorization where


0.622704
U2
0.782457
and

12.7958
T2
0

0.447214i
0.894427
6
3 4i

0.782457
0.622704

6
.
3.20417

On the following pages is the output from Maple.

|\^/|
Maple 9.5 (IBM INTEL LINUX)
._|\|
|/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2004
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
|
Type ? for help.
# Kincaid and Cheney Problem 5.2#1
# Schurs factorization for 2x2 matrix
>
# Written December 2012 by Eric Olson for Math 701
>
> restart;
> kernelopts(printbytes=false);
true
> with(LinearAlgebra):
> A:=Matrix([[3,8],[-2,3]]);
[ 3
A := [
[-2

8]
]
3]

> Lambda,X:=Eigenvectors(A);
[3 + 4 I] [-2 I
Lambda, X := [
], [
[3 - 4 I] [ 1
> lambda:=Lambda[1];
lambda := 3 + 4 I
> x:=X[1..2,1];
[-2 I]
x := [
]
[ 1 ]
> x:=x/Norm(x,2);
[
1/2]
[-2/5 I 5
]
[
]
x := [
1/2
]
[
5
]
[
---]
[
5
]
> e1:=Vector([1,0]);
[1]
e1 := [ ]
[0]
> beta:=x[1]/abs(x[1]);
beta := -I
> y:=beta*e1;
[-I]
y := [ ]
[0 ]

2 I]
]
1 ]

> alpha:=simplify(rationalize(sqrt(2)/Norm(x-y,2)));
1/2
1/2 1/2
1/2
2
(50 - 20 5
)
(5 + 2 5
)
alpha := ----------------------------------10
> v:=simplify(alpha*(x-y));
[
1/2
1/2 1/2
]
[
1/10 I 2
(50 - 20 5
)
]
[
]
v := [ 1/2
1/2 1/2
1/2
1/2]
[2
(50 - 20 5
)
(5 + 2 5
) 5
]
[----------------------------------------]
[
50
]
> U:=simplify(IdentityMatrix(2)-v.HermitianTranspose(v));
[
1/2
]
[ 2 5
1/2]
[ ------1/5 I 5
]
[
5
]
U := [
]
[
1/2 ]
[
1/2
2 5
]
[1/5 I 5
- ------ ]
[
5
]
> T:=simplify(U.A.HermitianTranspose(U));
[3 + 4 I
T := [
[
0

-6

]
]
3 - 4 I]

> U:=evalf(U);
[ 0.8944271908
U := [
[0.4472135954 I

-0.4472135954 I]
]
-0.8944271908 ]

> T:=evalf(T);
[3. + 4. I
T := [
[
0.

-6.

]
]
3. - 4. I]

>
> matprint:=proc(outfile,M)
>
local fd,i,j;
>
fd:=open(outfile,WRITE);
>
for i from 1 to 2 do
>
for j from 1 to 2 do
>
if(type(M[i,j],realcons)) then
>
fprintf(fd,"%g",M[i,j])
>
elif(type(M[i,j],imaginary)) then
>
fprintf(fd,"%gI",Im(M[i,j]))
>
else
>
fprintf(fd,"%Zg",M[i,j])
>
fi;
>
if(j<2) then

>
fprintf(fd,"&")
>
else
>
fprintf(fd,"\\cr\n")
>
fi
>
od
>
od;
>
close(fd)
> end proc;
matprint := proc(outfile, M)
local fd, i, j;
fd := open(outfile, WRITE);
for i to 2 do for j to 2 do
if type(M[i, j], realcons) then fprintf(fd, "%g", M[i, j])
elif type(M[i, j], imaginary) then
fprintf(fd, "%gI", Im(M[i, j]))
else fprintf(fd, "%Zg", M[i, j])
end if;
if j < 2 then fprintf(fd, "&")
else fprintf(fd, "\cr
")
end if
end do
end do;
close(fd)
end proc
>
> matprint("U.tex",U);
> matprint("T.tex",T);
> quit
bytes used=5311532, alloc=3210676, time=0.73

|\^/|
Maple 9.5 (IBM INTEL LINUX)
._|\|
|/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2004
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
|
Type ? for help.
# Kincaid and Cheney Problem 5.2#1
# Schurs factorization for 2x2 matrix
>
# Written December 2012 by Eric Olson for Math 701
>
> restart;
> kernelopts(printbytes=false);
true
> with(LinearAlgebra):
> A:=Matrix([[4,7],[1,12]]);
[4
A := [
[1
> Lambda,X:=Eigenvectors(A);
[
1/2]
[8 + 23
]
Lambda, X := [
],
[
1/2]
[8 - 23
]

7]
]
12]

[
7
[--------[
1/2
[4 + 23
[
[
1

7
]
---------]
1/2]
4 - 23
]
]
1
]

> lambda:=Lambda[1];
1/2
lambda := 8 + 23
> x:=X[1..2,1];
[
7
]
[---------]
x := [
1/2]
[4 + 23
]
[
]
[
1
]
> x:=x/Norm(x,2);
[
7
]
[---------------------------------]
[/
49
\1/2
1/2 ]
[|1 + ------------|
(4 + 23
)]
[|
1/2 2|
]
[\
(4 + 23
) /
]
x := [
]
[
1
]
[
--------------------]
[
/
49
\1/2
]
[
|1 + ------------|
]
[
|
1/2 2|
]
[
\
(4 + 23
) /
]

> e1:=Vector([1,0]);
[1]
e1 := [ ]
[0]
> beta:=x[1]/abs(x[1]);
beta := 1
> y:=beta*e1;
[1]
y := [ ]
[0]
> alpha:=simplify(rationalize(sqrt(2)/Norm(x-y,2)));
1/2 3/4
1/2 1/2
1/2
1/2 1/4
1/2
1/2
16 (11 + 23
)
%1
2
23
8 (11 + 23
)
%1
23
alpha := - ----------------------------------- - ----------------------------49
7
1/2 3/4
1/2 1/2
1/2
1/2 1/4
78 (11 + 23
)
%1
2
39 %1
(11 + 23
)
+ ----------------------------- + -----------------------49
7
1/2
1/2 1/2
%1 := -14 2
+ 8 (11 + 23
)
> v:=simplify(alpha*(x-y));
v :=
[
1/2 3/4
1/2 1/2
1/2 3/4
1/2 1/2
1/2]
[ 11 (11 + 23
)
%1
2
(11 + 23
)
%1
2
23
]
[- ----------------------------- + --------------------------------]
[
392
392
]
[
1/2 3/4
1/2 1/2
1/2
1/2 1/4
1/2
1/2
[15 (11 + 23
)
%1
2
23
(11 + 23
)
%1
23
[----------------------------------- + --------------------------[
392
7
1/2 3/4
1/2 1/2
1/2
1/2 1/4]
67 (11 + 23
)
%1
2
4 %1
(11 + 23
)
]
- ----------------------------- - -----------------------]
392
7
]
1/2
1/2 1/2
%1 := -14 2
+ 8 (11 + 23
)
> U:=simplify(IdentityMatrix(2)-v.HermitianTranspose(v));
U :=
[
1/2
1/2 1/2
1/2 1/2 1/2
1/2
[11 2
(11 + 23
)
(11 + 23
)
2
23
[----------------------- - -------------------------- ,
[
56
56

1/2 1/2 1/2

1/2
1/2
1/2 1/2]
(11 + 23
)
2
23
3 2
(11 + 23
)
]
-------------------------- + ----------------------]
56
56
]
[
1/2 1/2 1/2
1/2
1/2
1/2 1/2
[(11 + 23
)
2
23
3 2
(11 + 23
)
[-------------------------- + ---------------------- ,
[
56
56
1/2
1/2 1/2
1/2 1/2 1/2
1/2]
11 2
(11 + 23
)
(11 + 23
)
2
23
]
- ----------------------- + --------------------------]
56
56
]
> T:=simplify(U.A.HermitianTranspose(U));
[
1/2
[8 + 23
T := [
[
[
0

]
]
]
1/2]
8 - 23
]
-6

> U:=evalf(U);
[0.6227042073
U := [
[0.7824573279

0.7824573279 ]
]
-0.6227042073]

> T:=evalf(T);
[12.79583152
T := [
[
0.

]
]
3.204168477]

>
> matprint:=proc(outfile,M)
>
local fd,i,j;
>
fd:=open(outfile,WRITE);
>
for i from 1 to 2 do
>
for j from 1 to 2 do
>
if(type(M[i,j],realcons)) then
>
fprintf(fd,"%g",M[i,j])
>
elif(type(M[i,j],imaginary)) then
>
fprintf(fd,"%gI",Im(M[i,j]))
>
else
>
fprintf(fd,"%Zg",M[i,j])
>
fi;
>
if(j<2) then
>
fprintf(fd,"&")
>
else
>
fprintf(fd,"\\cr\n")
>
fi
>
od
>
od;
>
close(fd)
> end proc;

-6.

matprint := proc(outfile, M)
local fd, i, j;
fd := open(outfile, WRITE);
for i to 2 do for j to 2 do
if type(M[i, j], realcons) then fprintf(fd, "%g", M[i, j])
elif type(M[i, j], imaginary) then
fprintf(fd, "%gI", Im(M[i, j]))
else fprintf(fd, "%Zg", M[i, j])
end if;
if j < 2 then fprintf(fd, "&")
else fprintf(fd, "\cr
")
end if
end do
end do;
close(fd)
end proc
>
> matprint("U.tex",U);
> matprint("T.tex",T);
> quit
bytes used=35862804, alloc=6814496, time=5.74

2. [Kincaid and Cheney Problem 5.2#16] Prove that if (I vv )x = y for three vectors
v, x and y then hx, yi is real.
Computing yields
hx, yi = hx, (I vv )xi = hx, x (
v x)vi = kxk2 (v x
)hx, vi
= kxk2 hv, cihx, vi = kxk2 |hx, vi|2 R.

10

3. [Kincaid and Cheney Problem 5.2#34] Let U = I uu where u is a given vector.

Find all the complex values of for which U is unitary.
Write = a + ib. Since U is unitary then
)
I = U U = (I uu )(I uu ) = (I uu )(I uu

= I ( + )uu
+ ||2 uu uu = I 2auu + (a2 + b2 )|u|2 uu
implies 2a + (a2 + b2 )|u|2 = 0. Completing the square yields


1 2
1
2
+
b
=
|u|2
|u|4

which is a circle of radius 1/|u|2 centered at 1/|u|2 . Therefore any on the circle
n
ensures that U is unitary.


1 o
1 

C :  2  =
|u|
|u|2

11

4. [Kincaid and Cheney Problem 5.2#39] Prove that det(I + xx ) = 1 + x x.

Let x Rn . Define

= kxk

1
x1 /|x1 |

if x1 = 0
otherwise

and let y = e1 . Then

kyk2 = ||2 = kxk2

and

hx, yi = x1 = |x1 | R

imply by The Second Lemma on Unitary Matrices from page 267 of Kincaid and Cheney
that there exists a unitary matrix U of the form I vv such that U x = y. Now


det(I + xx ) = det U (I + xx )U = det(U U + U xx U )

= det(I + yy ) = det(I + ||2 e1 e1 ) = 1 + ||2 = 1 + x x.

12

5. [Kincaid and Cheney Problem 5.3#7] What is the determinant of a unitary matrix?
What is the determinant of an orthogonal matrix?
Let U Cnn be a unitary matrix. By definition U U = I. Therefore
1 = det(I) = det(U U ) = det(U ) det(U ) = det(U )det(U T ) = | det(U )|2 .
It follows that det(U ) = ei for some R. Conversely, suppose R. Let U Cnn be
the diagonal matrix such that u11 = ei and ujj = 1 for j > 1. Then det U = ei . Moreover
U U is the diagonal matrix with diagonal entries [U U ]jj = ujj u
jj = 1 for j = 1, . . . , n.

Thus U U = I and so U is unitary.

Now let O Rnn be an orthogonal matrix. By definition OOT = I. Therefore
1 = det(I) = det(OOT ) = det(O) det(OT ) = det(O)2 .
Since det(O) R it follows that det(O) = 1. Conversely, given = 1 define O Rnn
to be the diagonal matrix such that o11 = and ojj = 1 for j > 1. Clearly OOT = I and
so O is orthogonal.

13

6. [Kincaid and Cheney Problem 5.3#11] A matrix A such that A2 = I is called an

involution or is said to be involutory. Find the necessary and sufficient conditions of
u and v in order that I uv be an involution.

Let A Cnn be such that A2 = I. Then

I = A2 = (I uv )2 = I 2uv + uv uv
= I 2uv + (
v u)uv = I 2uv + hu, viuv
implies hu, vi = 2 or that u = 0 or v = 0. On the other hand, if either hu, vi = 2 or u = 0
or v = 0, then the same equality shows that I uv is an involution.
Therefore, it is a necessary and sufficient condition that either hu, vi = 2 or u = 0 or
v = 0 in order for I uv to be an involution.

14

7. [Kincaid and Cheney Problem 6.2#4] Prove if f is a polynomial of degree k, then for
n > k that f [x0 , x1 , . . . , xn ] = 0.
Let n > k and pn (x) be the interpolating polynomial of the form
pn (x) =

n
X

ci

i=0

i1
Y

(x xj )

j=0

passing through the points jk , f (xj ) where j = 0, . . . , n. By definition
ci = f [x0 , x1 , . . . , xi ].
By the Theorem of Polynomial Interpolation Error from page 315 in Kincaid and Cheney
we have that
n
Y
1
f (n+1) () (x xi ).
f (x) pn (x) =
(n + 1)!
i=0
Since f (x) is a polynomial of degree k then f (n+1) (x) = 0. Therefore f (x) = pn (x). This
implies that pn (x) is a polynomial of degree k. Thus cn = f [x0 , x1 , . . . , xn ] = 0.

15

8. [Kincaid and Cheney Problem 6.2#25] Establish an iteration method for solving
f (x) = 0 as follows:
Let qn be the
quadratic interpolating



polynomial through the
points xn , f (xn ) , xn1 , f (xn1 ) and xn2 , f (xn2 ) and define xn+1 to be the
zero of qn closest to xn .
As an aside this method was established by D.E. M
uller in 1956 and is discussed in detail
in A Survey of Numerical Mathematics by D. Young and R. Gregory.
Newtons divided difference formula yields
qn (x) = f [xn ] + f [xn , xn1 ](x xn ) + f [xn , xn1 , xn2 ](x xn )(x xn1 )
where the f [x], f [x, y] and f [x, y, z] are given recursively as
f [x] = f (x)
f [x, y] = (f [y] f [x])/(y x)

f [x, y, z] = (f [y, z] f [x, y])/(z x).

Writing c0 = f [xn ], c1 = f [xn , xn1 ] and c2 = f [xn , xn1 , xn2 ] we see that


qn (x) = c0 + c1 + c2 (xn xn1 ) (x xn ) + c2 (x xn )2
= a(x xn )2 + b(x xn ) + c

where a = c2 , b = c1 + c2 (xn xn1 ) and c = c0 . Now use the quadratic formula to solve
for xn+1 . First define
M1 = b +

p
b2 4ac

Then
xn+1 = xn +

and


2c/M2
2c/M1

M2 = b

p
b2 4ac.

if |M1 | < |M2 |

otherwise.

Note that we have written the quadratic formula as

2c

b b2 4ac
with the square root in the denominator tominimize loss of precision. For example, if
|M1 | is small, this means the terms b and b2 4ac nearly cancelled. This cancellation
causes the value of M1 to lose precision. Therefore, we compute xn+1 using M2 . Similarly,
if M2 is small, we compute xn+1 using M1 .

16

9. [Kincaid and Cheney Problem 6.2#1] For n = 5, 10 and 15 find the Newton interpolating polynomial pn for the function f (x) = 1/(1 + x2 ) on the interval [5, 5]. Use
equally spaced nodes. In each case, compute f (x) pn (x) for 30 equally spaced points
in [5, 5] in order to see the divergence of pn from f .
The program is
1
2
3

/*
Kincaid and Cheney Programming Problem 6.2#1
Polynomial approximation of degree n

4
5
6
7
8

Written December 2012 by Eric Olson for Math 701

*/
#include <stdio.h>
#include <stdlib.h>

9
10
11
12
13
14
15
16
17
18
19
20
21
22

void divdif(int n,double x[n+1],double c[n+1]){

for(int k=1;k<=n;k++) for(int i=n;i>=k;i--) {
c[i]=(c[i-1]-c[i])/(x[i-k]-x[i]);
}
}
double interp(int n,double x[n+1],double c[n+1],double z){
double y=0;
for(int i=n;i>0;i--) y=(z-x[i-1])*(y+c[i]);
return y+c[0];
}
double f(double x){
return 1/(1+x*x);
}

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

int main(){
int n;
scanf("%d",&n);
printf("# Kincaid and Cheney Programming Problem 6.2#1\n"
"# Polynomial approximation of degree %d\n\n",n);
double x[n+1],c[n+1];
double a=-5,b=5,h=(b-a)/n;
for(int i=0;i<=n;i++){
x[i]=a+h*i;
c[i]=f(x[i]);
}
divdif(n,x,c);
h=(b-a)/30;
printf("#%17s %18s %18s %18s\n","z","f(z)","Pn(z)","f(z)-Pn(z)");
for(int i=0;i<=30;i++){
double z=a+h*i, y=f(z), p=interp(n,x,c,z);
printf("%18.10e %18.10e %18.10e %18.10e\n",z,y,p,y-p);
}
exit(0);
}

and my plotting script is

1
2
3
4
5
6

set terminal pstex

set key outside
set key width 1
set samples 1000
set size 1.1,0.75
plot datafile using 1:4 ti "$f(x)-p_n(x)$"

17

Output for n = 5
0.5

f (x) pn (x)

0.4
0.3
0.2
0.1
0
-0.1
-0.2
-6

-4

-2

# Kincaid and Cheney Programming Problem 6.2#1

# Polynomial approximation of degree 5
#
z
-5.0000000000e+00
-4.6666666667e+00
-4.3333333333e+00
-4.0000000000e+00
-3.6666666667e+00
-3.3333333333e+00
-3.0000000000e+00
-2.6666666667e+00
-2.3333333333e+00
-2.0000000000e+00
-1.6666666667e+00
-1.3333333333e+00
-1.0000000000e+00
-6.6666666667e-01
-3.3333333333e-01
-2.7755575616e-16
3.3333333333e-01
6.6666666667e-01
1.0000000000e+00
1.3333333333e+00
1.6666666667e+00
2.0000000000e+00
2.3333333333e+00
2.6666666667e+00
3.0000000000e+00
3.3333333333e+00
3.6666666667e+00
4.0000000000e+00
4.3333333333e+00
4.6666666667e+00
5.0000000000e+00

f(z)
3.8461538462e-02
4.3902439024e-02
5.0561797753e-02
5.8823529412e-02
6.9230769231e-02
8.2568807339e-02
1.0000000000e-01
1.2328767123e-01
1.5517241379e-01
2.0000000000e-01
2.6470588235e-01
3.6000000000e-01
5.0000000000e-01
6.9230769231e-01
9.0000000000e-01
1.0000000000e+00
9.0000000000e-01
6.9230769231e-01
5.0000000000e-01
3.6000000000e-01
2.6470588235e-01
2.0000000000e-01
1.5517241379e-01
1.2328767123e-01
1.0000000000e-01
8.2568807339e-02
6.9230769231e-02
5.8823529412e-02
5.0561797753e-02
4.3902439024e-02
3.8461538462e-02

Pn(z)
3.8461538462e-02
-2.8323836657e-02
-5.4605887939e-02
-4.8076923077e-02
-1.5859449193e-02
3.5493827160e-02
1.0000000000e-01
1.7224596391e-01
2.4738841406e-01
3.2115384615e-01
3.8983855651e-01
4.5030864198e-01
5.0000000000e-01
5.3691832858e-01
5.5963912631e-01
5.6730769231e-01
5.5963912631e-01
5.3691832858e-01
5.0000000000e-01
4.5030864198e-01
3.8983855651e-01
3.2115384615e-01
2.4738841406e-01
1.7224596391e-01
1.0000000000e-01
3.5493827160e-02
-1.5859449193e-02
-4.8076923077e-02
-5.4605887939e-02
-2.8323836657e-02
3.8461538462e-02

18

f(z)-Pn(z)
0.0000000000e+00
7.2226275682e-02
1.0516768569e-01
1.0690045249e-01
8.5090218424e-02
4.7074980179e-02
0.0000000000e+00
-4.8958292680e-02
-9.2216000262e-02
-1.2115384615e-01
-1.2513267415e-01
-9.0308641975e-02
-1.1102230246e-16
1.5538936372e-01
3.4036087369e-01
4.3269230769e-01
3.4036087369e-01
1.5538936372e-01
2.2204460493e-16
-9.0308641975e-02
-1.2513267415e-01
-1.2115384615e-01
-9.2216000262e-02
-4.8958292680e-02
2.7755575616e-17
4.7074980179e-02
8.5090218424e-02
1.0690045249e-01
1.0516768569e-01
7.2226275682e-02
5.3429483060e-16

Output for n = 10
0.5

f (x) pn (x)

0
-0.5
-1
-1.5
-2
-6

-4

-2

# Kincaid and Cheney Programming Problem 6.2#1

# Polynomial approximation of degree 10
#
z
-5.0000000000e+00
-4.6666666667e+00
-4.3333333333e+00
-4.0000000000e+00
-3.6666666667e+00
-3.3333333333e+00
-3.0000000000e+00
-2.6666666667e+00
-2.3333333333e+00
-2.0000000000e+00
-1.6666666667e+00
-1.3333333333e+00
-1.0000000000e+00
-6.6666666667e-01
-3.3333333333e-01
-2.7755575616e-16
3.3333333333e-01
6.6666666667e-01
1.0000000000e+00
1.3333333333e+00
1.6666666667e+00
2.0000000000e+00
2.3333333333e+00
2.6666666667e+00
3.0000000000e+00
3.3333333333e+00
3.6666666667e+00
4.0000000000e+00
4.3333333333e+00
4.6666666667e+00
5.0000000000e+00

f(z)
3.8461538462e-02
4.3902439024e-02
5.0561797753e-02
5.8823529412e-02
6.9230769231e-02
8.2568807339e-02
1.0000000000e-01
1.2328767123e-01
1.5517241379e-01
2.0000000000e-01
2.6470588235e-01
3.6000000000e-01
5.0000000000e-01
6.9230769231e-01
9.0000000000e-01
1.0000000000e+00
9.0000000000e-01
6.9230769231e-01
5.0000000000e-01
3.6000000000e-01
2.6470588235e-01
2.0000000000e-01
1.5517241379e-01
1.2328767123e-01
1.0000000000e-01
8.2568807339e-02
6.9230769231e-02
5.8823529412e-02
5.0561797753e-02
4.3902439024e-02
3.8461538462e-02

Pn(z)
3.8461538462e-02
1.9446662941e+00
1.0042671440e+00
5.8823529412e-02
-2.5560610028e-01
-1.3136348377e-01
1.0000000000e-01
2.3655582000e-01
2.4523053137e-01
2.0000000000e-01
1.9862072522e-01
3.0030433349e-01
5.0000000000e-01
7.3724559916e-01
9.2749142537e-01
1.0000000000e+00
9.2749142537e-01
7.3724559916e-01
5.0000000000e-01
3.0030433349e-01
1.9862072522e-01
2.0000000000e-01
2.4523053137e-01
2.3655582000e-01
1.0000000000e-01
-1.3136348377e-01
-2.5560610028e-01
5.8823529412e-02
1.0042671440e+00
1.9446662941e+00
3.8461538462e-02

19

f(z)-Pn(z)
0.0000000000e+00
-1.9007638551e+00
-9.5370534628e-01
0.0000000000e+00
3.2483686951e-01
2.1393229111e-01
0.0000000000e+00
-1.1326814877e-01
-9.0058117580e-02
0.0000000000e+00
6.6085157131e-02
5.9695666510e-02
0.0000000000e+00
-4.4937906849e-02
-2.7491425367e-02
1.1102230246e-16
-2.7491425367e-02
-4.4937906849e-02
1.1102230246e-16
5.9695666510e-02
6.6085157131e-02
1.0269562978e-15
-9.0058117580e-02
-1.1326814877e-01
2.5118795932e-15
2.1393229111e-01
3.2483686951e-01
6.9527716917e-15
-9.5370534628e-01
-1.9007638551e+00
-4.7066517350e-14

Output for n = 15
0.4

f (x) pn (x)

0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.2
-1.4
-1.6
-6

-4

-2

# Kincaid and Cheney Programming Problem 6.2#1

# Polynomial approximation of degree 15
#
z
-5.0000000000e+00
-4.6666666667e+00
-4.3333333333e+00
-4.0000000000e+00
-3.6666666667e+00
-3.3333333333e+00
-3.0000000000e+00
-2.6666666667e+00
-2.3333333333e+00
-2.0000000000e+00
-1.6666666667e+00
-1.3333333333e+00
-1.0000000000e+00
-6.6666666667e-01
-3.3333333333e-01
-2.7755575616e-16
3.3333333333e-01
6.6666666667e-01
1.0000000000e+00
1.3333333333e+00
1.6666666667e+00
2.0000000000e+00
2.3333333333e+00
2.6666666667e+00
3.0000000000e+00
3.3333333333e+00
3.6666666667e+00
4.0000000000e+00
4.3333333333e+00
4.6666666667e+00
5.0000000000e+00

f(z)
3.8461538462e-02
4.3902439024e-02
5.0561797753e-02
5.8823529412e-02
6.9230769231e-02
8.2568807339e-02
1.0000000000e-01
1.2328767123e-01
1.5517241379e-01
2.0000000000e-01
2.6470588235e-01
3.6000000000e-01
5.0000000000e-01
6.9230769231e-01
9.0000000000e-01
1.0000000000e+00
9.0000000000e-01
6.9230769231e-01
5.0000000000e-01
3.6000000000e-01
2.6470588235e-01
2.0000000000e-01
1.5517241379e-01
1.2328767123e-01
1.0000000000e-01
8.2568807339e-02
6.9230769231e-02
5.8823529412e-02
5.0561797753e-02
4.3902439024e-02
3.8461538462e-02

Pn(z)
3.8461538462e-02
1.6149814985e+00
5.0561797753e-02
-1.5893928816e-01
6.9230769231e-02
1.3917382108e-01
1.0000000000e-01
9.9622123022e-02
1.5517241379e-01
2.1502247843e-01
2.6470588235e-01
3.4583594891e-01
5.0000000000e-01
7.1094460164e-01
9.0000000000e-01
9.7624707634e-01
9.0000000000e-01
7.1094460164e-01
5.0000000000e-01
3.4583594891e-01
2.6470588235e-01
2.1502247843e-01
1.5517241379e-01
9.9622123022e-02
1.0000000000e-01
1.3917382108e-01
6.9230769231e-02
-1.5893928816e-01
5.0561797753e-02
1.6149814985e+00
3.8461538464e-02

20

f(z)-Pn(z)
0.0000000000e+00
-1.5710790595e+00
0.0000000000e+00
2.1776281757e-01
0.0000000000e+00
-5.6605013742e-02
0.0000000000e+00
2.3665548211e-02
0.0000000000e+00
-1.5022478430e-02
0.0000000000e+00
1.4164051091e-02
0.0000000000e+00
-1.8636909330e-02
0.0000000000e+00
2.3752923656e-02
0.0000000000e+00
-1.8636909330e-02
4.4408920985e-16
1.4164051091e-02
3.0531133177e-15
-1.5022478430e-02
7.6882944455e-15
2.3665548211e-02
4.6976311729e-14
-5.6605013742e-02
8.2156503822e-14
2.1776281757e-01
-4.5050768671e-13
-1.5710790595e+00
-2.2520735277e-12

10. [Kincaid and Cheney Problem 6.2#2] Program M

ullers method for finding roots as
described in Problem 6.2.25 above and test it.
The program is
1
2
3

/*
Kincaid and Cheney Programming Problem 6.2#2
Solve F=0 by Mullers Method

4
5
6
7
8
9
10

Written December 2012 by Eric Olson for Math 701

*/
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

#define N 40
#define STR1(s) #s
#define STR2(s) STR1(s)
typedef double complex Complex;
Complex f(Complex z){
return F;
}
Complex muller(Complex z[3],Complex w[3]){
Complex c[3];
for(int k=0;k<3;k++) c[k]=w[k];
for(int k=1;k<3;k++) for(int i=0;i<3-k;i++)
c[i]=(c[i+1]-c[i])/(z[i+k]-z[i]);
Complex b=c[1]+c[0]*(z[2]-z[1]);
Complex d=csqrt(b*b-4*c[2]*c[0]);
Complex M1=-b+d,M2=-b-d;
if(cabs(M1)<cabs(M2)) return z[2]+2*c[2]/M2;
else return z[2]+2*c[2]/M1;
}
int main(){
Complex w[N],z[N]={0.25,0.5,0.75};
printf("# Kincaid and Cheney Programming Problem 6.2#2\n"
"# Solve %s=0 by Mullers Method\n\n",STR2(F));
printf("#%2s %17s %17s %17s %17s\n",
"i","real(z[i])","imag(z[i])","real(w[i])","imag(w[i])");
for(int i=0;i<3;i++) {
w[i]=f(z[i]);
printf("%3d %17.10e %17.10e %17.10e %17.10e\n",
i,creal(z[i]),cimag(z[i]),creal(w[i]),cimag(w[i]));
}
for(int i=3;i<N;i++){
z[i]=muller(&z[i-3],&w[i-3]);
w[i]=f(z[i]);
printf("%3d %17.10e %17.10e %17.10e %17.10e\n",
i,creal(z[i]),cimag(z[i]),creal(w[i]),cimag(w[i]));
if(z[i]==z[i-1]) exit(0);
}
printf("# Didnt converge!\n");
exit(0);
}

The program can be compiled to solve roots for different functions as

gcc -std=gnu99 -DF="z*z+1" -o prog10a prog10.c -lm
gcc -std=gnu99 -DF="cexp(z)+z" -o prog10b prog10.c -lm
gcc -std=gnu99 -DF="z*z*z" -o prog10c prog10.c -lm

21

The output is
# Kincaid and Cheney Programming Problem 6.2#2
# Solve z*z+1=0 by Mullers Method
# i
0
1
2
3
4

real(z[i])
imag(z[i])
2.5000000000e-01 0.0000000000e+00
5.0000000000e-01 0.0000000000e+00
7.5000000000e-01 0.0000000000e+00
0.0000000000e+00 -1.0000000000e+00
0.0000000000e+00 -1.0000000000e+00

real(w[i])
imag(w[i])
1.0625000000e+00 0.0000000000e+00
1.2500000000e+00 0.0000000000e+00
1.5625000000e+00 0.0000000000e+00
0.0000000000e+00 -0.0000000000e+00
0.0000000000e+00 -0.0000000000e+00

and
# Kincaid and Cheney Programming Problem 6.2#2
# Solve cexp(z)+z=0 by Mullers Method
# i
0
1
2
3
4
5
6
7
8
9
10
11

real(z[i])
imag(z[i])
real(w[i])
imag(w[i])
2.5000000000e-01 0.0000000000e+00 1.5340254167e+00 0.0000000000e+00
5.0000000000e-01 0.0000000000e+00 2.1487212707e+00 0.0000000000e+00
7.5000000000e-01 0.0000000000e+00 2.8670000166e+00 0.0000000000e+00
-1.1085844094e+00 -7.3847789802e-02 -7.7945808573e-01 -9.8197320828e-02
-4.8978113442e-01 7.9411254946e-03 1.2296003658e-01 1.2807082314e-02
-5.6306863605e-01 9.0034282357e-04 6.3900509273e-03 1.4130510042e-03
-5.6713463412e-01 4.0762660781e-06 1.3565659639e-05 6.3881130461e-06
-5.6714329055e-01 -1.3198551104e-10 -2.1671109351e-10 -2.0684020805e-10
-5.6714329041e-01 4.5536763180e-19 0.0000000000e+00 7.1362632885e-19
-5.6714329041e-01 -3.2710376473e-27 0.0000000000e+00 -5.1261847016e-27
-5.6714329041e-01 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
-5.6714329041e-01 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00

and
# Kincaid and Cheney Programming Problem 6.2#2
# Solve z*z*z=0 by Mullers Method
# i
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

real(z[i])
2.5000000000e-01
5.0000000000e-01
7.5000000000e-01
2.2916666667e-01
1.8242733623e-01
1.1320131803e-01
4.6023561609e-02
1.1607687387e-02
-9.9103148828e-03
-2.0770500757e-02
-2.4036155562e-02
-2.2906696118e-02
-1.9475254483e-02
-1.5209134709e-02
-1.1002140304e-02
-7.3440644873e-03
-4.4365844552e-03
-2.2967330808e-03
-8.4064006431e-04

imag(z[i])
0.0000000000e+00
0.0000000000e+00
0.0000000000e+00
-9.9913156736e-02
-1.3330418665e-01
-1.3486746718e-01
-1.1245265817e-01
-9.1451368313e-02
-6.7765499686e-02
-4.5953326202e-02
-2.8722430507e-02
-1.5664398875e-02
-6.5590742174e-03
-7.7341135865e-04
2.5012212938e-03
3.9925003511e-03
4.3213996534e-03
3.9724821450e-03
3.2941275908e-03

real(w[i])
1.5625000000e-02
1.2500000000e-01
4.2187500000e-01
5.1721643519e-03
-3.6540716340e-03
-4.7265130502e-03
-1.6485010517e-03
-2.8967355195e-04
1.3555600737e-04
1.2262303100e-04
4.5601328675e-05
4.8425545711e-06
-4.8731210290e-06
-3.4908505335e-06
-1.1252853490e-06
-4.4909733064e-08
1.6122638091e-07
9.6616349733e-08
2.6771993616e-08

22

imag(w[i])
0.0000000000e+00
0.0000000000e+00
0.0000000000e+00
-1.4744129032e-02
-1.0940166999e-02
-2.7316571328e-03
7.0745042963e-04
7.2787401925e-04
2.9122358509e-04
3.7565331570e-05
-2.6086639360e-05
-2.0814489640e-05
-7.1811050500e-06
-5.3624876460e-07
8.9264877630e-07
5.8236921990e-07
1.7447802315e-07
1.7611658246e-10
-2.8761849654e-08

19 6.1640117010e-05
20 5.4938226284e-04
21 7.4979191331e-04
22 7.6791662094e-04
23 6.8376017065e-04
24 5.5377150915e-04
25 4.1455403644e-04
26 2.8733132646e-04
27 1.8228363969e-04
28 1.0228172709e-04
29 4.5832133906e-05
30 9.2287945586e-06
31 -1.1994753859e-05
32 -2.2175025194e-05
33 -2.5045773302e-05
34 -2.3552319884e-05
35 -1.9846253763e-05
36 -1.5379790863e-05
37 -1.1043277129e-05
38 -7.3095232165e-06
39 -4.3642015246e-06
# Didnt converge!

2.5180731351e-03
1.7834926724e-03
1.1616988086e-03
6.7801675860e-04
3.2942105091e-04
9.7697631033e-05
-4.1474291277e-05
-1.1278707803e-04
-1.3811091152e-04
-1.3508677373e-04
-1.1687225098e-04
-9.2551441047e-05
-6.7859165502e-05
-4.5988948670e-05
-2.8347969344e-05
-1.5190386228e-05
-6.1062869208e-06
-3.7146558373e-07
2.8220586809e-06
4.2292375149e-06
4.4833962082e-06

-1.1722888466e-09
-5.0766862725e-09
-2.6141078567e-09
-6.0621199341e-10
9.7075677322e-11
1.5396420833e-10
6.9103960168e-11
1.2756509906e-11
-4.3741778071e-12
-4.5294192110e-12
-1.7818059998e-12
-2.3636915123e-13
1.6397698034e-13
1.2979526675e-13
4.4669921557e-14
3.2391912834e-15
-5.5969128459e-15
-3.6315378443e-15
-1.0829250043e-15
1.6828000373e-18
1.8005067527e-16

-1.5937624687e-08
-4.0581298001e-09
3.9151528855e-10
8.8778231918e-10
4.2629247760e-10
8.8948204981e-11
-2.1311358469e-11
-2.6500107687e-11
-1.1132754862e-11
-1.7745293966e-12
8.5987234105e-13
7.6912628792e-13
2.8319285408e-13
2.9423185432e-14
-3.0566620755e-14
-2.1773709465e-14
-6.9876355555e-15
-2.6354588463e-16
1.0100088588e-15
6.0224739337e-16
1.6605568481e-16

M
ullers method quickly finds the roots in the first two cases. In the last case, however,
the method converges slowly because the derivative vanishes at the root.

23