Problem 1
edit
Consider the system
A
x
=
b
{\displaystyle Ax=b\,\!}
. The GMRES method starts with a point
x
0
{\displaystyle \,\!x_{0}}
and normalizes the residual
r
0
=
b
−
A
x
0
{\displaystyle \,\!r_{0}=b-Ax_{0}}
so that
v
1
=
r
0
ν
{\displaystyle \,\!v_{1}={\frac {r_{0}}{\nu }}}
has 2-norm one. It then constructs orthonormal Krylov bases
V
k
=
(
v
1
v
2
⋯
v
m
)
{\displaystyle \,\!V_{k}=(v_{1}\,v_{2}\,\cdots v_{m})}
satisfying
A
V
k
=
V
k
+
1
H
k
{\displaystyle \,\!AV_{k}=V_{k+1}H_{k}}
where
H
k
{\displaystyle \,\!H_{k}}
is a
(
k
+
1
)
×
k
{\displaystyle \,\!(k+1)\times k}
upper Hessenberg matrix. One then looks for an approximation to
x
{\displaystyle \displaystyle x}
of the form
x
(
c
)
=
x
0
+
V
k
c
{\displaystyle \,\!x(c)=x_{0}+V_{k}c}
choosing
c
k
{\displaystyle \,\!c_{k}}
so that
‖
r
(
c
)
‖
=
‖
b
−
A
x
(
c
)
‖
{\displaystyle \,\!\|r(c)\|=\|b-Ax(c)\|}
is minimized, where
‖
⋅
‖
{\displaystyle \,\!\|\cdot \|}
is the usual Euclidean norm.
Problem 1a
edit
Show that
c
k
{\displaystyle \,\!c_{k}}
minimizes
‖
ν
e
1
−
H
k
c
‖
{\displaystyle \|\nu e_{1}-H_{k}c\|}
.
Solution 1a
edit
We wish to show that
‖
b
−
A
x
(
c
)
‖
=
‖
ν
e
1
−
H
k
c
‖
{\displaystyle \,\!\|b-Ax(c)\|=\|\nu e_{1}-H_{k}c\|}
‖
b
−
A
x
(
c
)
‖
=
‖
b
−
A
(
x
0
+
V
k
c
)
‖
=
‖
b
−
A
x
0
−
A
V
k
c
‖
=
‖
r
0
−
A
V
k
c
‖
=
‖
r
0
−
V
k
+
1
H
k
c
‖
=
‖
ν
v
1
−
V
k
+
1
H
k
c
‖
=
‖
V
k
+
1
(
ν
e
1
−
H
k
c
)
⏟
h
c
‖
=
(
V
k
+
1
h
c
,
V
k
+
1
h
c
)
1
2
=
(
(
V
k
+
1
h
c
)
T
V
k
+
1
h
c
)
1
2
=
(
h
c
T
V
k
+
1
T
V
k
+
1
h
c
)
1
2
=
(
h
c
T
h
c
)
1
2
=
‖
h
c
‖
=
‖
ν
e
1
−
H
k
c
‖
{\displaystyle {\begin{aligned}\|b-Ax(c)\|&=\|b-A(x_{0}+V_{k}c)\|\\&=\|b-Ax_{0}-AV_{k}c\|\\&=\|r_{0}-AV_{k}c\|\\&=\|r_{0}-V_{k+1}H_{k}c\|\\&=\|\nu v_{1}-V_{k+1}H_{k}c\|\\&=\|V_{k+1}\underbrace {(\nu e_{1}-H_{k}c)} _{h_{c}}\|\\&=(V_{k+1}h_{c},V_{k+1}h_{c})^{\frac {1}{2}}\\&=((V_{k+1}h_{c})^{T}V_{k+1}h_{c})^{\frac {1}{2}}\\&=(h_{c}^{T}V_{k+1}^{T}V_{k+1}h_{c})^{\frac {1}{2}}\\&=(h_{c}^{T}h_{c})^{\frac {1}{2}}\\&=\|h_{c}\|\\&=\|\nu e_{1}-H_{k}c\|\end{aligned}}}
Problem 1b
edit
Suppose we choose to solve the least squares problem in part a for
c
k
{\displaystyle c_{k}\!\,}
by the method of orthogonal triangularization (QR). What is the order of the floating point operation for this method. Give reasons.
Solution 1b
edit
We would like to transform
H
k
{\displaystyle H_{k}\!\,}
, a
(
k
+
1
)
×
k
{\displaystyle (k+1)\times k\!\,}
upper Hessenberg matrix, into QR form.
The cost is on the order of
4
k
2
+
1
2
k
2
=
4.5
k
2
{\displaystyle 4k^{2}+{\frac {1}{2}}k^{2}=4.5k^{2}}
from the cost of Given's Rotations and backsolving.
Given's Rotations Cost
edit
We need
k
{\displaystyle k\!\,}
Given's Rotation multiplies to zero out each of the
k
{\displaystyle k\!\,}
subdiagonal entries and hence transform
H
k
{\displaystyle H_{k}\!\,}
into upper triangular form
R
{\displaystyle R\!\,}
,
Each successive Given's Rotations multiply on an upper Hessenberg matrix requires four fewer multiplies because each previous subdiagonal entry has been zeroed out by a Given's Rotation multiply.
Hence the cost of
k
{\displaystyle k\!\,}
Given's Rotations multiplies is
4
k
+
4
(
k
−
1
)
+
…
+
4
⋅
1
<
4
k
⋅
k
=
4
k
2
{\displaystyle 4k+4(k-1)+\ldots +4\cdot 1<4k\cdot k=4k^{2}}
Back Solving Cost
edit
R
{\displaystyle R\!\,}
is a
(
k
+
1
)
×
k
{\displaystyle (k+1)\times k\!\,}
upper triangular matrix with last row zero. Hence, we need to backsolve
k
{\displaystyle k\!\,}
upper triangular rows.
1
+
2
+
…
+
k
=
k
(
k
+
1
)
2
=
k
2
2
+
k
2
{\displaystyle 1+2+\ldots +k={\frac {k(k+1)}{2}}={\frac {k^{2}}{2}}+{\frac {k}{2}}}
Problem 2
edit
We wish to approximate the integral
I
(
f
)
≡
∫
a
b
f
(
x
)
d
x
{\displaystyle I(f)\equiv \int _{a}^{b}f(x)dx}
Problem 2a
edit
Solution 2a
edit
The composite trapezoidal rule is
Q
T
,
n
(
f
)
=
h
2
(
f
(
a
)
+
f
(
b
)
)
+
h
∑
k
=
1
n
−
1
f
(
x
k
)
{\displaystyle Q_{T,n}(f)={\frac {h}{2}}(f(a)+f(b))+h\sum _{k=1}^{n-1}f(x_{k})\!\,}
The error is
I
(
f
)
−
Q
T
,
n
(
f
)
=
−
(
b
−
a
)
f
′
′
(
ξ
)
12
h
2
{\displaystyle I(f)-Q_{T,n}(f)=-{\frac {(b-a)f^{\prime \prime }(\xi )}{12}}h^{2}}
where
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]\!\,}
.
Derivation of Composite Trapezoidal Error
edit
The local error, the error on one interval, is
−
h
3
12
f
′
′
(
η
)
{\displaystyle -{\frac {h^{3}}{12}}f^{\prime \prime }(\eta )\!\,}
.
Observe that
n
min
η
∈
[
a
,
b
]
f
′
′
(
η
i
)
≤
∑
i
=
1
n
f
′
′
(
η
i
)
≤
n
max
η
∈
[
a
,
b
]
f
′
′
(
η
i
)
{\displaystyle n\min _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\leq \sum _{i=1}^{n}f^{\prime \prime }(\eta _{i})\leq n\max _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\!\,}
which implies
min
η
∈
[
a
,
b
]
f
′
′
(
η
i
)
≤
∑
i
=
1
n
f
′
′
(
η
i
)
n
≤
max
η
∈
[
a
,
b
]
f
′
′
(
η
i
)
{\displaystyle \min _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\leq \sum _{i=1}^{n}{\frac {f^{\prime \prime }(\eta _{i})}{n}}\leq \max _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\!\,}
Hence, the Intermediate Value Theorem implies there exists a
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]\!\,}
such that
f
′
′
(
ξ
)
=
∑
i
=
1
n
f
′
′
(
η
i
)
n
{\displaystyle f^{\prime \prime }(\xi )=\sum _{i=1}^{n}{\frac {f^{\prime \prime }(\eta _{i})}{n}}\!\,}
.
Multiplying both sides of the equation by
n
{\displaystyle n\!\,}
,
n
f
′
′
(
ξ
)
=
∑
i
=
1
n
f
′
′
(
η
i
)
{\displaystyle nf^{\prime \prime }(\xi )=\sum _{i=1}^{n}f^{\prime \prime }(\eta _{i})\!\,}
Using this relationship, we have
I
(
f
)
−
Q
T
,
n
(
f
)
=
∑
i
=
1
n
I
i
(
f
)
−
Q
i
(
f
)
=
∑
i
=
1
n
[
−
h
3
12
f
′
′
(
η
i
)
]
=
−
h
3
12
f
′
′
(
ξ
)
n
=
−
h
3
12
f
′
′
(
ξ
)
(
b
−
a
h
)
=
−
h
2
12
f
′
′
(
ξ
)
(
b
−
a
)
{\displaystyle {\begin{aligned}I(f)-Q_{T,n}(f)&=\sum _{i=1}^{n}I_{i}(f)-Q_{i}(f)\\&=\sum _{i=1}^{n}\left[-{\frac {h^{3}}{12}}f^{\prime \prime }(\eta _{i})\right]\\&=-{\frac {h^{3}}{12}}f^{\prime \prime }(\xi )n\\&=-{\frac {h^{3}}{12}}f^{\prime \prime }(\xi )({\frac {b-a}{h}})\\&=-{\frac {h^{2}}{12}}f^{\prime \prime }(\xi )(b-a)\end{aligned}}}
Derivation of Local Error
edit
The error in polynomial interpolation can be found by using the following theorem:
Assume
f
n
+
1
{\displaystyle f^{n+1}\!\,}
exists on
[
a
,
b
]
{\displaystyle [a,b]\!\,}
and
{
x
0
,
x
1
,
…
,
x
n
|
x
∈
[
a
,
b
]
}
.
{\displaystyle \{x_{0},x_{1},\ldots ,x_{n}|x\in [a,b]\}.\!\,}
P
n
{\displaystyle P_{n}\!\,}
interpolates
f
{\displaystyle f\!\,}
at
{
x
j
}
j
=
0
n
{\displaystyle \{x_{j}\}_{j=0}^{n}\!\,}
.
Then there is a
ξ
x
{\displaystyle \xi _{x}\!\,}
(
ξ
{\displaystyle \xi \!\,}
is dependent on
x
{\displaystyle x\!\,}
) such that
f
(
x
)
−
p
n
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
.
.
.
(
x
−
x
n
)
(
n
+
1
)
!
f
(
n
+
1
)
(
ξ
x
)
{\displaystyle f(x)-p_{n}(x)={\frac {(x-x_{0})(x-x_{1})...(x-x_{n})}{(n+1)!}}f^{(n+1)}(\xi _{x})\!\,}
where
ξ
x
{\displaystyle \xi _{x}\!\,}
lies in
(
m
,
M
)
{\displaystyle (m,M)\!\,}
,
m
=
min
{
x
0
,
x
1
,
…
,
x
n
,
X
}
{\displaystyle m=\min\{x_{0},x_{1},\ldots ,x_{n},X\}\!\,}
,
M
=
max
{
x
0
,
x
1
,
…
,
x
n
,
X
}
{\displaystyle M=\max\{x_{0},x_{1},\ldots ,x_{n},X\}\!\,}
Applying the theorem yields,
f
(
x
)
−
p
1
(
x
)
=
(
x
−
a
)
(
x
−
b
)
f
′
′
(
ξ
x
)
2
{\displaystyle f(x)-p_{1}(x)=(x-a)(x-b){\frac {f^{\prime \prime }(\xi _{x})}{2}}}
Hence,
E
(
f
)
=
∫
a
b
f
(
x
)
−
p
1
(
x
)
d
x
=
∫
a
b
(
x
−
a
)
⏟
>
0
(
x
−
b
)
⏟
<
0
⏟
w
(
x
)
f
′
′
(
ξ
x
)
2
d
x
{\displaystyle {\begin{aligned}E(f)&=\int _{a}^{b}f(x)-p_{1}(x)dx\\&=\int _{a}^{b}\underbrace {\underbrace {(x-a)} _{>0}\underbrace {(x-b)} _{<0}} _{w(x)}{\frac {f^{\prime \prime }(\xi _{x})}{2}}dx\end{aligned}}}
Since
a
{\displaystyle a\!\,}
is the start of the interval,
x
−
a
{\displaystyle x-a\!\,}
is always positive. Conversely, since
b
{\displaystyle b\!\,}
is the end of the interval,
x
−
b
{\displaystyle x-b\!\,}
is always negative. Hence,
w
(
x
)
<
0
{\displaystyle w(x)<0\!\,}
is always of one sign. Hence from the mean value theorem of integration , there exists a
ζ
∈
[
a
,
b
]
{\displaystyle \zeta \in [a,b]\!\,}
such that
E
(
f
)
=
f
′
′
(
ζ
)
2
∫
a
b
(
x
−
a
)
(
x
−
b
)
d
x
{\displaystyle E(f)={\frac {f^{\prime \prime }(\zeta )}{2}}\int _{a}^{b}(x-a)(x-b)dx\!\,}
Note that
f
′
′
(
ζ
)
{\displaystyle f^{\prime \prime }(\zeta )}
is a constant and does not depend on
x
{\displaystyle x\!\,}
.
Integrating
∫
a
b
(
x
−
a
)
(
x
−
b
)
d
x
{\displaystyle \int _{a}^{b}(x-a)(x-b)dx\!\,}
, yields,
E
(
f
)
=
f
′
′
(
ζ
)
2
(
−
(
b
−
a
)
6
)
=
−
f
′
′
(
ζ
)
(
b
−
a
)
12
{\displaystyle {\begin{aligned}E(f)&={\frac {f^{\prime \prime }(\zeta )}{2}}\left(-{\frac {(b-a)}{6}}\right)\\&=-{\frac {f^{\prime \prime }(\zeta )(b-a)}{12}}\!\,\end{aligned}}}
Problem 2b
edit
Use the error formula to derive a new quadrature rule obtained by performing one step of extrapolation on the composite trapezoidal rule. What is this rule, and how does its error depend on
h
{\displaystyle h\!\,}
? You may assume here that
f
{\displaystyle f\!\,}
is as smooth as you need it to be.
Solution 2b
edit
STOER AND BUESCH pg 162
INTRODUCTION TO APPLIED NUMERICAL ANALYSIS by RICHARD HAMMING pg 178
The error for the composite trapezoidal rule at
n
{\displaystyle n\!\,}
points in
[
a
,
b
]
{\displaystyle [a,b]\!\,}
is
E
n
=
I
(
f
)
−
Q
T
,
n
=
−
(
b
−
a
)
h
2
12
+
O
(
h
3
)
{\displaystyle E_{n}=I(f)-Q_{T,n}={\frac {-(b-a)h^{2}}{12}}+{\mathcal {O}}(h^{3})\!\,}
With
2
n
{\displaystyle 2n\!\,}
points, there are twice as many intervals, but the intervals are half as wide. Hence, the error for the composite trapezoidal rule at
2
n
{\displaystyle 2n\!\,}
points in
[
a
,
b
]
{\displaystyle [a,b]\!\,}
is
E
2
n
=
I
(
f
)
−
Q
T
,
2
n
=
2
⋅
−
(
b
−
a
)
(
h
2
)
2
12
+
O
(
h
3
)
=
−
(
b
−
a
)
h
2
24
+
O
(
h
3
)
{\displaystyle E_{2n}=I(f)-Q_{T,2n}=2\cdot {\frac {-(b-a)({\frac {h}{2}})^{2}}{12}}+{\mathcal {O}}(h^{3})={\frac {-(b-a)h^{2}}{24}}+{\mathcal {O}}(h^{3})\!\,}
We can eliminate the
h
2
{\displaystyle h^{2}\!\,}
term by choosing using an appropriate linear combination of
E
n
{\displaystyle E_{n}\!\,}
and
E
2
n
{\displaystyle E_{2n}\!\,}
. This gives a new error rule with
h
3
{\displaystyle h^{3}\!\,}
error.
E
n
−
2
E
2
n
=
−
(
b
−
a
)
h
2
12
−
2
⋅
−
(
b
−
a
)
h
2
24
+
O
(
h
3
)
=
O
(
h
3
)
{\displaystyle {\begin{aligned}E_{n}-2E_{2n}&={\frac {-(b-a)h^{2}}{12}}-2\cdot {\frac {-(b-a)h^{2}}{24}}+{\mathcal {O}}(h^{3})\\&={\mathcal {O}}(h^{3})\end{aligned}}}
Substituting our equations for
E
n
{\displaystyle E_{n}\!\,}
and
E
2
n
{\displaystyle E_{2n}\!\,}
on the left had side gives
I
(
f
)
−
Q
T
,
n
−
2
I
(
f
)
−
2
Q
T
,
2
n
=
O
(
h
3
)
{\displaystyle I(f)-Q_{T,n}-2I(f)-2Q_{T,2n}={\mathcal {O}}(h^{3})\!\,}
I
(
f
)
=
2
Q
T
,
2
n
−
Q
T
,
n
+
O
(
h
3
)
{\displaystyle I(f)=2Q_{T,2n}-Q_{T,n}+{\mathcal {O}}(h^{3})\!\,}
If we call our new rule
Q
^
{\displaystyle {\hat {Q}}\!\,}
we have
Q
^
=
2
Q
T
,
2
n
−
Q
T
,
n
{\displaystyle {\hat {Q}}=2Q_{T,2n}-Q_{T,n}\!\,}
whose error is on the order of
O
(
h
3
)
{\displaystyle {\mathcal {O}}(h^{3})\!\,}
Problem 3
edit
Consider the shifted QR iteration for computing a 2 x 2 matrix
A
{\displaystyle A\!\,}
: starting with
A
0
=
A
{\displaystyle A_{0}=A\!\,}
, compute
A
k
−
σ
k
I
=
Q
k
R
k
A
k
+
1
=
R
k
Q
k
+
σ
k
I
{\displaystyle A_{k}-\sigma _{k}I=Q_{k}R_{k}\quad \quad \quad A_{k+1}=R_{k}Q_{k}+\sigma _{k}I}
where
σ
k
{\displaystyle \sigma _{k}\!\,}
is a scalar shift.
Problem 3a
edit
If
A
k
=
[
a
11
a
12
a
21
a
22
]
,
{\displaystyle A_{k}={\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}},}
specify the orthogonal matrix
Q
k
{\displaystyle Q_{k}\!\,}
used to perform this step when Givens rotations are used. The matrix should be described in terms of the entries of the shifted matrix.
Solution 3a
edit
Let
A
~
k
{\displaystyle {\tilde {A}}_{k}\!\,}
be the
σ
k
{\displaystyle \sigma _{k}\!\,}
-shifted matrix of
A
k
{\displaystyle A_{k}\!\,}
i.e.
A
~
k
=
A
k
−
σ
k
I
=
[
a
11
−
σ
k
a
12
a
21
a
22
−
σ
k
]
,
{\displaystyle {\tilde {A}}_{k}=A_{k}-\sigma _{k}I={\begin{bmatrix}a_{11}-\sigma _{k}&a_{12}\\a_{21}&a_{22}-\sigma _{k}\end{bmatrix}},}
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
is an orthogonal 2x2 Given's rotations matrix.
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
's entries are chosen such that when
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
is pre-multiplied against the 2x2 matrix
A
~
k
{\displaystyle {\tilde {A}}_{k}\!\,}
,
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
will zero out the (2,1) entry of
A
~
k
{\displaystyle {\tilde {A}}_{k}\!\,}
and scale
A
~
k
{\displaystyle {\tilde {A}}_{k}\!\,}
's remaining three entries i.e.
Q
k
T
A
~
k
=
[
∗
∗
0
∗
]
=
R
k
{\displaystyle Q_{k}^{T}{\tilde {A}}_{k}={\begin{bmatrix}*&*\\0&*\end{bmatrix}}=R_{k}}
where
∗
{\displaystyle *\!\,}
denotes calculable scalar values we are not interested in and
R
k
{\displaystyle R_{k}\!\,}
is our desired upper triangular matrix sought by the QR algorithm.
Since
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
is orthogonal, the above equation implies
Q
k
T
A
~
k
=
R
k
Q
k
Q
k
T
A
~
k
=
Q
k
R
k
A
~
k
=
Q
k
R
k
{\displaystyle {\begin{aligned}Q_{k}^{T}{\tilde {A}}_{k}&=R_{k}\\Q_{k}Q_{k}^{T}{\tilde {A}}_{k}&=Q_{k}R_{k}\\{\tilde {A}}_{k}&=Q_{k}R_{k}\end{aligned}}}
The Given's rotation
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
is given by
Q
k
T
=
[
c
s
−
s
c
]
{\displaystyle Q_{k}^{T}={\begin{bmatrix}c&s\\-s&c\end{bmatrix}}}
where
c
=
a
11
−
σ
k
r
s
=
a
21
r
r
=
(
a
11
−
σ
k
)
2
+
a
21
2
{\displaystyle {\begin{aligned}c&={\frac {a_{11}-\sigma _{k}}{r}}\\s&={\frac {a_{21}}{r}}\\r&={\sqrt {(a_{11}-\sigma _{k})^{2}+a_{21}^{2}}}\end{aligned}}}
Taking the transpose of
Q
k
T
{\displaystyle Q_{k}^{T}\!\,}
yields
Q
k
=
[
c
−
s
s
c
]
{\displaystyle Q_{k}={\begin{bmatrix}c&-s\\s&c\end{bmatrix}}}
Problem 3b
edit
Solution 3b
edit
From hypothesis and part (a),
A
k
+
1
=
R
k
Q
k
+
σ
k
I
=
[
∗
∗
0
r
22
]
[
c
−
s
s
c
]
+
[
σ
k
0
0
σ
k
]
{\displaystyle {\begin{aligned}A_{k+1}&=R_{k}Q_{k}+\sigma _{k}I\\&={\begin{bmatrix}*&*\\0&r_{22}\end{bmatrix}}{\begin{bmatrix}c&-s\\s&c\end{bmatrix}}+{\begin{bmatrix}\sigma _{k}&0\\0&\sigma _{k}\end{bmatrix}}\end{aligned}}}
Let
A
k
+
1
(
2
,
1
)
{\displaystyle A_{k+1}(2,1)\!\,}
be the (2,1) entry of
A
k
{\displaystyle A_{k}\!\,}
. Using the above equation, we can find
A
k
+
1
(
2
,
1
)
{\displaystyle A_{k+1}(2,1)\!\,}
by finding the inner product of the second row of
R
k
{\displaystyle R_{k}\!\,}
and first column
Q
k
{\displaystyle Q_{k}\!\,}
and adding the (2,1) entry of
σ
k
I
{\displaystyle \sigma _{k}I\!\,}
i.e.
A
k
+
1
(
2
,
1
)
=
(
0
⋅
c
+
r
22
⋅
s
)
+
0
=
r
22
s
{\displaystyle {\begin{aligned}A_{k+1}(2,1)&=(0\cdot c+r_{22}\cdot s)+0\\&=r_{22}s\end{aligned}}}
We need to find the value of
r
22
{\displaystyle r_{22}\!\,}
so we need to calculate
R
k
{\displaystyle R_{k}\!\,}
.
From hypothesis and the orthogonality of
Q
k
{\displaystyle Q_{k}\!\,}
, we have
A
k
−
σ
k
I
=
Q
k
R
k
Q
k
T
(
A
k
−
σ
k
I
)
=
Q
k
T
Q
k
R
k
Q
k
T
(
A
k
−
σ
k
I
)
=
R
k
R
k
=
Q
k
T
(
A
k
−
σ
k
I
)
=
[
c
s
−
s
c
]
[
a
11
−
σ
k
a
12
δ
a
22
−
σ
k
]
{\displaystyle {\begin{aligned}A_{k}-\sigma _{k}I&=Q_{k}R_{k}\\Q_{k}^{T}(A_{k}-\sigma _{k}I)&=Q_{k}^{T}Q_{k}R_{k}\\Q_{k}^{T}(A_{k}-\sigma _{k}I)&=R_{k}\\R_{k}&=Q_{k}^{T}(A_{k}-\sigma _{k}I)\\&={\begin{bmatrix}c&s\\-s&c\end{bmatrix}}{\begin{bmatrix}a_{11}-\sigma _{k}&a_{12}\\\delta &a_{22}-\sigma _{k}\end{bmatrix}}\end{aligned}}}
From
R
k
{\displaystyle R_{k}\!\,}
, we can find its (2,2) entry
r
2
2
{\displaystyle r_{2}2\!\,}
by using inner products.
r
22
=
−
s
⋅
a
12
+
c
⋅
(
a
22
−
σ
k
)
{\displaystyle r_{22}=-s\cdot a_{12}+c\cdot (a_{22}-\sigma _{k})\!\,}
Now that we have
r
22
{\displaystyle r_{22}\!\,}
we can calculate
A
k
+
1
(
2
,
1
)
{\displaystyle A_{k+1}(2,1)\!\,}
by using appropriate substitutions
A
k
+
1
(
2
,
1
)
=
r
22
⋅
s
=
(
−
s
a
12
+
c
(
a
22
−
σ
k
)
)
s
=
−
s
2
a
12
+
c
s
(
a
22
−
σ
k
)
=
−
δ
2
a
12
(
a
11
−
σ
k
)
2
+
δ
2
+
(
a
11
−
σ
k
)
(
a
22
−
σ
k
)
δ
(
a
11
−
σ
k
)
2
+
δ
2
=
−
δ
2
a
12
+
δ
(
a
11
−
σ
k
)
(
a
22
−
σ
k
)
(
a
11
−
σ
k
)
2
+
δ
2
≈
−
δ
2
a
12
+
δ
(
a
11
−
σ
k
)
(
a
22
−
σ
k
)
(
a
11
−
σ
k
)
2
{\displaystyle {\begin{aligned}A_{k+1}(2,1)&=r_{22}\cdot s\\&=(-sa_{12}+c(a_{22}-\sigma _{k}))s\\&=-s^{2}a_{12}+cs(a_{22}-\sigma _{k})\\&={\frac {-\delta ^{2}a_{12}}{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}+{\frac {(a_{11}-\sigma _{k})(a_{22}-\sigma _{k})\delta }{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}\\&={\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}\\&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}}}\end{aligned}}}
since
δ
{\displaystyle \delta \!\,}
is a small number.
Let our shift
σ
k
=
a
22
{\displaystyle \sigma _{k}=a_{22}\!\,}
. Then the above equation yields,
A
k
+
1
(
2
,
1
)
≈
−
δ
2
a
12
+
δ
(
a
11
−
σ
k
)
(
a
22
−
σ
k
)
(
a
11
−
σ
k
)
2
≈
−
δ
2
a
12
+
δ
(
a
11
−
a
22
)
(
a
22
−
a
22
)
(
a
11
−
a
22
)
2
=
−
δ
2
a
12
(
a
11
−
a
22
)
2
{\displaystyle {\begin{aligned}A_{k+1}(2,1)&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}}}\\&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-a_{22})(a_{22}-a_{22})}{(a_{11}-a_{22})^{2}}}\\&={\frac {-\delta ^{2}a_{12}}{(a_{11}-a_{22})^{2}}}\end{aligned}}}
Hence,
|
A
k
+
1
(
2
,
1
)
|
≈
|
δ
2
a
12
(
a
11
−
a
22
)
2
|
≤
|
δ
2
|
{\displaystyle {\begin{aligned}|A_{k+1}(2,1)|&\approx \left|{\frac {\delta ^{2}a_{12}}{(a_{11}-a_{22})^{2}}}\right|\\&\leq |\delta ^{2}|\end{aligned}}}
since
|
a
12
|
≤
(
a
11
−
a
22
)
2
{\displaystyle |a_{12}|\leq (a_{11}-a_{22})^{2}\!\,}
Hence we have shown that
A
k
+
1
(
2
,
1
)
{\displaystyle A_{k+1}(2,1)\!\,}
is of order
δ
2
{\displaystyle \delta ^{2}\!\,}
.
If
δ
{\displaystyle \delta \!\,}
is small, the QR convergence rate will be quadratic.