\magnification=\magstep1
\vsize=9truein
\hsize=6.5truein
%\parindent=0.8truein
\def\C{{\bf C}}
\def\R{{\bf R}}
\def\N{{\bf N}}
\def\span{{\rm span}}
\font\sc=cmcsc10
\font\smalltt=cmtt10 scaled840
\font\smseven=cmr7 scaled840

%verbatim environment
\def\uncatcodespecials{\def\do##1{\catcode`##1=12 }\dospecials}
\newcount\lineno
%\def\listing#1{\par\begingroup\setupverbatim #1\endgroup}
\def\setupverbatimc{\parindent=1truecm\tt
	\smalltt\baselineskip=10pt plus 1pt minus 1pt
	\def\par{\leavevmode\endgraf}\catcode`\`=\active\obeylines
	\uncatcodespecials\obeyspaces
	\everypar{\global\advance\lineno by 1\llap{\smseven\the\lineno\ \ }}}
\def\setupverbatimb{\parindent=1truecm\global\lineno=0%
	\smalltt\baselineskip=10pt plus 1pt minus 1pt
	\def\par{\leavevmode\endgraf}\catcode`\`=\active\obeylines
	\uncatcodespecials\obeyspaces
	\everypar{\advance\lineno by 1{}}}
\def\setupverbatim{\global\lineno=0\setupverbatimc}

{\obeyspaces\global\let =\ }

\def\beginverbatim{\par\begingroup\setupverbatim\doverbatim}
{\catcode`\|=0\catcode`\\=12\catcode`|`=|active|gdef`{|relax|lq}%
|obeylines|gdef|doverbatim^^M#1\endverbatim{#1|endgroup}}

\def\beginverbatimb{\par\begingroup\setupverbatimb\doverbatim}
\def\listing#1{\medskip
	\begingroup\setupverbatim\input#1\endgroup\medskip
	\noindent\ignorespaces}
\def\listout#1{\medskip
	\begingroup\setupverbatimb\input#1\endgroup\medskip
	\noindent\ignorespaces}
\def\macinput#1{{\catcode`\$=10\catcode`@=3 \input #1 }}

\newcount\qnum\qnum=0
\def\qn{\advance\qnum by1\medskip\par\noindent\item{\the\qnum.}}

