Derive the one-, two-, and three-point Gaussian quadrature formulas such that
∫
−
1
1
f
(
x
)
x
4
d
x
≈
∑
j
=
1
n
f
(
x
j
)
w
j
{\displaystyle \int _{-1}^{1}f(x)x^{4}dx\approx \sum _{j=1}^{n}f(x_{j})w_{j}\!\,}
Give Bounds on the error of these formulas.
Find Orthogonal Polynomials w.r.t. Weighting Function
edit
For this problem, we need the first three orthogonal polynomials
P
1
,
P
2
,
and
P
3
{\displaystyle P_{1},P_{2},{\mbox{ and }}P_{3}\!\,}
with respect to a weighted inner product
⟨
f
,
g
⟩
=
∫
−
1
1
f
(
x
)
g
(
x
)
w
(
x
)
d
x
{\displaystyle \left\langle f,g\right\rangle =\int _{-1}^{1}f(x)g(x)w(x)dx\!\,}
where
w
(
x
)
=
x
4
{\displaystyle w(x)=x^{4}\!\,}
in our case. This can be done using Gram-Schmidt Orthogonalization on the basis
{
1
,
x
,
x
2
,
x
3
}
{\displaystyle \left\{1,x,x^{2},x^{3}\right\}\!\,}
Zero Order Polynomial
edit
Let
P
0
(
x
)
=
1
{\displaystyle P_{0}(x)=1\!\,}
First Order Polynomial
edit
P
1
(
x
)
=
x
−
⟨
1
,
x
⟩
⟨
1
,
1
⟩
⋅
1
=
x
−
∫
−
1
1
x
5
d
x
∫
−
1
1
x
4
d
x
⏟
0
2
/
5
=
x
{\displaystyle {\begin{aligned}P_{1}(x)&=x-{\frac {\left\langle 1,x\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1\\&=x-\underbrace {\frac {\int _{-1}^{1}x^{5}dx}{\int _{-1}^{1}x^{4}dx}} _{\frac {0}{2/5}}\\&=x\end{aligned}}}
Second Order Polynomial
edit
P
2
(
x
)
=
x
2
−
⟨
1
,
x
2
⟩
⟨
1
,
1
⟩
⋅
1
−
⟨
x
,
x
2
⟩
⟨
x
,
x
⟩
⋅
x
=
x
2
−
∫
−
1
1
x
6
d
x
∫
−
1
1
x
4
d
x
−
x
⋅
∫
−
1
1
x
7
d
x
∫
−
1
1
x
6
d
x
⏟
0
2
/
7
=
x
2
−
5
7
{\displaystyle {\begin{aligned}P_{2}(x)&=x^{2}-{\frac {\left\langle 1,x^{2}\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1-{\frac {\left\langle x,x^{2}\right\rangle }{\left\langle x,x\right\rangle }}\cdot x\\&=x^{2}-{\frac {\int _{-1}^{1}x^{6}dx}{\int _{-1}^{1}x^{4}dx}}-x\cdot \underbrace {\frac {\int _{-1}^{1}x^{7}dx}{\int _{-1}^{1}x^{6}dx}} _{\frac {0}{2/7}}\\&=x^{2}-{\frac {5}{7}}\end{aligned}}}
Third Order Polynomial
edit
P
3
(
x
)
=
x
3
−
⟨
1
,
x
3
⟩
⟨
1
,
1
⟩
⋅
1
−
⟨
x
,
x
3
⟩
⟨
x
,
x
⟩
⋅
x
−
⟨
x
2
−
5
7
,
x
3
⟩
⟨
x
2
−
5
7
,
x
2
−
5
7
⟩
⋅
(
x
2
−
5
7
)
=
x
3
−
0
−
7
2
x
⋅
∫
−
1
1
x
8
d
x
∫
−
1
1
x
6
d
x
−
0
=
x
3
−
7
9
x
{\displaystyle {\begin{aligned}P_{3}(x)&=x^{3}-{\frac {\left\langle 1,x^{3}\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1-{\frac {\left\langle x,x^{3}\right\rangle }{\left\langle x,x\right\rangle }}\cdot x-{\frac {\left\langle x^{2}-{\frac {5}{7}},x^{3}\right\rangle }{\left\langle x^{2}-{\frac {5}{7}},x^{2}-{\frac {5}{7}}\right\rangle }}\cdot (x^{2}-{\frac {5}{7}})\\&=x^{3}-0-{\frac {7}{2}}x\cdot {\frac {\int _{-1}^{1}x^{8}dx}{\int _{-1}^{1}x^{6}dx}}-0\\&=x^{3}-{\frac {7}{9}}x\end{aligned}}}
Find Roots of Polynomials
edit
The roots of these polynomials will be the interpolation nodes used in Gauss Quadrature.
Polynomial
Roots
P
1
(
x
)
=
x
0
P
2
(
x
)
=
x
2
−
5
7
±
5
7
P
3
(
x
)
=
x
3
−
7
9
x
0
,
±
7
9
{\displaystyle {\begin{array}{|c|c|}{\mbox{Polynomial}}&{\mbox{Roots}}\\\hline P_{1}(x)=x&0\\\\P_{2}(x)=x^{2}-{\frac {5}{7}}&\pm {\sqrt {\frac {5}{7}}}\\\\P_{3}(x)=x^{3}-{\frac {7}{9}}x&0,\pm {\sqrt {\frac {7}{9}}}\\\\\hline \end{array}}}
In Gaussian quadrature we have
w
j
=
∫
−
1
1
l
i
(
x
)
x
4
d
x
{\displaystyle w_{j}=\int _{-1}^{1}l_{i}(x)x^{4}dx\!\,}
, where
l
i
(
x
)
=
∏
j
=
1
,
j
≠
i
n
x
−
x
j
x
i
−
x
j
{\displaystyle l_{i}(x)=\prod _{j=1,j\neq i}^{n}{\frac {x-x_{j}}{x_{i}-x_{j}}}\!\,}
l
1
(
x
)
=
1
{\displaystyle l_{1}(x)=1\!\,}
w
1
=
∫
−
1
1
x
4
d
x
=
2
5
{\displaystyle w_{1}=\int _{-1}^{1}x^{4}dx={\frac {2}{5}}\!\,}
∫
−
1
1
f
(
x
)
x
4
d
x
≈
2
5
f
(
x
1
)
{\displaystyle \int _{-1}^{1}f(x)x^{4}dx\approx {\frac {2}{5}}f(x_{1})\!\,}
where
x
1
{\displaystyle x_{1}\!\,}
is the root of
P
1
(
x
)
{\displaystyle P_{1}(x)\!\,}
.
l
1
(
x
)
=
x
−
x
1
x
1
−
x
2
{\displaystyle l_{1}(x)={\frac {x-x_{1}}{x_{1}-x_{2}}}\!\,}
l
2
(
x
)
=
x
−
x
1
x
2
−
x
1
{\displaystyle l_{2}(x)={\frac {x-x_{1}}{x_{2}-x_{1}}}\!\,}
w
1
=
∫
−
1
1
x
−
x
1
x
1
−
x
2
x
4
d
x
=
2
x
2
5
(
x
2
−
x
1
)
{\displaystyle w_{1}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{1}-x_{2}}}x^{4}dx={\frac {2x_{2}}{5(x_{2}-x_{1})}}\!\,}
w
2
=
∫
−
1
1
x
−
x
1
x
2
−
x
1
x
4
d
x
=
2
x
1
5
(
x
1
−
x
2
)
{\displaystyle w_{2}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{2}-x_{1}}}x^{4}dx={\frac {2x_{1}}{5(x_{1}-x_{2})}}\!\,}
∫
−
1
1
f
(
x
)
x
4
d
x
≈
w
1
f
(
x
1
)
+
w
2
f
(
x
2
)
{\displaystyle \int _{-1}^{1}f(x)x^{4}dx\approx w_{1}f(x_{1})+w_{2}f(x_{2})\!\,}
where
x
1
and
x
2
{\displaystyle x_{1}{\mbox{ and }}x_{2}\!\,}
are the roots of
P
2
(
x
)
{\displaystyle P_{2}(x)\!\,}
.
l
1
(
x
)
=
x
−
x
2
x
1
−
x
2
x
−
x
3
x
1
−
x
3
{\displaystyle l_{1}(x)={\frac {x-x_{2}}{x_{1}-x_{2}}}{\frac {x-x_{3}}{x_{1}-x_{3}}}\!\,}
l
2
(
x
)
=
x
−
x
1
x
2
−
x
1
x
−
x
3
x
2
−
x
3
{\displaystyle l_{2}(x)={\frac {x-x_{1}}{x_{2}-x_{1}}}{\frac {x-x_{3}}{x_{2}-x_{3}}}\!\,}
l
3
(
x
)
=
x
−
x
1
x
3
−
x
1
x
−
x
2
x
3
−
x
2
{\displaystyle l_{3}(x)={\frac {x-x_{1}}{x_{3}-x_{1}}}{\frac {x-x_{2}}{x_{3}-x_{2}}}\!\,}
w
1
=
∫
−
1
1
x
−
x
2
x
1
−
x
2
x
−
x
3
x
1
−
x
3
x
4
d
x
=
1
(
x
1
−
x
2
)
(
x
1
−
x
3
)
(
2
7
+
2
x
2
x
3
5
)
{\displaystyle w_{1}=\int _{-1}^{1}{\frac {x-x_{2}}{x_{1}-x_{2}}}{\frac {x-x_{3}}{x_{1}-x_{3}}}x^{4}dx={\frac {1}{(x_{1}-x_{2})(x_{1}-x_{3})}}\left({\frac {2}{7}}+{\frac {2x_{2}x_{3}}{5}}\right)\!\,}
w
2
=
∫
−
1
1
x
−
x
1
x
2
−
x
1
x
−
x
3
x
2
−
x
3
x
4
d
x
=
1
(
x
2
−
x
1
)
(
x
2
−
x
3
)
(
2
7
+
2
x
1
x
3
5
)
{\displaystyle w_{2}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{2}-x_{1}}}{\frac {x-x_{3}}{x_{2}-x_{3}}}x^{4}dx={\frac {1}{(x_{2}-x_{1})(x_{2}-x_{3})}}\left({\frac {2}{7}}+{\frac {2x_{1}x_{3}}{5}}\right)\!\,}
w
3
=
∫
−
1
1
x
−
x
1
x
3
−
x
1
x
−
x
2
x
3
−
x
2
x
4
d
x
=
1
(
x
3
−
x
1
)
(
x
3
−
x
2
)
(
2
7
+
2
x
1
x
2
5
)
{\displaystyle w_{3}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{3}-x_{1}}}{\frac {x-x_{2}}{x_{3}-x_{2}}}x^{4}dx={\frac {1}{(x_{3}-x_{1})(x_{3}-x_{2})}}\left({\frac {2}{7}}+{\frac {2x_{1}x_{2}}{5}}\right)\!\,}
∫
−
1
1
f
(
x
)
x
4
d
x
≈
w
1
f
(
x
1
)
+
w
2
f
(
x
2
)
+
w
3
f
(
x
3
)
{\displaystyle \int _{-1}^{1}f(x)x^{4}dx\approx w_{1}f(x_{1})+w_{2}f(x_{2})+w_{3}f(x_{3})\!\,}
where
x
1
,
x
2
and
x
3
{\displaystyle x_{1}{\mbox{, }}x_{2}{\mbox{ and }}x_{3}\!\,}
are the roots of
P
3
(
x
)
{\displaystyle P_{3}(x)\!\,}
.
We know that this quadrature is exact for all polynomials of degree at most
2
n
−
1
{\displaystyle 2n-1\!\,}
.
We choose a polynomial
p
{\displaystyle p\!\,}
of degree at most
2
n
−
1
{\displaystyle 2n-1\!\,}
that Hermite interpolates i.e.
p
(
x
i
)
=
f
(
x
i
)
p
′
(
x
i
)
=
f
′
(
x
i
)
(
0
≤
i
≤
n
−
1
)
{\displaystyle p(x_{i})=f(x_{i})\quad p'(x_{i})=f'(x_{i})\quad (0\leq i\leq n-1)\!\,}
The error for this interpolation is
f
(
x
)
−
p
(
x
)
=
f
(
2
n
)
(
ζ
(
x
)
)
(
2
n
)
!
∏
i
=
0
n
−
1
(
x
−
x
i
)
2
{\displaystyle f(x)-p(x)={\frac {f^{(2n)}(\zeta (x))}{(2n)!}}\prod _{i=0}^{n-1}(x-x_{i})^{2}\!\,}
Compute the error of the quadrature as follows:
E
=
∫
−
1
1
f
(
x
)
x
4
d
x
−
∑
i
=
0
n
−
1
w
i
f
(
x
i
)
=
∫
−
1
1
f
(
x
)
x
4
d
x
−
∑
i
=
0
n
−
1
w
i
p
(
x
i
)
=
∫
−
1
1
f
(
x
)
x
4
d
x
−
∫
−
1
1
p
(
x
)
x
4
d
x
=
∫
−
1
1
(
f
(
x
)
−
p
(
x
)
)
x
4
d
x
=
∫
−
1
1
f
(
2
n
)
(
ζ
(
x
)
)
(
2
n
)
!
∏
i
=
0
n
−
1
(
x
−
x
i
)
2
⋅
x
4
d
x
=
f
(
2
n
)
(
ξ
)
(
2
n
)
!
∫
−
1
1
x
4
⋅
∏
i
=
0
n
−
1
(
x
−
x
i
)
2
⏟
(
p
n
(
x
)
)
2
d
x
{\displaystyle {\begin{aligned}E&=\int _{-1}^{1}f(x)x^{4}dx-\sum _{i=0}^{n-1}w_{i}f(x_{i})\\&=\int _{-1}^{1}f(x)x^{4}dx-\sum _{i=0}^{n-1}w_{i}p(x_{i})\\&=\int _{-1}^{1}f(x)x^{4}dx-\int _{-1}^{1}p(x)x^{4}dx\\&=\int _{-1}^{1}(f(x)-p(x))x^{4}dx\\&=\int _{-1}^{1}{\frac {f^{(2n)}(\zeta (x))}{(2n)!}}\prod _{i=0}^{n-1}(x-x_{i})^{2}\cdot x^{4}dx\\&={\frac {f^{(2n)}(\xi )}{(2n)!}}\int _{-1}^{1}x^{4}\cdot \underbrace {\prod _{i=0}^{n-1}(x-x_{i})^{2}} _{(p_{n}(x))^{2}}dx\end{aligned}}\!\,}
where the last line follows from the mean value theorem of integrals.
Notice that
p
n
(
x
)
{\displaystyle p_{n}(x)\!\,}
is simply the polynomials orthogonal with respect to weight function since
x
i
{\displaystyle x_{i}\!\,}
are the roots of the polynomials.
Hence, the error for 1 point gaussian quadrature is
E
=
f
(
2
)
(
ξ
)
(
2
)
!
∫
−
1
1
x
4
⋅
x
2
d
x
=
f
(
2
)
(
ξ
)
7
{\displaystyle {\begin{aligned}E&={\frac {f^{(2)}(\xi )}{(2)!}}\int _{-1}^{1}x^{4}\cdot x^{2}dx\\&={\frac {f^{(2)}(\xi )}{7}}\end{aligned}}\!\,}
The error for 2 point quadrature:
E
=
f
(
4
)
(
ξ
)
(
4
)
!
∫
−
1
1
x
4
⋅
(
x
2
−
5
7
)
2
d
x
{\displaystyle {\begin{aligned}E&={\frac {f^{(4)}(\xi )}{(4)!}}\int _{-1}^{1}x^{4}\cdot (x^{2}-{\frac {5}{7}})^{2}dx\end{aligned}}\!\,}
The error for 3 point quadrature:
E
=
f
(
6
)
(
ξ
)
(
6
)
!
∫
−
1
1
x
4
⋅
(
x
3
−
7
9
x
)
2
d
x
{\displaystyle {\begin{aligned}E&={\frac {f^{(6)}(\xi )}{(6)!}}\int _{-1}^{1}x^{4}\cdot (x^{3}-{\frac {7}{9}}x)^{2}dx\end{aligned}}\!\,}
We wish to solve
A
x
=
b
{\displaystyle Ax=b\!\,}
iteratively where
A
=
[
1
1
−
1
1
2
1
1
1
1
]
{\displaystyle A={\begin{bmatrix}1&1&-1\\1&2&1\\1&1&1\end{bmatrix}}\!\,}
Show that for this
A
{\displaystyle A\!\,}
the Jacobi method and the Gauss-Seidel method both converge. Explain why for this
A
{\displaystyle A\!\,}
one of these methods is better than the other.
Decompose
A
{\displaystyle A\!\,}
into its diagonal, lower, and upper parts i.e.
A
=
D
+
(
L
+
U
)
{\displaystyle A=D+(L+U)\!\,}
Derive Jacobi iteration by solving for x as follows:
A
x
=
b
(
D
+
(
L
+
U
)
)
x
=
b
D
x
+
(
L
+
U
)
x
=
b
D
x
=
b
−
(
L
+
U
)
x
x
=
D
−
1
(
b
−
(
L
+
U
)
x
)
{\displaystyle {\begin{aligned}Ax&=b\\(D+(L+U))x&=b\\Dx+(L+U)x&=b\\Dx&=b-(L+U)x\\x&=D^{-1}(b-(L+U)x)\end{aligned}}\!\,}
Decompose
A
{\displaystyle A\!\,}
into its diagonal, lower, and upper parts i.e.
A
=
(
D
+
L
)
+
U
{\displaystyle A=(D+L)+U\!\,}
Derive Gauss-Sediel iteration by solving for x as follows:
A
x
=
b
(
(
D
+
L
)
+
U
)
x
=
b
(
D
+
L
)
x
+
U
x
=
b
(
D
+
L
)
x
=
b
−
U
x
x
=
(
D
+
L
)
−
1
(
b
−
U
x
)
{\displaystyle {\begin{aligned}Ax&=b\\((D+L)+U)x&=b\\(D+L)x+Ux&=b\\(D+L)x&=b-Ux\\x&=(D+L)^{-1}(b-Ux)\end{aligned}}\!\,}
Convergence occurs for the Jacobi iteration if the spectral radius of
D
−
1
(
L
+
U
)
{\displaystyle D^{-1}(L+U)\!\,}
is less than 1 i.e.
ρ
(
D
−
1
(
L
+
U
)
)
<
1
{\displaystyle \rho (D^{-1}(L+U))<1\!\,}
The eigenvalues of the matrix
D
−
1
(
L
+
U
)
=
[
0
1
−
1
1
2
0
1
2
1
1
0
]
{\displaystyle D^{-1}(L+U)={\begin{bmatrix}0&1&-1\\{\frac {1}{2}}&0&{\frac {1}{2}}\\1&1&0\end{bmatrix}}\!\,}
are
λ
=
0
,
0
,
0
{\displaystyle \lambda =0,0,0\!\,}
i.e. the eigenvalue
0
{\displaystyle 0\!\,}
has order 3.
Therefore, the spectral radius is
ρ
=
0
<
1
{\displaystyle \rho =0<1\!\,}
.
Gauss-Seidel Convergence
edit
Convergence occurs for the Gauss-Seidel iteration if the spectral radius of
(
D
+
L
)
−
1
U
{\displaystyle (D+L)^{-1}U\!\,}
is less than 1 i.e.
ρ
(
(
D
+
L
)
−
1
U
)
<
1
{\displaystyle \rho ((D+L)^{-1}U)<1\!\,}
The eigenvalues of the matrix
(
D
+
L
)
−
1
U
=
[
0
1
−
1
0
−
1
2
1
0
−
1
2
0
]
{\displaystyle (D+L)^{-1}U={\begin{bmatrix}0&1&-1\\0&-{\frac {1}{2}}&1\\0&-{\frac {1}{2}}&0\end{bmatrix}}\!\,}
are
λ
=
0
,
−
1
4
+
7
i
4
,
−
1
4
−
7
i
4
{\displaystyle \lambda =0,-{\frac {1}{4}}+{\frac {{\sqrt {7}}i}{4}},-{\frac {1}{4}}-{\frac {{\sqrt {7}}i}{4}}\!\,}
Therefore, the spectral radius is
ρ
=
1
16
+
7
16
=
2
2
≈
.7
<
1
{\displaystyle \rho ={\sqrt {{\frac {1}{16}}+{\frac {7}{16}}}}={\frac {\sqrt {2}}{2}}\approx .7<1\!\,}
.
Comparison of Methods
edit
In general, the Gauss-Seidel iteration converges faster than the Jacobi iteration since Gauss-Seidel uses the new components of
x
i
{\displaystyle x_{i}\!\,}
as they become available, but in this case
ρ
J
a
c
o
b
i
<
ρ
G
a
u
s
s
−
S
e
i
d
e
l
{\displaystyle \rho _{Jacobi}<\rho _{Gauss-Seidel}\!\,}
, so the Jacobi iteration is faster.
Consider the three-term polynomial recurrence
p
k
+
1
(
x
)
=
(
x
−
μ
k
)
p
k
(
x
)
−
ν
k
2
p
k
−
1
(
x
)
for
k
=
1
,
2
,
…
,
{\displaystyle p_{k+1}(x)=(x-\mu _{k})p_{k}(x)-\nu _{k}^{2}p_{k-1}(x)\quad {\mbox{ for }}k=1,2,\dots ,\!\,}
initialized by
p
0
(
x
)
=
1
and
p
1
(
x
)
=
x
−
μ
0
{\displaystyle p_{0}(x)=1{\mbox{ and }}p_{1}(x)=x-\mu _{0}\!\,}
, where each
μ
k
and
ν
k
{\displaystyle \mu _{k}{\mbox{ and }}\nu _{k}\!\,}
is real and each
ν
k
{\displaystyle \nu _{k}\!\,}
is nonzero.
Prove that each
p
k
(
x
)
{\displaystyle p_{k}(x)\!\,}
is a monic polynomial of degree
k
{\displaystyle k\!\,}
, and for every
n
=
0
,
1
,
…
{\displaystyle n=0,1,\dots \!\,}
, one has
span
{
p
0
(
x
)
,
p
1
(
x
)
,
…
,
p
n
(
x
)
}
=
span
{
1
,
x
,
x
2
,
…
,
x
n
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{n}(x)\right\}=\operatorname {span} \left\{1,x,x^{2},\dots ,x^{n}\right\}\!\,}
We prove the claim by by induction.
p
0
(
x
)
=
1
{\displaystyle p_{0}(x)=1\!\,}
is monic with degree zero, and
span
{
p
0
(
x
)
}
=
span
{
1
}
{\displaystyle \operatorname {span} \left\{p_{0}(x)\right\}=\operatorname {span} \left\{1\right\}\!\,}
.
p
1
(
x
)
=
x
−
μ
0
{\displaystyle p_{1}(x)=x-\mu _{0}\!\,}
is monic with degree one, and
span
{
p
0
(
x
)
,
p
1
(
x
)
}
=
span
{
1
,
x
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x)\right\}=\operatorname {span} \left\{1,x\right\}\!\,}
.
p
2
(
x
)
=
(
x
−
μ
0
)
(
x
−
μ
1
)
−
ν
1
2
{\displaystyle p_{2}(x)=(x-\mu _{0})(x-\mu _{1})-\nu _{1}^{2}\!\,}
is monic with degree 2, and
span
{
p
0
(
x
)
,
p
1
(
x
)
,
p
2
(
x
)
}
=
span
{
1
,
x
,
x
2
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x),p_{2}(x)\right\}=\operatorname {span} \left\{1,x,x^{2}\right\}\!\,}
.
Assume
p
k
−
1
(
x
)
{\displaystyle p_{k-1}(x)\!\,}
is monic with degree
k
−
1
{\displaystyle k-1\!\,}
, and
span
{
p
0
(
x
)
,
p
1
(
x
)
,
…
,
p
k
−
1
(
x
)
}
=
span
{
1
,
x
,
…
,
x
k
−
1
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k-1}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k-1}\right\}\!\,}
.
Also, assume
p
k
(
x
)
{\displaystyle p_{k}(x)\!\,}
is monic with degree
k
{\displaystyle k\!\,}
, and
span
{
p
0
(
x
)
,
p
1
(
x
)
,
…
,
p
k
(
x
)
}
=
span
{
1
,
x
,
…
,
x
k
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k}\right\}\!\,}
.
Then from the hypothesis,
p
k
+
1
(
x
)
=
(
x
−
μ
k
)
p
k
(
x
)
−
ν
k
2
p
k
−
1
(
x
)
{\displaystyle p_{k+1}(x)=(x-\mu _{k})p_{k}(x)-\nu _{k}^{2}p_{k-1}(x)\!\,}
is monic with degree
k
+
1
{\displaystyle k+1\!\,}
.
Also,
span
{
p
0
(
x
)
,
p
1
(
x
)
,
…
,
p
k
+
1
(
x
)
}
=
span
{
1
,
x
,
…
,
x
k
+
1
}
{\displaystyle \operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k+1}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k+1}\right\}\!\,}
since
p
k
+
1
(
x
)
{\displaystyle p_{k+1}(x)\!\,}
is just
x
n
+
1
{\displaystyle x^{n+1}\!\,}
plus a linear combination of lower order terms already in
span
{
1
,
x
,
…
,
x
k
}
{\displaystyle \operatorname {span} \left\{1,x,\dots ,x^{k}\right\}\!\,}
.
We prove the claim by induction.
Let's consider the case
k
=
2
{\displaystyle k=2\!\,}
. We know that
p
2
(
x
)
=
(
x
−
μ
0
)
(
x
−
μ
1
)
−
ν
1
2
.
{\displaystyle p_{2}(x)=(x-\mu _{0})(x-\mu _{1})-\nu _{1}^{2}.\!\,}
The quadratic formula shows that
p
2
(
x
)
{\displaystyle p_{2}(x)\!\,}
has two simple real roots.
Let
r
21
{\displaystyle r_{21}\!\,}
and
r
22
{\displaystyle r_{22}\!\,}
be the zeros of
p
2
(
x
)
{\displaystyle p_{2}(x)\!\,}
. Then, we have (because of sign changes) that
p
2
(
μ
0
)
=
−
ν
1
2
{\displaystyle p_{2}(\mu _{0})=-\nu _{1}^{2}\!\,}
and
lim
x
→
∞
p
2
(
x
)
>
0
{\displaystyle \lim _{x\rightarrow \infty }p_{2}(x)>0\!\,}
⇒
{\displaystyle \Rightarrow \!\,}
there exists
r
22
∈
(
μ
0
,
∞
)
{\displaystyle r_{22}\in (\mu _{0},\infty )\!\,}
such that
p
2
(
r
22
)
=
0
{\displaystyle p_{2}(r_{22})=0\!\,}
.
p
2
(
μ
0
)
=
−
ν
1
2
{\displaystyle p_{2}(\mu _{0})=-\nu _{1}^{2}\!\,}
and
lim
x
→
−
∞
p
2
(
x
)
>
0
{\displaystyle \lim _{x\rightarrow -\infty }p_{2}(x)>0\!\,}
⇒
{\displaystyle \Rightarrow \!\,}
there exists
r
21
∈
(
−
∞
,
μ
0
)
{\displaystyle r_{21}\in (-\infty ,\mu _{0})\!\,}
such that
p
2
(
r
21
)
=
0
{\displaystyle p_{2}(r_{21})=0\!\,}
.
In conclusion,
r
21
<
μ
0
⏟
r
11
<
r
22
.
{\displaystyle r_{21}<\underbrace {\mu _{0}} _{r_{11}}<r_{22}.\!\,}
.
Let
r
k
−
1
,
1
,
…
,
r
k
−
1
,
k
−
1
{\displaystyle r_{k-1,1},\dots ,r_{k-1,k-1}\!\,}
and
r
k
,
1
,
…
,
r
k
,
k
{\displaystyle r_{k,1},\dots ,r_{k,k}\!\,}
be the simple real roots of
p
k
−
1
(
x
)
{\displaystyle p_{k-1}(x)\!\,}
and
p
k
(
x
)
{\displaystyle p_{k}(x)\!\,}
, respectively, such that the roots of
p
k
{\displaystyle p_{k}\!\,}
are interlaced with the roots of
p
k
−
1
{\displaystyle p_{k-1}\!\,}
, i.e.,
r
k
,
1
<
r
k
−
1
,
1
<
r
k
,
2
<
r
k
−
1
,
2
<
⋯
<
r
k
−
1
,
k
−
1
<
r
k
,
k
.
{\displaystyle r_{k,1}<r_{k-1,1}<r_{k,2}<r_{k-1,2}<\dots <r_{k-1,k-1}<r_{k,k}.\!\,}
Then, we have that
p
k
+
1
(
r
k
,
1
)
=
(
r
k
,
1
−
μ
k
)
p
k
(
r
k
,
1
)
⏟
0
−
ν
k
2
p
k
−
1
(
r
k
,
1
)
=
−
ν
k
2
p
k
−
1
(
r
k
,
1
)
,
{\displaystyle p_{k+1}(r_{k,1})=(r_{k,1}-\mu _{k})\underbrace {p_{k}(r_{k,1})} _{0}-\nu _{k}^{2}p_{k-1}(r_{k,1})=-\nu _{k}^{2}p_{k-1}(r_{k,1}),\!\,}
p
k
+
1
(
r
k
,
2
)
=
(
r
k
,
2
−
μ
k
)
p
k
(
r
k
,
2
)
⏟
0
−
ν
k
2
p
k
−
1
(
r
k
,
2
)
=
−
ν
k
2
p
k
−
1
(
r
k
,
2
)
.
{\displaystyle p_{k+1}(r_{k,2})=(r_{k,2}-\mu _{k})\underbrace {p_{k}(r_{k,2})} _{0}-\nu _{k}^{2}p_{k-1}(r_{k,2})=-\nu _{k}^{2}p_{k-1}(r_{k,2}).\!\,}
From the induction hypothesis
p
k
−
1
(
r
k
,
1
)
{\displaystyle p_{k-1}(r_{k,1})\!\,}
and
p
k
−
1
(
r
k
,
2
)
{\displaystyle p_{k-1}(r_{k,2})\!\,}
have different signs since
r
k
,
1
<
r
k
−
1
,
1
<
r
k
,
2
{\displaystyle r_{k,1}<r_{k-1,1}<r_{k,2}\!\,}
. Then, there exists
r
k
+
1
,
1
∈
(
r
k
,
1
,
r
k
,
2
)
{\displaystyle r_{k+1,1}\in (r_{k,1},r_{k,2})\!\,}
such that
p
k
+
1
(
r
k
+
1
,
1
)
=
0
{\displaystyle p_{k+1}(r_{k+1,1})=0\!\,}
.
Proceeding in the same manner for all the intervals
(
r
k
,
j
,
r
k
,
j
+
1
)
{\displaystyle (r_{k,j},r_{k,j+1})\!\,}
, we obtain that there exists
r
k
+
1
,
j
∈
(
r
k
,
j
,
r
k
,
j
+
1
)
{\displaystyle r_{k+1,j}\in (r_{k,j},r_{k,j+1})\!\,}
such that
p
k
+
1
(
r
k
+
1
,
j
)
=
0
{\displaystyle p_{k+1}(r_{k+1,j})=0\!\,}
for
j
=
1
,
…
,
k
−
1.
{\displaystyle j=1,\dots ,k-1.\!\,}
We now consider the smallest and largest roots of
p
k
+
1
{\displaystyle p_{k+1}\!\,}
i.e.
r
k
+
1
,
0
,
r
k
+
1
,
k
+
1
{\displaystyle r_{k+1,0},r_{k+1,k+1}\!\,}
Since for
k
=
0
,
1
,
2
,
…
{\displaystyle k=0,1,2,\ldots \!\,}
,
p
k
{\displaystyle p_{k}\!\,}
is a monic polynomial,
lim
x
→
∞
p
k
(
x
)
=
∞
{\displaystyle \lim _{x\rightarrow \infty }p_{k}(x)=\infty \!\,}
Hence for any
x
>
r
k
{\displaystyle x>r_{k}\!\,}
(any
x
{\displaystyle x\!\,}
larger than the largest root of
p
k
{\displaystyle p_{k}\!\,}
)
p
k
(
x
)
>
0
{\displaystyle p_{k}(x)>0\!\,}
Therefore
lim
x
→
∞
p
k
+
1
(
x
)
>
0
{\displaystyle \lim _{x\rightarrow \infty }p_{k+1}(x)>0\!\,}
and
p
k
+
1
(
r
k
,
k
)
=
−
ν
k
2
p
k
−
1
(
r
k
,
k
)
<
0
,
{\displaystyle p_{k+1}(r_{k,k})=-\nu _{k}^{2}p_{k-1}(r_{k,k})<0,\!\,}
implies there exists
r
k
+
1
,
k
∈
(
r
k
,
k
,
∞
)
{\displaystyle r_{k+1,k}\in (r_{k,k},\infty )\!\,}
such that
p
k
+
1
(
r
k
+
1
,
k
)
=
0
{\displaystyle p_{k+1}(r_{k+1,k})=0\!\,}
.
If
k
+
1
{\displaystyle k+1\!\,}
is even, then by similar reasoning
lim
x
→
−
∞
p
k
+
1
(
x
)
>
0
{\displaystyle \lim _{x\rightarrow -\infty }p_{k+1}(x)>0\!\,}
and
p
k
+
1
(
r
k
,
1
)
=
−
ν
k
2
p
k
−
1
(
r
k
,
1
)
<
0
,
{\displaystyle p_{k+1}(r_{k,1})=-\nu _{k}^{2}p_{k-1}(r_{k,1})<0,\!\,}
⇒
{\displaystyle \Rightarrow \!\,}
there exists
r
k
+
1
,
0
∈
(
−
∞
,
r
k
,
1
)
{\displaystyle r_{k+1,0}\in (-\infty ,r_{k,1})\!\,}
such that
p
k
+
1
(
r
k
+
1
,
0
)
=
0
{\displaystyle p_{k+1}(r_{k+1,0})=0\!\,}
.
If
k
+
1
{\displaystyle k+1\!\,}
is odd,
lim
x
→
−
∞
p
k
+
1
(
x
)
<
0
{\displaystyle \lim _{x\rightarrow -\infty }p_{k+1}(x)<0\!\,}
and
p
k
+
1
(
r
k
,
1
)
=
−
ν
k
2
p
k
−
1
(
r
k
,
1
)
>
0
,
{\displaystyle p_{k+1}(r_{k,1})=-\nu _{k}^{2}p_{k-1}(r_{k,1})>0,\!\,}
⇒
{\displaystyle \Rightarrow \!\,}
there exists
r
k
+
1
,
0
∈
(
−
∞
,
r
k
,
1
)
{\displaystyle r_{k+1,0}\in (-\infty ,r_{k,1})\!\,}
such that
p
k
+
1
(
r
k
+
1
,
0
)
=
0
{\displaystyle p_{k+1}(r_{k+1,0})=0\!\,}
.
Show that for every
n
=
0
,
1
,
…
{\displaystyle n=0,1,\dots \!\,}
the roots of
p
k
+
1
(
x
)
{\displaystyle p_{k+1}(x)\!\,}
are the eigenvalues of the symmetric tridiagonal matrix
T
=
(
μ
0
ν
1
0
⋯
0
ν
1
μ
1
ν
2
⋱
⋮
0
ν
2
ν
2
⋱
0
⋮
⋱
⋱
⋱
ν
n
0
⋯
0
ν
n
μ
n
)
{\displaystyle T={\begin{pmatrix}\mu _{0}&\nu _{1}&0&\cdots &0\\\nu _{1}&\mu _{1}&\nu _{2}&\ddots &\vdots \\0&\nu _{2}&\nu _{2}&\ddots &0\\\vdots &\ddots &\ddots &\ddots &\nu _{n}\\0&\cdots &0&\nu _{n}&\mu _{n}\end{pmatrix}}}
Let
T
n
+
1
=
(
μ
0
ν
1
0
⋯
0
ν
1
μ
1
ν
2
⋱
⋮
0
ν
2
ν
2
⋱
0
⋮
⋱
⋱
⋱
ν
n
0
⋯
0
ν
n
μ
n
)
{\displaystyle T_{n+1}={\begin{pmatrix}\mu _{0}&\nu _{1}&0&\cdots &0\\\nu _{1}&\mu _{1}&\nu _{2}&\ddots &\vdots \\0&\nu _{2}&\nu _{2}&\ddots &0\\\vdots &\ddots &\ddots &\ddots &\nu _{n}\\0&\cdots &0&\nu _{n}&\mu _{n}\end{pmatrix}}}
Then
det
(
x
I
n
+
1
−
T
n
+
1
)
{\displaystyle \operatorname {det} (xI_{n+1}-T_{n+1})\!\,}
is a monic polynomial in
x
{\displaystyle x\!\,}
of degree
n
+
1
{\displaystyle n+1\!\,}
.
Expanding this determinant along the last row, we have
det
(
x
I
n
+
1
−
T
n
+
1
)
=
(
x
−
μ
n
)
det
(
x
I
n
−
T
n
)
−
ν
n
2
det
(
x
I
n
−
1
−
T
n
−
1
)
(
1
)
{\displaystyle \operatorname {det} (xI_{n+1}-T_{n+1})=(x-\mu _{n})\operatorname {det} (xI_{n}-T_{n})-\nu _{n}^{2}\operatorname {det} (xI_{n-1}-T_{n-1})\qquad (1)\!\,}
where
det
(
x
I
n
−
T
n
)
{\displaystyle \operatorname {det} (xI_{n}-T_{n})\!\,}
and
det
(
x
I
n
−
1
−
T
n
−
1
)
{\displaystyle \operatorname {det} (xI_{n-1}-T_{n-1})\!\,}
are monic polynomials of degree
n
{\displaystyle n\!\,}
and
n
−
1
{\displaystyle n-1\!\,}
, respectively.
Notice that
det
(
x
I
1
−
T
1
)
=
x
−
μ
0
≡
p
1
(
x
)
{\displaystyle \operatorname {det} (xI_{1}-T_{1})=x-\mu _{0}\equiv p_{1}(x)\!\,}
and if we let
p
0
(
x
)
≡
1
{\displaystyle p_{0}(x)\equiv 1\!\,}
, then (1) is equivalent to the three-term recurrence stated in the problem.
Thus,
p
n
+
1
(
x
)
≡
det
(
x
I
n
+
1
−
T
n
+
1
)
=
0
{\displaystyle p_{n+1}(x)\equiv \operatorname {det} (xI_{n+1}-T_{n+1})=0\!\,}
shows that the roots of
p
n
+
1
(
x
)
{\displaystyle p_{n+1}(x)\!\,}
are the eigenvalues of
T
n
+
1
{\displaystyle T_{n+1}\!\,}
.