\input epsf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\qn [Kincaid and Cheney Problem 5.2\#1] \ Find the Schur factorization
of $$A_1=\left[\matrix{3&8\cr-2&3}\right]\quad\hbox{and}\quad
	A_2=\left[\matrix{4&7\cr-1&2}\right].$$
\medskip\noindent
This problem can be worked by hand.  For accuracy we use Maple to
perform each step in the proof on Schur's theorem.  The Maple script
is
\listing{schur.mpl}
\vfill\eject
\noindent
This script is run with the commands
\medskip
\beginverbatimb
	maple -DA0="[[3,8],[-2,3]]" schur.mpl >schurA1.out
	maple -DA0="[[4,7],[1,12]]" schur.mpl >schurA2.out
\endverbatim
\medskip\noindent
The output indicates that $A_1$ has a factorization where
$$
	U_1\approx\left[\matrix{
\input U1.tex
	}\right]
$$
and
$$
	T_1\approx\left[\matrix{
\input T1.tex
	}\right]
$$
and that $A_2$ has a factorization where
$$
	U_2\approx\left[\matrix{
\input U2.tex
	}\right]
$$
and
$$
	T_2\approx\left[\matrix{
\input T2.tex
	}\right].
$$
On the following pages is the output from Maple.
\vfill\eject
\listout{schurA1.out}
\vfill\eject
\listout{schurA2.out}
\vfill\eject
\qn [Kincaid and Cheney Problem 5.2\#16] \ Prove that if $(I-vv^*)x=y$ for
three vectors $v$, $x$ and $y$ then $\langle x,y\rangle$ is real.
\medskip\noindent
Computing yields
$$\eqalign{
	\langle x,y\rangle
		&=\langle x, (I-vv^*)x\rangle
		=\langle x, x-(\bar v\cdot x)v\rangle
		=\|x\|^2-(v\cdot \bar x)\langle x,v\rangle\cr
		&=\|x\|^2-\langle v,c\rangle \langle x,v\rangle
		=\|x\|^2-|\langle x,v\rangle|^2\in\R.
}$$
\vfill\eject

\qn [Kincaid and Cheney Problem 5.2\#34] \ Let $U=I-\lambda uu^*$ where
$u$ is a given vector.  Find all the complex values of $\lambda$
for which $U$ is unitary.
\medskip\noindent
Write $\lambda=a+ib$.  
Since $U$ is unitary then 
$$\eqalign{
	I&=UU^*=(I-\lambda uu^*)(I-\lambda uu^*)^*
		=(I-\lambda uu^*)(I-\bar\lambda uu^*)\cr
		&=I-(\lambda+\bar\lambda) uu^* +|\lambda|^2uu^*uu^*
		=I-2a uu^* +(a^2+b^2)|u|^2 uu^*
}$$
implies $-2a+(a^2+b^2)|u|^2=0$.  Completing the square yields
%$$
%	a^2-{2\over |u|^2}a+b^2=0
%$$
%and completing the square yields
$$
	\Big(a-{1\over |u|^2}\Big)^2+b^2={1\over |u|^4}
$$
which is a circle of radius $1/|u|^2$ centered at $1/|u|^2$.
Therefore any $\lambda$ on the circle
$$
	\Big\{\,\lambda\in\C: \Big|\lambda-{1\over |u|^2}\Big| = {1\over |u|^2}\Big\}
$$
ensures that $U$ is unitary.
\vfill\eject

\qn [Kincaid and Cheney Problem 5.2\#39] \ Prove that $\det(I+xx^*)=1+x^* x$.
\medskip\noindent
Let $x\in\R^n$.  Define
$$
	\beta=\|x\|\cases{ 1&if $x_1=0$\cr
				x_1/|x_1|& otherwise\cr}
$$
and let $y=\beta e_1$.  Then
$$
	\|y\|^2=|\beta|^2=\|x\|^2
\quad\hbox{and}\quad
	\langle x,y\rangle = x_1 \bar \beta = |x_1|\in\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 $Ux=y$.  Now
$$\eqalign{
	\det(I+xx^*)&=\det\big(U(I+xx^*)U^*\big)
		=\det(UU^*+Uxx^*U^*)\cr
		&=\det(I+yy^*)=\det(I+|\beta|^2 e_1 e_1^*)=1+|\beta|^2=1+x^*x.
}$$
\vfill\eject

\qn [Kincaid and Cheney Problem 5.3\#7] \ What is the determinant
of a unitary matrix?  What is the determinant of an orthogonal matrix?
\medskip\noindent
Let $U\in \C^{n\times n}$ be a unitary matrix.  By definition $U U^*=I$.
Therefore
$$
	1=\det(I)=\det(UU^*)=\det(U)\det(U^*)=\det(U)\overline{\det(U^T)}
	=|\det(U)|^2.
$$
It follows that $\det(U)=e^{i\theta}$ for some $\theta\in\R$.
Conversely, suppose $\theta\in\R$.  
Let $U\in\C^{n\times n}$ be the diagonal matrix
such that $u_{11}=e^{i\theta}$ and $u_{jj}=1$ for $j>1$.
Then $\det U=e^{i\theta}$.  Moreover $UU^*$ is the diagonal
matrix with diagonal entries $[UU^*]_{jj}=u_{jj} \bar u_{jj}=1$
for $j=1,\ldots,n$.  Thus $UU^*=I$ and so $U$ is unitary.

Now let $O\in\R^{n\times n}$ be an orthogonal matrix.  By definition
$O O^T=I$.  Therefore
$$
	1=\det(I)=\det(OO^T)=\det(O)\det(O^T)
	=\det(O)^2.
$$
Since $\det(O)\in\R$ it follows that $\det(O)=\pm 1$.  
Conversely, given $\alpha=\pm 1$ define $O\in\R^{n\times n}$ to be
the diagonal matrix such that $o_{11}=\alpha$ and $o_{jj}=1$ for $j>1$.
Clearly $OO^T=I$ and so $O$ is orthogonal.
\vfill\eject

\qn [Kincaid and Cheney Problem 5.3\#11] \ A matrix $A$ such
that $A^2=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.
\medskip\noindent
Let $A\in\C^{n\times n}$ be such that $A^2=I$.  Then
$$\eqalign{
	I=A^2&=(I-uv^*)^2=I-2uv^*+uv^*uv^*\cr
		&=I-2uv^*+(\bar v\cdot u) uv^*
		=I-2uv^*+\langle u,v\rangle uv^*
}$$
implies $\langle u,v\rangle=2$ or that $u=0$ or $v=0$.  On the other hand,
if either $\langle u,v\rangle=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
$\langle u,v\rangle=2$ or $u=0$ or $v=0$ in order 
for $I-uv^*$ to be an involution.

\vfill\eject

\qn [Kincaid and Cheney Problem 6.2\#4] \ Prove if $f$ is
a polynomial of degree $k$, then for $n>k$ that
$f[x_0,x_1,\ldots,x_n]=0.$
\medskip\noindent
Let $n>k$ and $p_n(x)$ be the interpolating polynomial of the form
$$
	p_n(x)=\sum_{i=0}^n c_i\prod_{j=0}^{i-1}(x-x_j)
$$
passing through the points $\big(j_k,f(x_j)\big)$ where $j=0,\ldots,n$.
By definition
$$
	c_i=f[x_0,x_1,\ldots,x_i].
$$
By the Theorem of Polynomial Interpolation Error from page 315 in
Kincaid and Cheney we have that
$$
	f(x)-p_n(x)={1\over (n+1)!} f^{(n+1)}(\xi)\prod_{i=0}^n (x-x_i).
$$
Since $f(x)$ is a polynomial of degree $k$ then $f^{(n+1)}(x)=0$.
Therefore $f(x)=p_n(x)$.  This implies that $p_n(x)$ is
a polynomial of degree $k$.  Thus $c_n=f[x_0,x_1,\ldots,x_n]=0$.
\vfill\eject

\qn [Kincaid and Cheney Problem 6.2\#25] \ Establish an iteration
method for solving $f(x)=0$ as follows:  Let $q_n$ be the
quadratic interpolating polynomial through the points
$\big(x_n,f(x_n)\big)$, $\big(x_{n-1},f(x_{n-1})\big)$ and
$\big(x_{n-2},f(x_{n-2})\big)$ and define $x_{n+1}$ to be
the zero of $q_n$ closest to $x_n$.
\medskip\noindent
As an aside this method was established by D.E.~M\"uller in 1956
and is discussed in detail in {\it A Survey of Numerical Mathematics}
by D.~Young and R.~Gregory.

Newton's divided difference formula yields
$$
	q_n(x)=f[x_n]+f[x_n,x_{n-1}](x-x_n)
		+f[x_n,x_{n-1},x_{n-2}](x-x_n)(x-x_{n-1})$$
where the $f[x]$, $f[x,y]$ and $f[x,y,z]$ are given recursively as
$$\eqalign{
	f[x]&=f(x)\cr
	f[x,y]&=(f[y]-f[x])/(y-x)\cr
	f[x,y,z]&=(f[y,z]-f[x,y])/(z-x).\cr
}$$
Writing $c_0=f[x_n]$, $c_1=f[x_n,x_{n-1}]$ and $c_2=f[x_n,x_{n-1},x_{n-2}]$
we see that
$$\eqalign{
	q_n(x)&=c_0+\big(c_1+c_2(x_n-x_{n-1})\big)(x-x_n)+c_2(x-x_n)^2\cr
		&=a(x-x_n)^2+b(x-x_n)+c}
$$
where $a=c_2$, $b=c_1+c_2(x_n-x_{n-1})$ and $c=c_0.$
Now use the quadratic formula to solve for $x_{n+1}$.
First define
$$
	M_1=-b+\sqrt{b^2-4ac}
\qquad\hbox{and}\qquad
	M_2=-b-\sqrt{b^2-4ac}.
$$
Then
$$
	x_{n+1}=x_n+\cases{2c/ M_2 & if $|M_1|<|M_2|$ \cr
		2c/M_1 & otherwise. }
$$
Note that we have written the quadratic formula as
$$
	{2c\over -b\mp\sqrt{b^2-4ac}}
$$
with the square root
in the denominator to minimize loss of precision.  For example, if
$|M_1|$ is small, this means the terms $-b$ and $\sqrt{b^2-4ac}$
nearly cancelled.  This cancellation causes the value of $M_1$ to
lose precision.  Therefore, we compute $x_{n+1}$ using $M_2$.  Similarly,
if $M_2$ is small, we compute $x_{n+1}$ using $M_1$.
\vfill\eject
	
\qn [Kincaid and Cheney Problem 6.2\#1] \ For $n=5, 10$ and $15$ find
the Newton interpolating polynomial $p_n$ for the function
$f(x)=1/(1+x^2)$ on the interval $[-5,5]$.  Use equally spaced nodes.
In each case, compute $f(x)-p_n(x)$ for $30$ equally spaced
points in $[-5,5]$ in order to see the divergence of $p_n$ from $f$.
\medskip\noindent
The program is
\listing{prog9.c}
and my plotting script is
\listing{plot9.plt}
\vfill\eject

\noindent Output for $n=5$
\bigskip
\input p5.tex
\par
\listout{p5.dat}
\vfill\eject

\noindent Output for $n=10$
\bigskip
\input p10.tex
\par
\listout{p10.dat}
\vfill\eject

\noindent Output for $n=15$
\bigskip
\input p15.tex
\par
\listout{p15.dat}
\vfill\eject

\qn [Kincaid and Cheney Problem 6.2\#2] \ Program M\"uller's
method for finding roots as described in Problem 6.2.25 above and test it.
\medskip\noindent
The program is
\listing{prog10.c}
The program can be compiled to solve roots for different functions as
\medskip
\beginverbatimb
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
\endverbatim
\vfill\eject
\noindent
The output is
\listout{prog10a.out}
and
\listout{prog10b.out}
and
\listout{prog10c.out}
M\"uller's 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.

\end
