PrelimNum PDF
PrelimNum PDF
PrelimNum PDF
Wenqiang Feng †
Abstract
This note is intended to assist my prelim examination preparation. You can download and distribute
it. Please be aware, however, that the note contains typos as well as incorrect or inaccurate solutions . At
here, I also would like to thank Liguo Wang for his help in some problems. This note is based on the Dr.
Abner J. Salgado’s lecture note [4]. Some solutions are from Dr. Steven Wise’s lecture note [5].
1
Contents
List of Figures 4
List of Tables 4
1 Preliminaries 5
1.1 Linear Algebra Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.1 Common Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.2 Similar and diagonalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.3 Eigenvalues and Eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.4 Unitary matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.5 Hermitian matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.6 Positive definite matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.7 Normal matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.8 Common Theorems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2 Calculus Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Preliminary Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.4 Norms’ Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.4.1 Vector Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.4.2 Matrix Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2 Direct Method 33
2.1 For squared or rectangular matrices A ∈ Cm,n , m ≥ n . . . . . . . . . . . . . . . . . . . . . . . 33
2.1.1 Singular Value Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.1.2 Gram-Schmidt orthogonalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.1.3 QR Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.2 For squared matrices A ∈ Cn,n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.1 Condition number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.2 LU Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2.3 Cholesky Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.2.4 The Relationship of the Existing Decomposition . . . . . . . . . . . . . . . . . . . . 39
2.2.5 Regular Splittings[3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3 Iterative Method 43
3.1 Diagonal dominant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 General Iterative Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Stationary cases iterative method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3.1 Jacobi Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3.2 Gauss-Seidel Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3.3 Richardson Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.4 Successive Over Relaxation (SOR) Method . . . . . . . . . . . . . . . . . . . . . . . . 50
3.4 Convergence in energy norm for steady cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Dynamic cases iterative method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.5.1 Chebyshev iterative Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.5.2 Minimal residuals Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.5.3 Minimal correction iterative method . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.5.4 Steepest Descent Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.5.5 Conjugate Gradients Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
2
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 3
4 Eigenvalue Problems 63
4.1 Schur algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.2 QR algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.3 Power iteration algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4 Inverse Power iteration algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6 Euler Method 79
6.1 Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2 Trapezoidal Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.3 Theta Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.4 Midpoint Rule Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7 Multistep Methond 88
7.1 The Adams Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.2 The Order and Convergence of Multistep Methods . . . . . . . . . . . . . . . . . . . . . . . . 88
7.3 Method of A-stable verification for Multistep Methods . . . . . . . . . . . . . . . . . . . . . 89
7.4 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
8 Runge-Kutta Methods 95
8.1 Quadrature Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.2 Explicit Runge-Kutta Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.3 Implicit Runge-Kutta Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
8.4 Method of A-stable verification for Runge-Kutta Method . . . . . . . . . . . . . . . . . . . . 96
8.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
References 113
Appendices 114
Appendix 114
Page 3 of 236
A Numerical Mathematics Preliminary Examination Sample Question, Summer, 2013 114
A.1 Numerical Linear Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
A.2 Numerical Solutions of Nonlinear Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
A.3 Numerical Solutions of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
A.4 Numerical Solutions of PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
A.5 Supplemental Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
List of Figures
List of Tables
4
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 5
1 Preliminaries
1.1 Linear Algebra Preliminaries
1.1.1 Common Properties
Properties 1.1. (Structure of Matrices) Let A = [Aij ] be a square or rectangular matrix, A is called
Properties 1.2. (Type of Matrices) Let A = [Aij ] be a square or rectangular matrix, A is called
Properties 1.3. (Properties of invertible matrices) Let A be n × n square matrix. If A is invertible , then
Properties 1.4. (Properties of conjugate transpose) Let A, B be n × n square matrix and γ be a complex
constant, then
Page 5 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 6
• A∗ = A−1 , • A∗ = I,
• A∗ is unitary,
• A is an isometry.
• A is diagonalizable,
• A is unitarily similar to a diagonal matrix , • the column vectors of A form an orthonormal
• the row vectors of A form an orthonormal set, set.
Properties 1.8. (Properties of positive definite Matrices) Let A ∈ Cn×n be a positive definite Matrix and
B ∈ Cn×n , then
Properties 1.9. (Properties of determinants) Let A, B be n × n square matrix and α be a real constant,
then
Properties 1.10. (Properties of inverse) Let A, B be n × n square matrix and α be a real constant, then
Page 6 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 7
Theorem 1.1. (Similar) A is said to be similar to B, if there is a nonsingular matrix X, such that
A = XBX −1 , (A ∼ B).
Theorem 1.2. (Diagonalizablea ) A matrix is diagonalizable , if and only if there exist a nonsingular
matrix X and a diagonal matrix D such that A = XDX −1 .
a Being diagonalizable has nothing to do with being invertible.
Theorem 1.3. (Diagonalizable) A matrix is diagonalizable , if and only if all its eigenvalues are semisimple
.
Theorem 1.4. (Diagonalizable) Suppose dim(A) = n. A is said to be diagonalizable , if and only if A has
n linearly independent eigenvectors .
Corollary 1.1. (Sample question #2, summer, 2013 ) Suppose dim(A) = n. If A has n distinct eigenvalues
, then A is diagonalizable .
Proof. (Sketch) Suppose n = 2, and let λ1 , λ2 be distinct eigenvalues of A with corresponding eigenvectors
v1 , v2 . Now, we will use contradiction to show v1 , v2 are lineally independent. Suppose v1 , v2 are lineally
dependent, then
c1 v1 + c2 v2 = 0, (1)
c1 λ1 v1 + c2 λ1 v2 = 0. (3)
Since λ1 , λ2 and v2 , 0, then c2 = 0. Similarly, we can get c1 = 0. Hence, we get the contradiction.
A similar argument gives the result for n. Then we get A has n linearly independent eigenvectors .
Theorem 1.5. (Diagonalizable) Every Hermitian matrix is diagonalizable , In particular, every real
symmetric matrix is diagonalizable.
Theorem 1.7. The eigenvalues of a triangular matrix are the entries on its main diagonal.
Page 7 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 8
Theorem 1.8. Let A be square matrix with eigenvalue λ and the corresponding eigenvector x.
• λn , n ∈ Z is an eigenvalue of An with corresponding eigenvector x,
• if A is invertible, then 1/λ is an eigenvalue of A−1 with corresponding eigenvector x.
Theorem 1.9. Let A be n × n square matrix and let λ1 , λ2 , · · · , λm be distinct eigenvalues of A with corre-
sponding eigenvectors v1 , v2 , · · · , vm . Then v1 , v2 , · · · , vm are linear independent.
A∗ A = I.
a A matrix A ∈ Rn×n is said to be orthogonal , if
AT A = I.
Theorem 1.10. (Angle preservation) A matrix is unitary , then the transformation defined by A preserves
angles.
<x,y>
Proof. For any vectors x, y ∈ Cn that is angle θ is determined from the inner product via cos θ = kxkky k
.
Since A is unitary (and thus an isometry), then
Theorem 1.11. (Angle preservation) A matrix is real orthogonal , then A has the transformation form
T (θ ) for some θ " # " #
1 0 cos(θ ) sin(θ )
A= T (θ ) = (5)
0 −1 sin(θ ) − cos(θ )
Finally, we can easily establish the diagonalzableility of the unitary matrices.
Theorem 1.12. (Shur Decomposition) A matrix A ∈ Cn×n is similar to a upper triangular matrix and
A = U T U −1 , (6)
Theorem 1.13. (Spectral Theorem for Unitary matrices) A is unitary , then A is diagonalizable and A is
unitarily similar to a diagonal matrix .
A = U DU −1 = U DU ∗ , (7)
Page 8 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 9
Note: this representation is often called the Spectral Representation or Spectral Decomposition of A.
A∗ = A.
Definition 1.4. Let A be Hermitian , then the different eigenvector are orthogonal i.e.
Proof. Let λ1 , λ2 be the arbitrary two different eigenvalues with corresponding eigenvector v1 , v2 . Then
Theorem 1.15. (Spectral Theorem for Hermitian matrices) A is Hermitian , then A is unitary diagonal-
izable .
A = U DU −1 = U DU ∗ , (14)
Theorem 1.16. If A, B are unitarily similar , then A is Hermitian if and only if B is Hermitian .
Page 9 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 10
Proof. Since A, B are unitarily similar , then A = U BU −1 , where U is a unitary matrix . And
∗
A∗ = U −1 B∗ U ∗ = U ∗−1 B∗ U ∗ = U B∗ U −1 ,
U BU −1 = A = A∗ = U B∗ U −1 .
Hence, B = B∗ .
x = α1 e1 + α2 e2 + · · · + αn en .
since,
n
X n
X
(x, x ) = ( αi ei , αj ej )
i =1 j =1
n X
X n
= αi ᾱj ei ej
i =1 j =1
n
X
= |αi |2 .
i =1
Therefore,
i.e.
kAxk`2
kAk2 = sup ≤ ρ (A).
x∈Cn kxk` 2
Let k be the index, s.t: |λn | = ρ (A) and x = ek , Ax = Aek = λn ek , so kAxk`2 = |λn | = ρ (A) and
kAxk`2 kAxk`2
kAk2 = sup ≥ = ρ (A).
x∈Cn kxk`2 kxk`2
Page 10 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 11
xT Ax > 0, ∀x , 0.
x∗ Ax > 0, ∀x , 0.
Problem 1.1. (Sample question #1, summer, 2013 ) Suppose A ∈ Cn×n is hermitian and σ (A) ⊂ (0, ∞).
Prove A is Hermitian Positive Defined (HPD).
x∗ Ax = x∗ U DU −1 x = x∗ U DU ∗ x = (U ∗ x )∗ D (U ∗ x ). (15)
Moreover, since σ (A) ⊂ (0, ∞) then x̃∗ D x̃ > 0 for any nonzero x̃. Hence
A∗ A = AA∗ .
Corollary 1.2. Unitary matrix and Hermitian matrix are normal matrices.
Theorem 1.19. A ∈ Cn×n is normal if and only if every matrix unitarily equivalent to A is normal.
Theorem 1.20. A ∈ Cn×n is normal if and only if every matrix unitarily equivalent to A is normal.
Page 11 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 12
Theorem 1.21. (Spectral theorem for normal matrices) If A ∈ Cn×n has eigenvalues λ1 , · · · , λn , counted
according to multiplicity, the following statements are equivalent.
1. A is normal,
2. A is unitarily diagonalizable,
3. ni=1 nj=1 |aij |2 = ni=1 |λi |2 ,
P P P
Proof. 1. For any ỹ ∈ R(A)⊥ , then ỹ T y = 0, ∀y ∈ R(A). And ∀y ∈ R(A), there exists x, such that Ax = y.
Then
ỹ T Ax = (AT ỹ )T x = 0.
R ( A ) ⊥ ⊂ N ( AT )
Conversely, suppose y ∈ N (AT ), then AT y = 0 and hence (AT y )T x = y T Ax = 0 for any x ∈ Rn . So,
y ∈ R(AT )⊥ . Therefore
N (AT ) ⊂ R(A)⊥
R(A)⊥ = N (AT ),
2. Similarly, we can prove R(AT )⊥ = N (A),
Definition 1.8. (Taylor formula for one variable) Let f (x ) to be n-th differentiable at x0 , then there exists
a neighborhood B(x0 , ), ∀x ∈ B(x0 , ), s.t.
f 00 (x0 ) f ( n ) ( x0 )
f (x ) = f (x0 ) + f 0 (x0 )(x − x0 ) + ( x − x0 ) 2 + · · · + (x − x0 )n + O ((x − x0 )n+1 )
2! n!
f 00 (x0 ) 2 f ( n ) ( x0 ) n
= f (x0 ) + f 0 (x0 )∆x + ∆x + · · · + ∆x + O (∆xn+1 ). (17)
2! n!
Page 12 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 13
Definition 1.9. (Taylor formula for two variables) Let f (x, y ) ∈ C k +1 (B((x0 , y0 ), )), then ∀(x0 +
∆x, y0 + ∆y ) ∈ B((x0 , y0 ), )),
!
∂ ∂
f (x0 + ∆x, y0 + ∆y ) = f (x0 , y0 ) + ∆x + ∆y f (x0 , y0 )
∂x ∂y
!2
1 ∂ ∂
+ ∆x + ∆y f ( x0 , y 0 ) + · · · (18)
2! ∂x ∂y
!k
1 ∂ ∂
+ ∆x + ∆y f ( x 0 , y 0 ) + Rk
k! ∂x ∂y
where
!k +1
1 ∂ ∂
Rk = ∆x + ∆y f (x0 + θ∆x, y0 + θ∆y ), θ ∈ (0, 1).
(k + 1) ! ∂x ∂y
a2 b 2
ab ≤ + , for all a, b ∈ R. (24)
2 2
a2 2
Proof. Since (a − b )2 = a2 − 2ab + b2 ≥ 0, therefore ab ≤ 2 + b2 , for all a, b ∈ R.
b2
ab ≤ a2 + , for all a, b > 0 , > 0. (25)
4
√
Proof. Using Cauchy’s Inequality with 2a, √1 b in place of a, b, we can get the result.
2
Page 13 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 14
ap b q
ab ≤ + for all a, b > 0 . (26)
p q
tp 1
f (t ) = + − t.
p q
We know that the minimum value is at t = 1, since f 0 (t ) = t p−1 = 0 at t = 1. Now, we setting t = ab−q/p ,
we get
(ab−q/p )p 1
0 ≤ f (ab−q/p ) = + − ab−q/p
p q
ap b−q 1
= + − ab−q/p .
p q
So,
ap b−q 1
ab−q/p ≤ + .
p q
ap b q
abq−q/p ≤ + .
p q
1
Since, p + 1q = 1, so pq = p + q and q − q/p = 1. Hence
ap b q
ab ≤ + for all a, b > 0 .
p q
1 1/p
Proof. Using Young’s Inequality with (p )1/p a, p b in place of a, b, we can get the result.
Page 14 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 15
R R R
Proof. Suppose U |u|p dx , 0 and U |v|q dx , 0. Otherwise, if U |u|p dx = 0, then u ≡ 0 a.e. and the Hölder’s
Inequality is trivial. We can use the same argument for v. Now, we define f , g as following
|u| |v|
f = ,g = . (29)
kukLp kvkLq
Now applying Young’s inequality for f g, we have
|u| |v| 1 |u|p 1 |u|q
fg = ≤ p + . (30)
kukLp kvkLq p kukLp q kukq q
L
Hence
Z Z
|uv|dx ≤ |u||v|dx ≤ kukLp kvkLq . (33)
U U
So uv ∈ L1 (U ).
Z Z Z
|uv|dx ≤ |u||v|dx ≤ kvkL∞ (U ) |u|dx = kukL1 (U ) kvkL∞ (U ) . (36)
U U U
Page 15 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 16
1 1
1= + .
p1 /r p2 /r
2. Induction assumption: Assume the inequality holds for k = n − 1, i.e. Πnk=1 ui ∈ Lr (U ) and
Z
r
|u1 · · · un−1 |r dx ≤ Πkn−1
=1 kui kLpk (U ).
U
1 1 1
+ ··· + = .
p1 pn−1 p
so
1 1 1
+ = .
p pn r
From the Hölder’s inequality for n = 2 and the induction assumption, we have
Z Z ! pr Z ! pr
n
r p pn
|u1 · · · un | dx ≤ |u1 · · · un−1 | dx |un | dx
U U U
≤ ku1 krLp (U ) ku2 krLpn (U ) = Πnk=1 kui krLpk (U ).
Page 16 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 17
for k = 1, 2, · · · , n − 1.
Definition 1.17. (Discrete Hölder’s Inequality) Let 1 < p, q < ∞, p1 + 1q = 1. Then for all ak , bk ∈ Rn ,
n
n
1/p n 1/q
X X X
p q
ak bk ≤
|ak |
|bk | . (39)
k =1 k =1 k =1
P p P q
Proof.
P p The idea of proof is same to the integral version. Suppose |ak | , 0 and |bk | , 0. Otherwise, if
|ak | = 0, then ak ≡ 0 and the Hölder’s Inequality is trivial. We can use the same argument for bk . Now,
we define f , g as following
ak b
fk = , gk = k . (40)
kak`p kbk`q
Now applying Young’s inequality for f g, we have
p q
ak bk 1 ak 1 bk
fk gk = ≤ p + . (41)
kak`p kbk`q p kak p q kbkqq
` `
Taking summation yields
∞ ∞
X X ak bk
fk gk =
kak`p kbk`q
k =1 k =1
∞ p q
X 1 ak 1 bk
≤ + (42)
p kakpp q kbkqq
k =1 ` `
P∞ p P∞ q
1 k = 1 ak 1 k =1 bk
= +
p kakpp q kbkqq
` `
1 1
= + = 1.
p q
Therefore
∞
X
ak bk ≤ kak`p kbk`q . (43)
k =1
Page 17 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 18
Proof.
n n
X n X n X n X X
!
ak bk ≤ |ak bk | ≤ |ak ||bk | ≤ sup(|bk |)|ak | ≤
|ak | sup |bk | .
(45)
k =1 k∈N k∈N
k =1 k =1 k =1 k =1
R R R
Proof. Suppose U |u|p dx , 0 and U |v|q dx , 0. Otherwise, if U |u|p dx = 0, then u ≡ 0 a.e.. We can use the
same argument for v. Then the Minkowski’s Inequality is trivial. First, We have the following fact
|u + v|p ≤ (|u| + |v|)p ≤ 2p max(|u|p , |v|p ) ≤ 2p (|u|p + |v|p ) < ∞. (49)
Hence u + v ∈ Lp (U ) if u, v ∈ Lp (U ). Let
1 1 p
+ = 1 or q = . (50)
p q p−1
Then, we have the fact that if u + v ∈ Lp then |u + v|p−1 ∈ Lq , since |u + v|p−1 < ∞ and
Z ! 1q Z ! p1 ·(p−1)
p−1 q p−1
|u + v|p−1
q = |u + v| dx = |u + v|p
dx = ku + vkLp < ∞. (51)
L
U U
Page 18 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 19
p−1
Since ku + vkLp , 0, dividing ku + vkLp on both side of (52) yields
P 1− p1
n p
Diving k =1 |ak + bk | on both sides of the above equation, we get
n 1/p n 1/p n 1/p
X X X
p p p
|ak + bk |
≤ |ak | + |bk | .
k =1 k =1 k =1
Page 19 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 20
Definition 1.22. (Integral Minkowski’s Inequality) Let 1 ≤ p < ∞ and u (x, y ) ∈ Lp (U ). Then
Z Z p ! 1 Z Z ! p1
u (x, y )dx dy p ≤ p
|u (x, y )| dy dx. (55)
Z Z p
u (x, y )dx dy
Z Z !p
≤ u (x, y ) dx dy
Z Z !p−1 Z !
= u (x, y ) dx
u (x, y ) dx dy
| {z }
independent on x
Z Z Z !p−1
= u (x, y ) dx
u (x, y ) dxdy
Z Z Z !p−1
= u (x, y ) dx
u (x, y ) dydx (Fubini)
Z Z Z !p−1
= u (x, y ) dx
u (x, y ) dydx
1/q
Z Z Z !p Z !1/p
u (x, y ) dx dy u (x, y )p dy 1 1
= dx ( + = 1)
p q
| {z }
constant
Z Z !p !1/q Z Z !1/p
u (x, y ) dx dy p
= u (x, y ) dy dx
So, we get
Z Z !p Z Z !p !1−1/p Z Z !1/p
u (x, y ) dx dy ≤ u (x, y ) dx dy p
u (x, y ) dy dx.
R R p 1−1/p
dividing u (x, y ) dx dy on both sides of the above equation yields
Z Z !p !1/p Z Z !1/p
u (x, y ) dx dy p
≤ u (x, y ) dy dx.
Page 20 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 21
Definition 1.23. (Differential Version of Gronwall’s Inequality ) Let η (·) be a nonnegative, absolutely
continuous function on [0, T], which satisfies for a.e t the differential inequality
η 0 (t ) ≤ φ (t )η (t ) + ψ (t ), (57)
where φ(t ) and ψ (t ) are nonnegative, summable functions on [0, T]. Then
Rt " Zt #
φ ( s ) ds
η (t ) ≤ e 0 η (0) + ψ (s )ds , ∀0 ≤ t ≤ T . (58)
0
In particular, if
η (t ) = 0, ∀0 ≤ t ≤ T . (60)
Proof. Since
η 0 (t ) ≤ φ(t )η (t ) + ψ (t ), a.e.0 ≤ t ≤ T . (61)
then
η 0 (s ) − φ(s )η (s ) ≤ ψ (s ), a.e.0 ≤ s ≤ T . (62)
Let
Rs
φ(ξ )dξ
f (s ) = η (s )e − 0 . (63)
By product rule and chain rule, we have
df Rs
φ(ξ )dξ
Rs
φ(ξ )dξ
= η 0 (s )e − 0 − η (s )e − 0 φ (s ), (64)
ds Rs
φ(ξ )dξ
= (η 0 (s ) − η (s )φ(s ))e− 0 (65)
Rs
φ(ξ )dξ
≤ ψ (s )e − 0 , a.e.0 ≤ t ≤ T . (66)
Integral the above equation from 0 to t, then we get
Zt Rs Rt Zt Rs
− 0 φ(ξ )dξ − 0 φ(ξ )dξ
η (s )e ds = η (t )e − η (0) ≤ ψ (s )e− 0 φ(ξ )dξ ds,
0 0
i.e.
Rt Z t Rs
φ(ξ )dξ φ(ξ )dξ
η (t )e − 0 ≤ η (0) + ψ (s )e − 0 ds.
0
Therefore
Rt " Z t Rs #
φ(ξ )dξ − φ(ξ )dξ
η (t ) ≤ e 0 η (0) + ψ (s )e 0 ds .
0
Page 21 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 22
Definition 1.24. (Integral Version of Gronwall’s Inequality) Let ξ (·) be a nonnegative, summable func-
tion on [0, T], which satisfies for a.e t the integral inequality
Z t
ξ (t ) ≤ C1 ξ (s )ds + C2 , (67)
0
where C1 , C2 ≥ 0. Then
ξ (t ) ≤ C2 1 + C1 teC1 t , ∀a.e. 0 ≤ t ≤ T . (68)
In particular, if
Z t
ξ ( t ) ≤ C1 ξ (s )ds, ∀a.e. 0 ≤ t ≤ T , (69)
0
ξ (t ) = 0, a.e. (70)
Proof. Let
Z t
η (t ) : = ξ (s )ds, (71)
0
then
η 0 (t ) = ξ (t ). (72)
Since
Z t
ξ (t ) ≤ C1 ξ (s )ds + C2 , (73)
0
so
η 0 (t ) ≤ C1 η (t ) + C2 . (74)
i.e.
η (t ) ≤ C2 teC1 t . (76)
Therefore
Z t
ξ (s )ds ≤ C2 teC1 t . (77)
0
Page 22 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 23
then,
n
a0 X fk
an + 1 ≤ +
+ β . (80)
(1 + γ ) n 1 (1 + γ )n−k +1
k =0
so
a0 f0
a1 ≤ +β . (82)
(1 + γ ) (1 + γ )n−k
2. Induction Assumption: Assume the discrete Gronwall’s inequality is valid for k = n − 1, i.e.
n−1
a0 X fk
an ≤ + β . (83)
(1 + γ )n (1 + γ )n−k
k =0
( 1 + γ ) an + 1 ≤ an + βfn
n−1
a0 X fk
≤ n
+ β + βfn
(1 + γ ) (1 + γ )n−k
k =0
n−1
a0 X fk fn
≤ +β +β (84)
(1 + γ )n (1 + γ ) n−k (1 + γ )n−n
k =0
n
a0 X fk
= +β .
(1 + γ )n (1 + γ )n−k
k =0
1 θ 1−θ
= + . (86)
r p q
Page 23 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 24
1 1
Proof. If 1 ≤ p < r < q then q < r < p1 , hence there exists θ ∈ [0, 1] s.t. 1
r = θ p1 + (1 − θ ) 1q , therefore:
rθ r (1 − θ ) 1 1
1= + = p + q . (88)
p q rθ r (1−θ )
p q
And |u|rθ ∈ L rθ , |u|r (1−θ ) ∈ L r (1−θ ) , since
! rθp ! rθp
Z p Z
|u| rθ rθ
dx = p
|u| dx = kukrθ
Lp (U )
< ∞, (89)
U U
! r (1−θ ) ! r (1−θ )
Z q q
Z q
r (1−θ )
r (1−θ ) r (1−θ ) q
|u| dx = |u| dx = kukLq (U ) < ∞. (90)
U U
Now, we can use Hölder’s inequality for |u|r = |u|rθ |u|r (1−θ ) , i.e.
Z Z
|u|r dx = |u|rθ |u|r (1−θ ) dx
U U
! rθp Z ! r (1−θ )
Z p q q
rθ rθ r (1−θ ) r (1−θ )
≤ |u| dx |u| dx . (91)
U U
r (1−θ )
= kukrθ
Lp (U ) kukLq (U )
. (92)
Therefore
1 1
Proof. If 1 ≤ p < r < q then q < r < p1 , hence there exists θ ∈ [0, 1] s.t. 1
r = θ p1 + (1 − θ ) 1q , therefore:
rθ r (1 − θ ) 1 1
1= + = p + q . (95)
p q rθ r (1−θ )
p q
And |u|rθ ∈ L rθ , |u|r (1−θ ) ∈ L r (1−θ ) , since
! rθp ! rθp
Z p Z
|u| rθ rθ
dx = p
|u| dx = kukrθ
Lp (U )
< ∞, (96)
U U
! r (1−θ ) ! r (1−θ )
Z q q
Z q
r (1−θ )
r (1−θ ) r (1−θ ) q
|u| dx = |u| dx = kukLq (U ) < ∞. (97)
U U
Page 24 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 25
Now, we can use Hölder’s inequality for |u|r = |u|rθ |u|r (1−θ ) , i.e.
Z Z
r
|u| dx = |u|rθ |u|r (1−θ ) dx
U U
! rθp Z ! r (1−θ )
Z p q q
rθ rθ r (1−θ ) r (1−θ )
≤ |u| dx |u| dx . (98)
U U
r (1−θ )
= kukrθ
Lp (U ) kukLq (U )
. (99)
Therefore
1/p−1/r
Let θ = 1/p−1/q , then we get
1/p−1/r 1/r−1/q
1/p−1/q 1/p−1/q
kukLr (U ) ≤ kuk Lp (U )
kuk Lq (U )
. (101)
Theorem 1.23. (1D Dirichlet-Poincaré inequality) Let a > 0, u ∈ C 1 ([−a, a]) and u (−a) = 0, then the
1D Dirichlet-Poincaré inequality is defined as follows
Z a Z a
u (x )2 dx ≤ 4a2 u 0 (x )2 dx.
−a −a
Therefore
Z x
u (x ) 0
≤
u ( ξ ) dξ
Z −a
x
≤ u 0 (ξ ) dξ
Z−aa
≤ u 0 (ξ ) dξ (x ≤ a)
−a
Z a !1/2 Z a !1/2
2
2 0
u (ξ ) dξ
≤ 1 dξ (Cauchy-Schwarz inequality)
−a −a
Z a !1/2
= (2a)1/2 u 0 (ξ )2 dξ .
−a
Therefore
2 Z a
u (x ) ≤ 2a u 0 (ξ )2 dξ.
−a
Page 25 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 26
>a
Theorem 1.24. (1D Neumann-Poincaré inequality) Let a > 0, u ∈ C 1 ([−a, a]) and ū = −a
u (x )dx, then
the 1D Neumann-Poincaré inequality is defined as follows
Z a Z a
u (x ) − ū (x )2 dx ≤ 2a(a − c ) u 0 (x )2 dx.
−a −a
>a
Proof. Since, ū = −a
u (x )dx, then by intermediate value theorem, there exists a c ∈ [−a, a], s.t.
u (c ) = ū (x ).
Therefore
Z x
u (x ) − ū (x ) 0
≤
u (ξ )dξ
c
Z x
≤ u 0 (ξ ) dξ
Zc a
≤ u 0 (ξ ) dξ (x ≤ a)
c
Z a !1/2 Z a !1/2
2
2 0
≤ 1 dξ u (ξ ) dξ
(Cauchy-Schwarz inequality)
c c
Z a !1/2
= (a − c )1/2 u 0 (ξ )2 dξ .
−a
Therefore
2 Z a
u (x ) − ū (x ) ≤ (a − c ) u 0 (ξ )2 dξ.
−a
Page 26 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 27
Definition 1.28. (Vector Norms) A vector norm is a function k·k : Rn 7− R satisfying the following condi-
tions for all x, y ∈ Rn and α ∈ R
1. nonnegative : kxk ≥ 0, (kxk = 0 ⇔ x = 0),
2. homegenity : kαxk = |α| kxk,
3. triangle inequality :
x + y
≤ kxk +
y
, ∀x, y ∈ Rn ,
Definition 1.29. For x ∈ Rn , some of the most frequently used vector norms are
n
X 3. ∞-norm : kxk∞ = max |xi |,
1. 1-norm : kxk1 = |xi |, 1≤i≤n
i =1
v
t n n 1/p
X X
2. 2-norm : kxk2 = |2 , p
|xi 4. p-norm : kxkp = |xi | .
i =1 i =1
Theorem 1.25. (vector 2-norm invariance) Vector 2-norm is invariant under the orthogonal transforma-
tion, i.e., if Q is an n × n orthogonal matrix, then
Proof.
kQxk22 = (Qx )T (Qx ) = xT QT Qx = xT x = kxk22 .
Page 27 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 28
Definition 1.30. (Matrix Norms) A matrix norm is a function k·k : Rm×n 7− R satisfying the following
conditions for all A, B ∈ Rm×n and α ∈ R
1. nonnegative : kxk ≥ 0, (kxk = 0 ⇔ x = 0),
2. homegenity : kαxk = |α| kxk,
3. triangle inequality :
x + y
≤ kxk +
y
, ∀x, y ∈ Rn ,
Definition 1.31. For A ∈ Rm×n , some of the most frequently matrix vector norms are
v
t m X
n Xn
X
1. F-norm : kAkF = 2
|aij | , 3. ∞-norm : kAk ∞ = max |aij |,
1≤i≤m
i =1 i =1 j =1
m
X kAxkp
2. 1-norm : kAk1 = max |aij |, 4. induced-norm : kAkp = sup .
1≤j≤n kxkp
i =1 x∈Rn ,x,0
Proof.
Theorem 1.26. (Matrix 2-norm and Frobenius invariance) (Matrix 2-norm and Frobenius are invariant
under the orthogonal transformation, i.e., if Q is an n × n orthogonal matrix, then
Page 28 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 29
Theorem 1.27. (Neumann Series) Suppose that A ∈ Rn×n . If kAk < 1, then (I − A) is nonsingular and
∞
X
(I − A)−1 = Ak (113)
k =0
with
1 1
≤
(I − A)−1
≤
. (114)
1 + kAk 1 − kAk
N
X N
X +1
(I − A)SN = SN − ASN = Ak − Ak = A0 − AN + 1 = I − AN + 1 .
k =0 k =1
So
and
∞
X
(I − A)−1 = Ak
k =0
3. bounded norm
Since
1 = kIk =
(I − A) ∗ (I − A)−1
.
So,
(1 − kAk)
(I − A)−1
≤ 1 ≤ (1 + kAk)
(I − A)−1
.
Page 29 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 30
Therefore,
1 1
≤
(I − A)−1
≤
.
1 + kAk 1 − kAk
Theorem 1.28. Let A be a nonnegative matrix. then ρ (A) < 1 if only if I − A is nonsingular and (I − A)−1
is nonnegative.
Au = ρ (A)u
or
1
(I − A)−1 u = u.
1 − ρ (A)
since I − A is nonsingular and (I − A)−1 is nonnegative, this show that 1 − ρ (A) > 0, which implies
ρ (A) < 1.
1.5 Problems
Problem 1.2. (Prelim Jan. 2011#2) Let A ∈ Cm×n and b ∈ Cm . Prove that the vector x ∈ Cn is a least
squares solution of Ax = b if and only if r⊥ range(A), where r = b − Ax.
A∗ Ax = A∗ b.
and
(r, Ax ) = (Ax ) ∗ r = x∗ A∗ (b − Ax )
= x∗ (A∗ b − A∗ Ax ))
= 0.
Therefore, r⊥ range(A). The above way is invertible, hence we prove the result. J
Page 30 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 31
Problem 1.3. (Prelim Jan. 2011#3) Suppose A, B ∈ Rn×n and A is non-singular and B is singular. Prove
that
1 kA − Bk
≤ ,
κ (A) kAk
where κ (A) = kAk ·
A−1
, and k·k is an reduced matrix norm.
Solution. Since B is singular, then there exists a vector x , 0, s.t. Bx = 0. Since A is non-singular, then
A−1 is also non-singular. Moreover, A−1 Bx = 0. Then, we have
So
kxk =
(I − A−1 B)x
≤
A−1 A − A−1 B
kxk ≤
A−1
kA − Bk kxk .
Since x , 0, so
1 ≤
A−1
kA − Bk .
1 kA − Bk
≤ ,
A−1
kAk kAk
i.e.
1 kA − Bk
≤ .
κ (A) kAk
J
for any x ∈ Rn .
3. Let b ∈ Rn be given. Prove that x∗ ∈ Rn solves Ax = b if and only if x∗ minimizes the quadratic
function f : Rn → R defined by
1 T
f (x ) = x Ax − xT b.
2
√ √
Solution. √ 1. (a) Obviously, kxkA = xT Ax ≥ 0. When x = 0, then kxkA = xT Ax = 0; when kxkA =
xT Ax = 0, then we have (Ax, x ) = 0, since A is SPD, therefore, x ≡ 0.
√ √ √
(b) kλxkA = λxT Aλx
= λ2
xT Ax = |λ|
xT
Ax = |λ| kxkA .
(c) Next we will show
x + y
A ≤ kxkA +
y
A . First, we would like to show
y T Ax ≤ kxkA
y
.
A
Page 31 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 32
Therefore,
1
∇f (x ) =
2Ax − b = Ax − b.
2
If Ax∗ = b, then ∇f (x∗ ) = Ax∗ − b = 0, therefore x∗ minimizes the quadratic function f. Conversely,
when x∗ minimizes the quadratic function f, then ∇f (x∗ ) = Ax∗ − b = 0, therefore Ax∗ = b.
J
Page 32 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 33
2 Direct Method
2.1 For squared or rectangular matrices A ∈ Cm,n , m ≥ n
2.1.1 Singular Value Decomposition
A = Û Σ̂ V̂ ∗ .
m×n n×n n×n
A= U Σ V∗ .
m×m m×n n×n
Remark 2.1. 1. SVD works for any matrices, spectral decomposition only works for squared matrices.
2. The spectral decomposition A = XΛX −1 works only if A is non-defective matrices.
For a symmetric matrix the following decompositions are equivalent to SVD.
1. Eigen-value decomposition: i.e. A = XΣX −1 . When A is symmetric, the eigenvalues are real and the
eigenvectors can be chosen to be orthonormal and hence X T X = XX T = I i.e.X −1 = XT . The only
difference is that the singular values are the magnitudes of the eigenvalues and hence the column of X
needs to be multiplied by a negative sign if the eigenvalue turns out to be negative to get the singular
value decomposition. Hence, U=X and σi = |λi |.
2. Orthogonal decomposition: i.e. A = P DP T , where P is a unitary matrix and D is a diagonal matrix.
This exists only when matrix A is symmetric and is the same as eigenvalue decomposition.
3. Schur decomposition i.e. A = QSQT , where Q is a unitary matrix and S is an upper triangular
matrix. This can be done for any matrix. When A is symmetric, then S is a diagonal matrix and
again is the same as the eigenvalue decomposition and orthogonal decomposition.
Page 33 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 34
(u, v)
proju (v) = u,
(u, u)
proj0 (v) = 0.
Remark 2.2. 1. This operator projects the vector v orthogonally onto the line spanned by vector u.
2. the projection map proj0 is the zero map, sending every vector to the zero vector.
Definition 2.2. (Gram-Schmidt orthogonalization) The Gram-Schmidt process then works as follows:
u1
u1 = v1 , q1 =
ku1 k
u
u2 = v2 − proju1 (v2 ), q2 = 2
ku2 k
u
u3 = v3 − proju1 (v3 ) − proju2 (v3 ), q3 = 3
ku3 k
u
u4 = v4 − proju1 (v4 ) − proju2 (v4 ) − proju3 (v4 ), q4 = 4
ku4 k
.. ..
. .
k−1
X uk
uk = vk − projuj (vk ), qk = .
kuk k
j =1
P2 = P.
I −P
P = P ∗.
Page 34 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 35
Definition 2.6. (projection with orthonormal basis) If P is a orthogonal projector, then P = P ∗ and P
has SVD, i.e. P = QΣQ∗ . Since an orthogonal projector has some singular values equal to zero (except the
identity map P=I), it is natural to drop the silent columns of Q and use the reduced rather than full SVD,
i.e.
P = Q̂Q̂∗ .
P = I − Q̂Q̂∗ .
Definition 2.8. (Householder reflectors) The householder reflector F is a particular matrix which satisfies
vv ∗
F = I −2 .
kvk
2.1.3 QR Decomposition
A = Q̂ R̂ .
m×n n×n
A= Q R .
m×m m×n
Theorem 2.5. (Existence of QR Decomposition) Every A ∈ Cm×n has full and reduced QR decomposition.
Page 35 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 36
Theorem 2.6. (Uniqueness of QR Decomposition) Each A ∈ Cm×n of full rank has a unique reduced QR
decomposition A = Q̂R̂ with rjj > 0.
f : D → S
Data → Solution
Definition 2.9. (Well posedness ) We say that a problem is well- posed if the solution depends continuously
on the data, otherwise we say it is ill-posed.
Definition 2.10. (absolute condition number) The absolute condition number κ̂ = κ̂ (x ) of the problem f
at x is defined as
f (x + δx ) − f (x )
κ̂ = lim sup .
δ→0 kδxk≤δ kδxk
If f is (Freēchet) differentiable
κ̂ =
Df (x )
.
∂f ∂f
Df (x ) = [ , ] = [1, −1].
∂x1 ∂x2
and
κ̂ =
Df (x )
∞ = 1.
Definition 2.11. (relative condition number) The absolute condition number κ = κ (x ) of the problem f
at x is defined as
kf (x+δx)−f (x)k
kδxk
κ = lim sup
δ→0 kδxk≤δ kf (x)k
kxk
kxk
=
κ̂.
f (x )
Page 36 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 37
Definition 2.12. (condition number of Matrix- Vector Multiplication) The absolute condition of f (x ) =
Ax
kA(x+δx)−A(x)k
kδxk
κ = lim sup
δ→0 kδxk≤δ kA ( x ) k
kxk
kxk
= kAk
.
A(x )
Theorem 2.7. (condition of Matrix- Vector Multiplication) Since, kxk =
A−1 Ax
≤
A−1
kAxk, then
kxk
≤
A−1
. So
kA ( x ) k
κ ≤ kAk
A−1
.
Particularly,
κ = kAk2
A−1
2 .
Definition 2.13. (condition number of Matrix) Let A ∈ Cn×n , invertible, the condition number of A is
κ (A)k·k = kAk
A−1
.
particularly,
σ
κ2 (A) = kAk2
A−1
2 = 1 .
σn
2.2.2 LU Decomposition
Definition 2.14. (LU Decomposition without pivoting) Let A ∈ Cn×n . An LU factorization refers to the
factorization of A, with proper row and/or column orderings or permutations, into two factors, a lower
triangular matrix L and an upper triangular matrix U ,
A = LU .
In the lower triangular matrix all elements above the diagonal are zero, in the upper triangular matrix,
all the elements below the diagonal are zero. For example, for a 3-by-3 matrix A, its LU decomposition
looks like this:
Page 37 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 38
Definition 2.15. (LU Decomposition with partial pivoting) The LU factorization with Partial Pivoting
refers often to the LU factorization with row permutations only,
P A = LU ,
where L and U are again lower and upper triangular matrices, and P is a permutation matrix which, when
left-multiplied to A, reorders the rows of A.
Definition 2.16. (LU Decomposition with full pivoting) An LU factorization with full pivoting involves
both row and column permutations,
P AQ = LU ,
where L, U and P are defined as before, and Q is a permutation matrix that reorders the columns of A
A = L̃D Ũ ,
where D is a diagonal matrix and L and U are unit triangular matrices , meaning that all the entries on
the diagonals of L and U are one.
For
example, for a 3-by-3 matrix A, its LDU decomposition looks like this:
Theorem 2.8. (existence of Decomposition) Any square matrix A admits an LUP factorization. If A is
invertible, then it admits an LU (or LDU) factorization if and only if all its leading principal minors are
nonsingular. If A is a singular matrix of rank k , then it admits an LU factorization if the first k leading
principal minors are nonsingular, although the converse is not true.
Definition 2.18. (Cholesky Decomposition) In linear algebra, the Cholesky decomposition or Cholesky
factorization is a decomposition of a Hermitian, positive-definite matrix into the product of a lower trian-
gular matrix and its conjugate transpose,
A = LL∗ .
Definition 2.19. (LDM Decomposition) Let A ∈ Rn×n and all the leading principal minors det (A(1 :
k; 1 : k )) , 0; k = 1, · · · , n − 1. Then there exist unique unit lower triangular matrices L and M and a unique
diagonal matrix D = diag (d1 , · · · , dn ), such that
A = L̃DM T .
Page 38 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 39
Definition 2.20. (LDL Decomposition) A closely related variant of the classical Cholesky decomposition
is the LDL decomposition,
A = L̃D L̃∗ ,
Remark 2.3. This decomposition is related to the classical Cholesky decomposition, of the form LL∗ , as follows:
∗ 1 1 1 1
A = L̃D L̃ = L̃D 2 D 2 ∗ L̃∗ = L̃D 2 (L̃D 2 )∗ .
The LDL variant, if efficiently implemented, requires the same space and computational complexity to construct
and use but avoids extracting square roots. Some indefinite matrices for which no Cholesky decomposition exists
have an LDL decomposition with negative entries in D. For these reasons, the LDL decomposition may be pre-
ferred. For real matrices, the factorization has the form A = LDLT and is often referred to as LDLT decomposition
(or LDLT decomposition). It is closely related to the eigendecomposition of real symmetric matrices, A = QΛQT .
A = M − N.
The pair of matrices M, N is a regular splitting of A, if M is nonsingular and M −1 and N are nonnegative .
Theorem 2.9. (The eigenvalue radius estimation of Regular Splittings[3]) Let M, N be a regular splitting
of A. Then
ρ (M −1 N ) < 1
Page 39 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 40
2. ⇐: since A,M are nonsingular and A−1 is nonnegative, then A = M (I − G ) is nonsingular. Moreover
−1
A−1 N = M (I − M −1 N ) N
= (I − M −1 N )−1 M −1 N
= (I − G )−1 G.
Clearly, G = M −1 N is nonnegative by the assumptions, and as a result of the Perron-Frobenius theo-
rem, there is a nonnegative eigenvector x associated with ρ (G ) which is an eigenvalue, such that
Gx = ρ (G )x.
Therefore
ρ (G )
A−1 N x = x.
1 − ρ (G )
and this can be true only when 0 ≤ ρ (G ) ≤ 1. Since I − G is nonsingular, then ρ (G ) , 1, which implies
that ρ (G ) < 1.
2.3 Problems
Problem 2.1. (Prelim Aug. 2010#1) Prove that A ∈ Cm×n (m > n) and let A = Q̂R̂ be a reduced QR
factorization.
1. Prove that A has rank n if and only if all the diagonal entries of R̂ are non-zero.
2. Suppose rank(A) = n, and define P = Q̂Q̂∗ . Prove that range(P ) = range(A).
3. What type of matrix is P?
Solution. 1. From the properties of reduced QR factorization, we knowQthat Q̂ has orthonormal columns,
therefore det(Q̂ ) = 1 and R̂ is upper triangular matrix, so det(R̂) = ni=1 rii . Then
n
Y
det(A) = det(Q̂R̂) = det(Q̂ ) det(R̂) = rii .
i =1
Therefore, A has rank n if and only if all the diagonal entries of R̂ are non-zero.
2. (a) range(A) ⊆ range(P ): Let y ∈ range(A), that is to say there exists a x ∈ Cn s.t. Ax = y. Then by
reduced QR factorization we have y = Q̂R̂x. then
therefore y ∈ range(P ).
(b) range(P ) ⊆ range(A): Let v ∈ range(P ), that is to say there exists a v ∈ Cn , s.t. v = P v = Q̂Q̂∗ v.
Claim 2.1.
Page 40 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 41
Proof.
−1
A (A∗ A)−1 A∗
= Q̂R̂ R̂∗ Q̂∗ Q̂R̂ R̂∗ Q̂∗
−1
= Q̂R̂ R̂∗ R̂ R̂∗ Q̂∗
−1
= Q̂R̂R̂−1 R̂∗ R̂∗ Q̂∗
= Q̂Q̂∗ .
J
Therefore by the claim, we have
Problem 2.2. (Prelim Aug. 2010#4) Prove that A ∈ Rn×n is SPD if and only if it has a Cholesky factor-
ization.
Problem 2.3. (Prelim Aug. 2009#2) Prove that for any matrix A ∈ Cn×n , singular or nonsingular, there
exists a permutation matrix P ∈ Rn×n such that PA has an LU factorization, i.e. PA=LU.
Solution. J
Problem 2.4. (Prelim Aug. 2009#4) Let A ∈ Cn×n and σ1 ≥ σ2 ≥ · · · σn ≥ 0 be its singular values.
1. Let λ be an eigenvalue of A. Show that |λ| ≤ σ1 .
Q
2. Show that det(A) = n σ .
j =1 j
Solution. 1. Since σ1 = kAk2 (proof follows by induction), so we need to show |λ| ≤ kAk2 .
|λ| kxk2 = kλxk2 = kAxk ≤ kAk2 kxk2 .
Therefore,
|λ| ≤ σ1 .
Page 41 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 42
2.
Y n
det(A) = det(U ΣV ∗ ) det(U ) det(Σ) det(V ∗ ) = det(Σ) = σj .
j =1
Solution. J
Page 42 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 43
3 Iterative Method
3.1 Diagonal dominant
Definition 3.1. (Diagonal dominant of size δ) A ∈ Cn×n has diagonal dominant of size δ > 0 if
X
|aii | ≥ |aij | + δ.
j,i
Pn
Proof. 1. Let b = Ax and chose k ∈ (1, 2, 3, · · · , n) s.t kxk∞ = |xk |. Moreover, let bk = j =1 akj xj . Since
X
|aii | ≥ |aij | + δ,
j,i
and
X X X
|aij xj | ≤ |aij ||xj | ≤ kxk∞ |aij |.
j,i j,i j,i
Then
X n
|bk | = akj xj
j =1
X n
= akk xk + akj xj
j,k
X n
≥ |akk xk | − akj xj
j,k
X
≥ |akk xk | − kxk∞ |aij |
j,i
X
≥ |akk | kxk∞ − kxk∞ |aij |
j,i
= δ kxk∞ .
So, kAxk∞ = kbk∞ ≥ kbk k∞ ≥ δ kxk∞ . If Ax = 0, then x = 0. So, ker (A) = 0, and then, A−1 exists.
2. Since kAxk∞ = kbk∞ ≥ kbk k∞ ≥ δ kxk∞ , so kAxk ≥ δ kxk∞ and
A−1
∞ ≤ 1δ .
Ax = b, (115)
Page 43 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 44
is a sequence given by
xk +1 = φ(A, b, xk , · · · , xk−r ).
Definition 3.2. (General Iterative Scheme) A general linear two layer iterative scheme reads
x k +1 − x k
!
Bk + Axk = b.
αk
If xk → x0 , then x0 solves Ax = b. So
!
x0 − x0
Bk + Ax0 = b,
αk
i.e.
Ax0 = b.
x k +1 − x k
!
B + Axk = b.
α
Then we get
xk +1 = xk + αB−1 (b − Axk ).
Definition 3.3. (Error Transfer Operator) Let ek = x − xk , where x is exact solution and xk is the approx-
imate solution at k step. Then
x = x + αB−1 (b − Ax )
x k +1 = xk + αB−1 (b − Axk ).
So, we get
After we defined the error transfer operator , the iterative can be written as
xk +1 = T xk + αB−1 b.
Page 44 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 45
Theorem 3.1. (sufficient condition for converges) The sufficient condition for converges is
kT k < 1. (116)
Theorem 3.2. (sufficient & necessary condition for converges) The sufficient & necessary condition for
converges is
ρ (T ) < 1, (117)
A = L + D + U.
Definition 3.5. (Error Transfer Operator for Jacobi Method ) the error transfer operator for Jacobi Method
is as follows
T = I − D −1 A.
A = L + D + U.
and
D xk +1 − xk + Axk = b.
Then we have
D xk +1 − xk + (L + D + U )xk = Lxk + Dxk +1 + U xk = b.
or
1 X
xik +1 = bi − aij xjk .
aii
j,i
Page 45 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 46
Theorem 3.3. (convergence of the Jacobi Method) If A is diagonal dominant , then the Jacobi Method
convergences.
Proof. We want to show If A is diagonal dominant , then
TJ
< 1, then Jacobi Method convergences. From
the definition of T, we know that T for Jacobi Method is as follows
TJ = I − D −1 A.
So,
X X aij
kT k∞ = max |tij | = max | |.
i i aii
j i,j
Therefore,
X |aij | δ
1≥ + .
|aii | |aii |
j,i
Hence, kT k∞ < 1
A = L + D + U.
Definition 3.7. (Error Transfer Operator for Gauss-Seidel Method ) The error transfer operator for Gauss-
Seidel Method is as follows
T = I − (L + D )−1 A
= I − (L + D )−1 (L + D + U )
= −(L + D )−1 U .
Page 46 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 47
Remark 3.2. The Gauss-Seidel method is an iterative technique for solving a square system of n linear equations
with unknown x:
Ax = b.
where the matrix A is decomposed into alower triangular component L∗ , and a strictly upper triangular
component U: A = L∗ + U .
In more detail, write out A, x and b in their components:
Then the decomposition of A into its lower triangular component and its strictly upper triangular component is
given by:
L∗ x = b − U x
The Gauss-Seidel method now solves the left hand side of this expression for x, using previous value for x on the
right hand side. Analytically, this may be written as:
However, by taking advantage of the triangular form of L∗ , the elements of x(k+1) can be computed sequentially
using forward substitution:
(k +1) 1 X (k +1)
X (k )
xi = bi − aij xj − aij xj , i, j = 1, 2, . . . , n.
aii
j<i j>i
The procedure is generally continued until the changes made by an iteration are below some tolerance, such as a
sufficiently small residual.
Theorem 3.4. (convergence of the Gauss-Seidel Method) If A is diagonal dominant , then the Gauss-Seidel
Method convergences.
Proof. We want to show If A is diagonal dominant , then kTGS k < 1, then Gauss-Seidel Method conver-
gences. From the definition of T, we know that T for Gauss-Seidel Method is as follows
Page 47 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 48
So,
X X
|aii | − |aij | ≥ |aij | + δ,
j<i j>i
which implies
P
j>i |aij |
( )
γ = maxi P | ≤ 1.
|aii | − j<i |aij
Moreover
X X
X
X
|((L + D )y )i0 | = | ai0 j yj + ai0 i0 yj | ≥ |ai0 i0 yj | − | ai0 j yj | = |ai0 i0 |
y
∞ − | ai0 j yj | ≥ |ai0 i0 |
y
∞ − |ai0 j |
y
∞ .
j<i0 j<i0 j<i0 j<i0
Therefore, we have
X
X
|ai0 i0 |
y
∞ − |ai0 j |
y
∞ ≤ |ai0 j | kxk∞ ,
j<i0 j>i0
which implies
P
j>i0 |ai0 j |
y
≤
∞
P kxk∞ .
|ai0 i0 | − j<i0 |ai0 j |
So,
kTGS xk∞ ≤ γ kxk∞ ,
which implies
kTGS k∞ ≤ γ < 1.
A = L + D + U.
x k +1 − x k
!
I + Axk = b.
ω
Page 48 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 49
Definition 3.9. (Error Transfer Operator for Gauss-Seidel Method ) The error transfer operator for Gauss-
Seidel Method is as follows
Remark 3.3. Richardson iteration is an iterative method for solving a system of linear equations. Richardson
iteration was proposed by Lewis Richardson in his work dated 1910. It is similar to the Jacobi and Gauss-Seidel
method. We seek the solution to a set of linear equations, expressed in matrix terms as
Ax = b.
where α is a scalar parameter that has to be chosen such that the sequence x(k ) converges.
It is easy to see that the method has the correct fixed points, because if it converges, then x(k +1) ≈ x(k ) andx(k )
has to approximate a solution of Ax = b.
2
Theorem 3.5. (convergence of the Richardson Method) Let A = A∗ > 0 (SPD). If 0 < ω < λmax , then the
Richardson Method convergences. Moreover, the best acceleration parameter is given by
2
ωopt = ,
λmin + λmax
Proof. 1. From the above lemma, we know that the error transform operator is as follows
Let λ ∈ σ (A), then ν := 1−ωλ ∈ σ (T ). From the sufficient and & necessary condition for convergence,
we know if σ (T ) < 1, then Richardson Method convergences, i.e.
|1 − ωλ| < 1,
which implies
ωλmax − 1 = 1 − ωλmin .
Therefore, we get
2
ωopt = .
λmin + λmax
Page 49 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 50
|1 − λmax | |1 − λmin |
1
ω
1 ωopt 1
λmax λmin
A = L + D + U.
x k +1 − x k
!
(ωL + D ) + Axk = b.
ω
Page 50 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 51
Definition 3.11. (Error Transfer Operator for Gauss-Seidel Method ) The error transfer operator for SOR
Method is as follows
Theorem 3.6. (Necessary condition for convergence of the SOR Method) If SOR method convergences,
then 0 < ω < 2.
Proof. If SOR method convergences, then ρ (T ) < 1, i.e |λ| < 1. Let λi are the roots of characteristic polyno-
mial XT (λ) = det (λI − T ) = (−1)n Πni=1 (λ − λi ). Then,
Since λi < 1, so |det (TSOR )| < 1. Since TSOR = −(ωL + D )−1 ((1 − ω−1 )D + U ), then
Theorem 3.7. (convergence of the SOR Method for SPD) If A = A∗ , and 0 < ω < 2, then SOR converges.
Proof. Since
TSOR = −(L + ω−1 D )−1 ((1 − ω−1 )D + U ) = (L + ω−1 D )−1 ((ω−1 − 1)D − U ).
I − TSOR = Q−1 A.
Let (λ, x ) be the eigenpair of T, i.e. T x = λx and y = (I − TSOR )x = (1 − λ)x. So, we have
Moreover,
Page 51 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 52
So, we get
(2ω−1 − 1)(Dy, y ) = (1 − |λ|2 )(Ax, x ).
Since 0 < ω < 2, (Dy, y ) > 0 and (Ax, x ) > 0, so we have
(1 − |λ|2 ) > 0.
Then, we have |λ| < 1.
kxkA = (Ax, x );
Let v k +1 = ek +1 − ek , then
1 k +1
Bv + Aek = 0.
α
Then take the inner product of both sides with v k +1 ,
1
(Bv k +1 , v k +1 ) + (Aek , v k +1 ) = 0.
α
Since
1 k +1 1 1 1
ek = (e + e k ) − (e k +1 − e k ) = (e k +1 + e k ) − v k +1 .
2 2 2 2
Therefore,
1
0 = (Bv k +1 , v k +1 ) + (Aek , v k +1 )
α
1 1 1
= (Bv k +1 , v k +1 ) + (A(ek +1 + ek ), v k +1 ) − (Av k +1 , v k +1 )
α 2 2
1 α k +1 k +1 1 k +1
= ((B − A)v , v ) + (A(e + e ), v k +1 )
k
α 2 2
1 α 1
2
2
= ((B − A)v k +1 , v k +1 ) + (
ek +1
A −
ek
A )
α 2 2
Page 52 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 53
i.e.
2m
k +1
2
k +1
2
k
2
v
+
e
≤
e
.
α 2 A A
Hence
2
2
ek +1
≤
ek
.
A A
and
2
ek +1
→ 0.
A
Definition 3.13.
(Chebyshev
iterative Method) Chebyshev iterative Method is going to choose
α1 , α2 , · · · , αk , s.t.
ek
2 is minimal for
x k +1 − x k
!
+ Axk = b.
αk +1
Theorem 3.9. (convergence of Chebyshev iterative Method) If A = A∗ > 0, then for a given n,
ek
is
minimized by choosing
α0
αk = , t = 1, · · · , n.
1 + ρ0 tk
Where
2 κ (A) − 1 (2k + 1) ∗ 2π
α0 = , ρ0 = 2 , tk = cos( ).
λmin + λmax κ2 ( A ) + 1 2n
Moreover, we have
ρ1k
p
κ (A) − 1
e
, where ρ1 = p 2
ek
≤ 2 0
2 2
.
1 + ρ12k κ2 ( A ) + 1
Page 53 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 54
Definition 3.14. (Minimal residuals Method) Minimal residuals iterative Method is going to choose
α1 , α2 , · · · , αk , s.t. the residuals r k = b − Axk is minimal for
x k +1 − x k
!
+ Axk = b.
αk +1
Theorem 3.10. (optimal αk +1 of minimal residuals iterative Method) The optimal αk +1 of minimal
residuals iterative Method is as follows
(r k , Ar k )
αk +1 =
.
Ar k
2
2
x k +1 − x k
!
+ Axk = b,
αk +1
we get
xk +1 = xk + αk +1 r k .
r k +1 = r k − αk +1 Ar k .
Therefore,
2
r k +1
2
= (r k − αk +1 Ar k , r k − αk +1 Ar k )
2
2
=
r k
2 − 2αk +1 (r k , Ar k ) + αk2+1
Ar k
2 .
(r k , Ar k )
αk +1 =
.
Ar k
2
2
Corollary 3.1. The residual r k +1 of minimal residuals iterative Method is orthogonal to residual r k in
A-norm.
Proof.
(Ar k +1 , r k ) = (r k +1 , Ar k ) = (r k − αk +1 Ar k , Ar k ) = (r k , Ar k ) − αk +1 (Ar k , Ar k ) = 0.
Page 54 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 55
Theorem 3.11. (convergence of minimal residuals iterative Method) The minimal residuals iterative
Method converges for any x0 and
κ (A) − 1
Aek
≤ ρ0n
Ae0
, with ρ0 = 2
2 2
.
κ2 ( A ) + 1
Definition 3.15. (Minimal correction
Method) Minimal correction iterative Method is going to choose
α1 , α2 , · · · , αk , s.t. the correction
wk +1
B (wk = B−1 (b − Axk ) = B−1 r k ,A = A∗ > 0, B = B∗ > 0) is minimal
for
x k +1 − x k
!
B + Axk = b.
αk +1
Page 55 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 56
Theorem 3.12. (optimal αk +1 of minimal correction iterative Method) The optimal αk +1 of minimal
correction iterative Method is as follows
(wk , Awk )
wk
αk +1 = −1 =
A .
(B Awk , Awk )
Awk
−1
B
x k +1 − x k
!
B + Axk = b,
αk +1
we get
xk +1 = xk + αk +1 B−1 r k .
r k +1 = r k − αk +1 AB−1 r k .
(wk , Awk )
αk +1 = .
(B−1 Awk , Awk )
Remark 3.5. Most of time, it’s not easy to compute k·kA , k·kB−1 . We will use the following alternative way to
1
implement the algorithm. let v k = B 2 wk , then from the iterative scheme
x k +1 − x k
!
B + Axk = b,
αk +1
x k +1 − x k
!
+ B−1 Axk = B−1 b.
αk + 1
−Axk +1 + Axk
!
− AB−1 Axk = −AB−1 b.
αk +1
Page 56 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 57
therefore
b − Axk +1 − (b − Axk )
!
+ AB−1 (b − Axk ) = 0,
αk +1
i.e.
r k +1 − r k
!
+ AB−1 r k = 0.
αk +1
w k +1 − w k
!
B + Awk = 0.
αk + 1
Then, we have
w k +1 − w k
!
1 1 1 1
B B
2 2 + AB− 2 B 2 wk = b.
αk + 1
1
Multiplying by B− 2 on both side of the above equation yields
k +1 − w k
!
1 w 1 1 1 1
B2 + B− 2 AB− 2 B 2 wk = B− 2 b,
αk +1
i.e.
v k +1 − v k
!
1 1
B + B− 2 AB− 2 v k = 0.
αk +1
1 1
Since B− 2 AB− 2 > 0, then we minimize
v k +1
2 instead of
wk +1
B . But
2 1 1 1 1
2
wk +1
= (Bwk +1 , wk +1 ) = (B 2 B 2 wk +1 , wk +1 ) = (B 2 wk +1 , B 2 wk +1 ) =
v k +1
.
B 2
Theorem 3.13. (convergence of minimal correction iterative Method) The minimal correction iterative
Method converges for any x0 and
κ (B−1 A) − 1
≤ ρ0n
Ae0
B−1 , with ρ0 = 2 −1
Aek
.
B−1 κ2 ( B A ) + 1
• compute xk +1 = xk + αk +1 wk
Page 57 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 58
x k +1 − x k
!
+ Axk = b.
αk +1
Theorem 3.14. (optimal αk +1 of Steepest Descent iterative Method) The optimal αk +1 of Steepest Descent
iterative Method is as follows
2
2
Aek
r k
2 2
αk +1 =
=
2 .
Aek
2 k
r
A A
x k +1 − x k
!
+ Axk = b = Ax,
αk +1
we get
ek +1 = ek + αk +1 Aek .
Therefore
2
ek +1
= (Aek +1 , ek +1 )
A
= (Aek + αk +1 A2 ek , ek + αk +1 Aek )
2
2
2
=
ek
− 2α
Aek
+ α 2
Aek
k +1 k +1
A 2 A
2
2
Aek
r k
2 2
αk +1 =
=
2 .
Aek
2 k
r
A A
Theorem 3.15. (convergence of Steepest Descent iterative Method) The Steepest Descent iterative Method
converges for any x0 (A = A∗ > 0, B = B∗ > 0) and
κ (A) − 1
ek
≤ ρ0n
e0
, with ρ0 = 2
A A
.
κ2 (A) + 1
Page 58 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 59
Definition 3.17. (Conjugate Gradients Method) Conjugate Gradients Method iterative Method
is a three-
layer iterative method which is going to choose α1 , α2 , · · · , αk and τ1 , τ2 , · · · , τk , s.t. the error
ek +1
A is
minimal for
1
Φ (x ) = (Ax, x ) − (f , x ).
2
In fact, the minimum value of Φ is − 12 (A−1 f , f ) at x = A−1 f and the residual r k is the negative gradient of
Φ at xk , i.e.
r k = −∇Φ (xk ).
• Richardson method is always using the increment along the negative gradient of Φ to correct the
result, i.e.
x k + 1 = x k + αk r k .
• Conjugate Gradients Method is using the increment along the direction pk which is not parallel to
the gradient of Φ to correct the result.
Definition 3.18. (A-Conjugate) The direction {pk } is call A-Conjugate, if (pj , Apk ) = 0 when j , k. In
particular,
(pk +1 , Apk ) = 0, ∀k ∈ N.
Let p0 , p1 , · · · , pm be the linearly independent series and x0 be the initial guess, then we can construct
the following series
xk +1 = xk + αk pk , 0 ≤ k ≤ m.
(r k , p k )
αk = .
(pk +1 , Apk )
Page 59 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 60
• compute xk +1 = xk + αk pk
• compute r k +1 = r k − αk Apk
2
(r k +1 ,Apk )
r k +1
• compute βk +1 = − (pk ,Apk )
= − (pk ,Apk2)
• compute xk +1 = xk + βk +1 pk
Properties 3.2. (properties of {pk } and {r k }) the {pk } and {r k } come from the Conjugate Gradients method
have the following properties:
• (pj , r j ) = 0, 0 ≤ i < j ≤ k
• (pi , Apj ) = 0, i , j 0 ≤ i, j ≤ k
• (r i , r j ) = 0, i , j 0 ≤ i, j ≤ k
Theorem 3.16. (convergence of Conjugate Gradients iterative Method) The Conjugate Gradients iterative
Method converges for any x0 (A = A∗ > 0, B = B∗ > 0) and
p
κ (A) − 1
e
≤ 2ρ0
e
, with ρ0 = p 2
k n 0
A A
.
κ2 ( A ) + 1
Definition 3.19. (Krylov subspace) In linear algebra, the order-k Krylov subspace generated by an n-by-n
matrix A and a vector b of dimension n is the linear subspace spanned by the images of b under the first k-1
powers of A (starting from A0 = I), that is,
Theorem 3.17. (Conjugate Gradients iterative Method in Krylov subspace) For Conjugate Gradients
iterative Method, we have
Page 60 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 61
3.6 Problems
Problem 3.1. (Prelim Jan. 2011#1) Consider a linear system Ax = b with A ∈ Rn×n . Richardson’s method
is an iterative method
Mxk +1 = N xk + b
1
Solution. 1. Since M = w,N = M − A = w1 I − A, then we have
xk +1 = (I − wA)xk + bw.
So TR = I − wA, From the sufficient and & necessary condition for convergence, we should have
ρ (TR ) < 1. Since λi are the eigenvalue of A, then we have 1 − λi w are the eigenvalues of TR . Hence
Richardson’s method converges if and only if |1 − λi w| < 1, i.e
λn w − 1 = 1 − λ1 w,
i,e
2
w0 = .
λ1 + λn
J
|1 − λn | |1 − λ1 |
1
w
1 wopt 1
λn λ1
Page 61 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 62
Problem 3.2. (Prelim Aug. 2010#3) Suppose that A ∈ Rn×n is SPD and b ∈ Rn is given. Then nth Krylov
subspace us defined as
D E
Kn := b, Ab, A2 b, · · · , Ak−1 b .
n on−1
Let xj , x0 = 0, denote the sequence of vectors generated by the conjugate gradient algorithm. Prove
j =0
that if the method has not already converged after n − 1 iterations, i.e. rn−1 = b − Axn−1 , 0, then the nth
iterate xn us the unique vector in Kn that minimizes
2
φ(y ) =
x∗ − y
A ,
where x∗ = A−1 b.
Solution. J
Solution. J
Page 62 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 63
4 Eigenvalue Problems
Definition 4.1. (Ger šchgorin disks) Let A ∈ Cn×n , the Ger šchgorin disks of A are
X
Di = {ξ ∈ C : |ξ − aii | < Ri } where Ri = |aij |.
i,j
Theorem 4.1. Every eigenvalue of A lies within at least one of the Ger šchgorin discs Di
Proof. Let λ be an eigenvalue of A and let x = (xj ) be a corresponding eigenvector. Let i ∈ {1, · · · , n} be
chosen so that |xi | = maxj |xj |. (That is to say, choose i so that xi is the largest (in absolute value) number
in the vector x) Then |xi | > 0, otherwise x = 0. Since x is an eigenvector, Ax = λx, and thus:
X
aij xj = λxi ∀i ∈ {1, . . . , n}.
j
We may then divide both sides by xi (choosing i as we explained, we can be sure that xi , 0) and take the
absolute value to obtain
P
j,i aij xj X aij xj X
|λ − aii | = ≤ xi ≤
|aij | = Ri
xi
j,i j,i
xj ≤ 1
xi for j , i.
Corollary 4.1. The eigenvalues of A must also lie within the Ger šchgorin discs Di corresponding to the
columns of A.
(Ax, x )
R(x ) = .
(x, x )
(Ax, x )
R(x ) = = λ.
(x, x )
Page 63 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 64
Properties 4.1. (properties of Reyleigh Quotient) Reyleigh Quotient has the following properties:
1.
2
∇R(x ) = [Ax − R(x )x ]
(x, x )
2. R(x ) minimizes
f (α ) = kAx − αxk2 .
= 2xi .
Therefore, we have
∂r (x ) 2 (Ax )i xT Ax2xi
= −
∂xi xT x (x T x )2
2
= ((Ax )i − R(x )xi ) .
xT x
Page 64 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 65
Hence
2 2
∇R(x ) = (Ax − R(x )x ) = (Ax − R(x )x ) .
xT x (x, x )
2. let
g (α ) = kAx − αxk22 .
Then,
and
g 0 (α ) = −2(Ax, x ) + 2α (x, x ),
f (α ) = kAx − αxk2
, then g 0 (α ) = 0, i.e.
(Ax, x )
α= = R(x ).
(x, x )
4.2 QR algorithm
−1 −1
Proof. 1. Since Qk Rk = Ak−1 , so Rk = Qk Ak−1 and Ak = Rk Qk = Qk Ak−1 Qk .
Page 65 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 66
Similarly,
∗ −1 ∗ −1
Ak−1 = Qk Ak Qk = Q k Ak Q k = Ak−1 .
v 0 , Av 0 , A2 v 0 , A3 v 0 , · · · .
If we want to prove that this sequence converges to an eigenvector of A, the matrix needs to be such that it has a
unique largest eigenvalue λ1 ,
There is another technical assumption. The initial vector v 0 needs to be chosen such that q1T v 0 , 0 . Otherwise,
if v 0 is completely perpendicular to the eigenvector q1 , the algorithm will not converge.
• compute wk = Av k−1
wk
• compute v k =
kwk k2
• compute λk = R(v k )
Theorem 4.2. (Convergence of power algorithm) If A = A∗ , q1T v 0 , 0 and |λ1 | > |λ2 | ≥ · · · |λm | ≥ 0, then
the convergence to the eigenvector of improved Power iteration algorithm is linear, while the convergence to
the eigenvalue is still quadratic, i.e.
k !
v k − (±q1 )
= O λ2
2 λ1
λ2 2k
!
k
λ − λ1
= O .
2 λ1
Page 66 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 67
αj λ2j qj
P P
2 1
αj λ j Aq j
w = Av = qP = qP
αj2 λ2j αj2 λ2j
αj λ2j qj
P
2
v = qP
αj2 λ2·2
j
···
αj λkj qj
P
k
w = q
P 2 2·(k−1)
αj λj
αj λkj qj
P
k
v = qP .
αj2 λ2·k
j
v k can be rewritten as
αj λkj qj α1 λk1 q1 + j>1 αj λkj qj
P P
k
v = qP =q
αj2 λ2·k α12 λ2k 2 2k
1 + j>1 αj λj
P
j
α j λj k
+
P
α1 λ1k q1 j>1 α1 λ1 qj
= ·
|α1 λk1 |
r
P α 2 λj 2k
1 + j>1 α j λ1 1
k
αj λj
q1 +
P
j>1 α1 λ1 qj
= ±1 r .
αj 2 λj 2k
1+
P
j>1 α1 λ1
Therefore,
X α λ !k !k k !
v k − (±q1 )
≤
j j
qj ≤ C λ2 = O λ2 .
2 α1 λ1 λ1 λ1
j>1
Remark 4.3. This shows that the speed of convergence depends on the gap between the two largest eigenvalues
of A. In particular, if the largest eigenvalue of A were complex (which it can’t be for the real symmetric matrices
we are considering), then λ2 = λ̄1 and the algorithm would not converge at all.
Page 67 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 68
• compute λ0 = R(v 0 )
• compute wk = (A − λk−1 I )−1 v k−1
wk
• compute v k =
kwk k2
• compute λk = R(v k )
Theorem 4.3. (Convergence of power algorithm) If A = A∗ , q1T v 0 , 0 and |λ1 | > |λ2 | ≥ · · · |λm | ≥ 0, If we
update the estimate µ for the eigenvalue with the Rayleigh quotient at each iteration we can get a cubically
convergent algorithm, i.e.
3
v k +1 − (±qJ )
= O
v k − (±qJ )
2 2
k k
λ − λJ
= O |λ − λJ | . 3
2
4.5 Problems
Solution. J
Page 68 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 69
Definition 5.1. (convergence with Order p) An iterative scheme converges with order p>0 if there is a
constant C > 0, such that
|x − xk +1 | ≤ C|x − xk |p . (118)
Definition 5.2. (Bisection method) The method is applicable for solving the equation f(x) = 0 for the real
variable x, where f is a continuous function defined on an interval [a, b] and f(a) and f(b) have opposite
signs i.e. f (a)f (b ) < 0. In this case a and b are said to bracket a root since, by the intermediate value
theorem, the continuous function f must have at least one root in the interval (a, b).
Definition 5.3. (Chord method) The method is applicable for solving the equation f(x) = 0 for the real
variable x, where f is a continuous function defined on an interval [a, b] and f(a) and f(b) have opposite
signs i.e. f (a)f (b ) < 0. Instead of the [a,b] segment halving, we?ll divide it relation f (a) : f (b ), It gives
the approach of a root of the equation
h i−1
x k +1 = x k − η k f (x k ).
where
f (b ) − f (a)
ηk =
b−a
Page 69 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 70
Definition 5.4. (Secant method) The method is applicable for solving the equation f(x) = 0 for the real
variable x, where f is a continuous function defined on an interval [a, b] and f(a) and f(b) have opposite
signs i.e. f (a)f (b ) < 0. Instead of the [a,b] segment halving, we?ll divide it relation f (xk ) : f (xk−1 ), It
gives the approach of a root of the equation
h i−1
x k +1 = x k − η k f (x k ).
where
f (xk ) − f (xk−1 )
ηk =
xk − xk−1
Definition 5.5. (Newton’s method) The method is applicable for solving the equation f(x) = 0 for the real
variable x, where f is a continuous function defined on an interval [a, b] and f(a) and f(b) have opposite
signs i.e. f (a)f (b ) < 0. Instead of the [a,b] segment halving, we?ll divide it relation f 0 (xk ), It gives the
approach of a root of the equation
h i−1
x k +1 = x k − η k f (x k ).
where
η k = f 0 (x k )
Page 70 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 71
Theorem 5.1. (convergence of Newton’s method) Let f ∈ C2 , f (x∗ ) = 0, f 0 (x ) , 0 and f 00 (x∗ ) is bounded
in a neighborhood of x∗ . Provide x0 is sufficient close to x∗ , then newton’s method converges quadratically,
i.e.
2
xk +1 − x∗ ≤ C xk − x∗ .
So,
h i−1 1h i−1
ek +1 = ek + f 0 (xk ) f (xk ) = − f 0 (xk ) f 00 (θ )(ek )2 ,
2
i.e.
f 00 (θ )
e k +1 = − h i (e k )2 ,
0
2 f (x )k
This implies
2
xk +1 − x∗ ≤ C xk − x∗ .
Page 71 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 72
Proof. Let x∗ be the root of F (x ) i.e. F (x∗ )=0. From the Newton’s scheme, we have
h i−1
xk+1 = xk − J(xk ) F(xk )
x∗ = x∗
Therefore, we have
h i−1
x∗ − xk + 1 = x∗ − xk + J(xk ) (F(xk ) − F(x∗ ))
h i−1
= x∗ − xk + J(xk ) F(xk ) − F(x∗ )+J(x∗ )(x∗ − xk ) − J(x∗ )(x∗ − xk )
h i−1 h i−1
= I − J (xk ) J (x∗ ) x∗ − xk − J (xk ) F (xk ) − F(x∗ ) + J(x∗ )(x∗ − xk ) .
So,
h i−1
h i−1
x∗ − xk+1
≤
I − J (xk ) J (x∗ )
x∗ − xk
+
J (xk )
F (xk ) − F(x∗ ) + J(x∗ )(x∗ − xk )
.
(119)
Page 72 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 73
h i−1
Now, we will estimate
I − J (xk ) J (x∗ )
and
F (xk ) − F(x∗ ) + J(x∗ )(x∗ − xk )
.
I − J (xk ) −1 J (x∗ )
h i
h i−1 h i h i−1
=
J (xk ) J (xk ) − J (xk ) J (x∗ )
h i−1
=
J (xk ) J (xk ) − J(x∗ )
(120)
h i−1
≤
J (xk )
J (xk ) − J(x∗ )
h i−1
≤ L
J (xk )
x∗ − xk
.
In the last step of the above equation, we use the J is Lipschitz continuous(If J is not Lipschitz contin-
uous, we can only get the Newton method converges linearly to x∗ ). Since F : Rn → Rn is continuously
differentiable, therefore
Z1
F (b ) = F (a) + J (a + θ (b − a))(b − a)dθ.
0
So
Z 1
F ( xk ) = F (x∗ ) + J(x∗ + θ (xk − x∗ ))(xk − x∗ )dθ
0
Z 1
= F (x∗ ) + J(x∗ + θ (xk − x∗ ))(xk − x∗ ) + J(x∗ )(x∗ − xk ) − J(x∗ )(x∗ − xk )dθ
0
Z 1
= F (x∗ ) − J(x∗ )(x∗ − xk ) + J(x∗ + θ (xk − x∗ ))(xk − x∗ ) + J(x∗ )(x∗ − xk )dθ
0
Hence
Z 1
k k
∗ ∗
F (x ) − F(x ) + J(x )(x − x ) = ∗
J(x∗ + θ (xk − x∗ ))(xk − x∗ ) + J(x∗ )(x∗ − xk )dθ.
0
So,
Z
1
F (xk ) − F(x∗ ) + J(x∗ )(x∗ − xk )
= ∗ k ∗ k ∗ ∗ ∗
J (x + θ (x − x ))(x − x ) + J(x )(x − x )dθ
k
0
Z 1
≤
J (x∗ + θ (xk − x∗ ))(xk − x∗ ) + J(x∗ )(x∗ − xk )
dθ
0
Z 1
≤
J (x∗ + θ (xk − x∗ )) − J(x∗ )
(x∗ − xk )
dθ (121)
0
Z 1
2
≤ Lθ
(x∗ − xk )
dθ
0
1
∗
2
≤ L
(x − xk )
.
2
From (119), (120) and (121), we have
3
h i−1
2
2
x∗ − xk+1
≤ L
J (xk )
x∗ − xk
≤ 3L
[J (x∗ )]−1
x∗ − xk
.
(122)
2
Page 73 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 74
Remark 5.2. From the last step of the above proof process, we can get the condition of . such as, If
1
x∗ − xk
≤
,
L
[J (x∗ )]−1
then
1
x∗ − xk+1
≤
x∗ − xk
.
(123)
2
Theorem 5.7. x is a fixed point of φ and Uδ = {z : |x − z| ≤ δ}. If φ is differentiable on Uδ and q < 1 such
φ0 (z ) ≤ q < 1 for all z ∈ Uδ , then
1. φ(Uδ ) ⊂ Uδ
2. φ is contraction.
5.7 Problems
Problem 5.1. (Prelim Jan. 2011#4) Let f : Ω ⊂ Rn → Rn be twice continuously differentiable. Suppose
x∗ ∈ Ω is a solution of f (x ) = 0, and the Jacobian matrix of f, denoted Jf , is invertible at x∗ .
1. Prove that if x0 ∈ Ω is sufficiently close to x∗ , then the following iteration converges to x∗ :
Solution. Let x∗ be the root of f(x ) i.e. f(x∗ )=0. From the Newton’s scheme, we have
h i−1
xk +1 = xk − J (x0 ) f(xk )
x∗ = x∗
Therefore, we have
h i−1
x∗ − xk +1 = x∗ − xk + J (x0 ) (f(xk ) − f(x∗ ))
h i−1
= x∗ − xk − J(x0 ) J(ξ )(x∗ − xk ).
therefore
J(ξ ) ∗
x∗ − xk +1 ≤ 1 −
0
x − xk
J(x )
From theorem
Page 74 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 75
Theorem 5.8. Suppose J : Rm → Rn×n is a continuous matrix-valued function. If J(x*) is nonsingular, then
there exists δ > 0 such that, for all x ∈ Rm with kx − x∗ k < δ, J(x) is nonsingular and
J (x )−1
< 2
J (x∗ )−1
.
we get
1
x∗ − xk +1 ≤ x∗ − xk .
2
Which also shows the convergence is typically only linear.
J
Problem 5.2. (Prelim Aug. 2010#5) Assume that f : R → R, f ∈ C 2 (R), f 0 (x ) > 0 for all x ∈ R, and
f 00 (x ) > 0, for all x ∈ R.
1. Suppose that a root ξ ∈ R exists. Prove that it is unique. Exhibit a function satisfying the assump-
tions above that has no root.
2. Prove that for any starting guess x0 ∈ R, Newton’s method converges, and the convergence rate is
quadratic.
Solution. 1. Let x1 and x2 are the two different roots. So, f (x1 ) = f (x2 ) = 0, then by Mean value
theorem, we have that there exists η ∈ [x1 , x2 ], such f 0 (η ) = 0 which contradicts f 0 (x ) > 0.
2. example f (x ) = ex .
3. Let x∗ be the root of f (x ). From the Taylor expansion, we know
1
0 = f (x∗ ) = f (xk ) + f 0 (xk )(x∗ − xk ) + f 00 (θ )(x∗ − xk )2 ,
2
where θ is between x∗ and xk . Define ek = x∗ − xk , then
1
0 = f (x∗ ) = f (xk ) + f 0 (xk )(ek ) + f 00 (θ )(ek )2 .
2
so
h i−1 1h i−1
f 0 (xk ) f (xk ) = −(ek ) − f 0 (xk ) f 00 (θ )(ek )2 .
2
From the Newton’s scheme, we have
h i−1
xk +1 = xk − f 0 (xk ) f (xk )
x∗ = x∗
So,
h i−1 1h i−1
ek +1 = ek + f 0 (xk ) f (xk ) = − f 0 (xk ) f 00 (θ )(ek )2 ,
2
i.e.
f 00 (θ )
e k +1 = − h i (e k )2 ,
2 f 0 (x k )
Page 75 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 76
Therefore,
k + 1
f 00 (θ ) k 2 C1 k 2
e ≤ h i (e ) ≤ e .
2 f 0 (xk ) 2C2
This implies
2
xk +1 − x∗ ≤ C xk − x∗ .
Problem 5.3. (Prelim Aug. 2010#4) Let f : Rn → R be twice continuously differentiable. Suppose x∗ is
a isolated root of f and the Jacobian of f at x∗ (J (x∗ )) is non-singular. Determine conditions on so that if
kx0 − x∗ k2 < then the following iteration converges to x∗ :
Solution. J
Problem 5.4. (Prelim Aug. 2009#5) Consider the two-step Newton method
f ( xk ) f (y )
y k = xk − , xk + 1 = y k − 0 k
f 0 ( xk ) f ( xk )
xk + 1 − x ∗ f 00 (x )
lim ∗ ∗
= 0 k ,
k→∞ (yk − x )(xk − x ) f ( xk )
xk +1 − x∗ 1 f 00 (xk )
!
lim = .
k→∞ (xk − x∗ )3 2 f 0 ( xk )
3. Would you say that this method is faster than Newton’s method given that its convergence is cubic?
Solution. 1. First, we will show that if xk ∈ [x − h, x + h], then yk ∈ [x − h, x + h]. By Taylor expansion
formula, we have
1 00
0 = f (x∗ ) = f (xk ) + f 0 (xk )(x∗ − xk ) + f (ξk )(x∗ − xk )2 ,
2!
where ξ is between x and xk . Therefore, we have
1 00
f (xk ) = −f 0 (xk )(x∗ − xk ) − f (ξk )(x∗ − xk )2 .
2!
Plugging the above equation to the first step of the Newton’s method, we have
1 f 00 (ξk ) ∗
y k = xk + ( x ∗ − xk ) + ( x − xk ) 2 .
2! f 0 (xk )
Page 76 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 77
then
1 f 00 (ξk ) ∗
yk − x ∗ = ( x − xk ) 2 . (124)
2! f 0 (xk )
Therefore,
1 f 00 (ξk )
1 f 00 (ξk ) ∗
∗
yk − x = ∗ 2 ∗
( x − x k ) ≤ ( x − x k ) ( x − x k ) .
2! f 0 (xk ) 2 f 0 (x )
k
Since we can choose the initial value very close to x∗ , shah that
00
f (ξ ) (x∗ − x ) ≤ 1
f 0 (x ) k
k
Then, we have that
1
yk − x∗ ≤ (x∗ − xk ) .
2
Hence, we proved the result, that is to say, if xk → x∗ , then yk , ξk → x∗ .
2. Next, we will show if xk ∈ [x − h, x + h], then xk +1 ∈ [x − h, x + h]. From the second step of the Newton’s
Method, we have that
f (yk )
xk + 1 − x ∗ = yk − x∗ −
f 0 ( xk )
1
= ((yk − x∗ )f 0 (xk ) − f (yk ))
f 0 ( xk )
1
= [(yk − x∗ )(f 0 (xk ) − f 0 (x∗ )) − f (yk ) + (yk − x∗ )f 0 (x∗ )]
f 0 ( xk )
By mean value theorem, we have there exists ηk between x∗ and xk , such that
f 0 (xk ) − f 0 (x∗ ) = f 00 ηk (xk − x∗ ),
and by Taylor expansion formula, we have
(yk − x∗ )2 00
f ( yk ) = f (x∗ ) + f 0 (x∗ )(yk − x∗ ) + f ( γk )
2
(y − x∗ )2 00
= f 0 (x∗ )(yk − x∗ ) + k f ( γk ) ,
2
where γ is between yk and x∗ . Plugging the above two equations to the second step of the Newton’s
method, we get
(yk − x∗ )2 00
" #
∗ 1 00 ∗ ∗ 0 ∗ ∗ ∗ 0 ∗
xk + 1 − x = f ηk (xk − x )(yk − x ) − f (x )(yk − x ) − f ( γk ) + ( y k − x ) f ( x )
f 0 ( xk ) 2
(yk − x∗ )2 00
" #
1 00 ∗ ∗
= f η (
k kx − x )( y k − x ) − f ( γ )
k . (125)
f 0 ( xk ) 2
Taking absolute values of the above equation, then we have
#
(yk − x∗ )2 00
"
1
xk +1 − x∗ = 0 00 ∗ ∗
f ηk (xk − x )(yk − x ) − f (γk )
f ( xk ) 2
A
≤ A |xk − x∗ | yk − x∗ + yk − x∗ yk − x∗
2
1 1 5
≤ |xk − x∗ | + |xk − x∗ | = |xk − x∗ | .
2 8 8
Hence, we proved the result, that is to say, if yk → x∗ , then xk +1 , ηk , γk → x∗ .
Page 77 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 78
3. Finally, we will prove the convergence order is cubic. From (215), we can get that
xk + 1 − x ∗ f 00 ηk (yk − x∗ )f 00 (γk )
= − .
(xk − x∗ )(yk − x∗ ) f 0 (xk ) 2(xk − x∗ )f 0 (xk )
xk + 1 − x ∗ f 00 ηk 1 f 00 (ξk ) ∗ f 00 (γk )
= − ( x − x k ) .
(xk − x∗ )(yk − x∗ ) f 0 (xk ) 4 f 0 (xk ) f 0 ( xk )
xk + 1 − x ∗ f 00 (x∗ )
lim = .
k→∞ (xk − x∗ )(yk − x∗ ) f 0 (x ∗ )
1 2 f 0 ( xk )
= .
yk − x∗ (x∗ − xk )2 f 00 (ξk )
Hence
!2
x − x∗ 1 f 00 (x∗ )
lim k +1 ∗ 3 = .
k→∞ (xk − x ) 2 f 0 (x ∗ )
J
Page 78 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 79
6 Euler Method
In this section, we focus on
y 0 = f (t, y ),
y (t0 ) = y0 .
is of order of p ≥ 1 , if
is convergent , if
lim max
y (tn ) − yn
= 0. (130)
h→0 n
is of order 1 .
a You can also use multi-step theorem to derive it.
Page 79 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 80
Theorem 6.2. (The convergence of Forward Euler Method) Forward Euler Method
is convergent.
Therefore,
en+1
≤ ken k + hλ ken k + ch2
= (1 + hλ) ken k + ch2 . (139)
Claim:[2]
c
ken k ≤ h[(1 + hλ)n − 1], n = 0, 1, · · · (140)
λ
Proof for Claim (221): The proof is by induction on n.
1. when n = 0, en = 0, hence ken k ≤ λc h[(1 + hλ)n − 1],
2. Induction assumption:
c
ken k ≤ h[(1 + hλ)n − 1]
λ
3. Induction steps:
en+1
≤(1 + hλ) ken k + ch2 (141)
c
≤ (1 + hλ) h[(1 + hλ)n − 1] + ch2 (142)
λ
c
= h[(1 + hλ)n+1 − 1]. (143)
λ
So, from the claim (221), we get ken k → 0, when h → 0. Therefore Forward Euler Method is convergent
.
Page 80 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 81
yn+1 = yn + hf (tn , yn ),
ξ1 = yn
yn + 1 = yn + hf (tn + 0h, ξ1 ).
is of order 1 .
a You can also use multi-step theorem to derive it.
So,
Theorem 6.4. (The convergence of Backward Euler Method) Backward Euler Method
is convergent.
Page 81 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 82
|f (tn+1 , yn+1 ) − f (tn+1 , y (tn+1 ))| ≤ λ|yn+1 − y (tn+1 )|, λ > 0. (152)
Therefore,
en+1
≤ ken k + hλ
en+1
+ ch2 . (153)
So,
(1 − hλ)
en+1
≤ ken k + ch2 . (154)
≤ cheT T .
So, from the claim (155), we get ken k → 0, when h → 0. Therefore Backward Euler Method is convergent
.
ξ1 = yn
ξ2 = yn + h [0f (tn + 0h, ξ1 ) + 1f (tn + 1h, ξ2 )]
yn + 1 = yn + hf (tn + h, ξ2 ).
1
yn+1 = yn + h[f (tn , yn ) + f (tn+1 , yn+1 )], n = 0, 1, 2, · · · . (156)
2
a Trapezoidal Method Method is a combination of Forward Euler Method and Backward Euler Method.
Page 82 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 83
1
y (tn+1 ) = y (tn ) + h[f (tn , y (tn )) + f (tn+1 , y (tn+1 ))], (157)
2
is of order 2 .
a You can also use multi-step theorem to derive it.
1 2 00
y (tn+1 ) = y (tn ) + hy 0 (tn ) + h y (tn ) + O (h3 ) (158)
2!
y 0 (tn+1 ) = y 0 (tn ) + hy 00 (tn ) + O (h2 ). (159)
So,
1
y (tn+1 ) − y (tn ) + h[f (tn , y (tn )) + f (tn+1 , y (tn+1 ))]
2
1
= y (tn+1 ) − y (tn ) + h[y 0 (tn ) + y 0 (tn+1 )]
2
1 1
= y (tn ) + hy (tn ) + h2 y 00 (tn ) + O (h3 ) − y (tn ) + h[y 0 (tn ) + y 0 (tn ) + hy 00 (tn ) + O (h2 )]
0
(160)
2! 2
= O (h3 ).
1
y (tn+1 ) = y (tn ) + h[f (tn , y (tn )) + f (tn+1 , y (tn+1 ))], (161)
2
is convergent.
1
y (tn+1 ) = y (tn ) + h[f (tn , y (tn )) + f (tn+1 , y (tn+1 ))] + O (h3 ), (162)
2
Subtracting (162) from (156), we get
1
en+1 = en + h[f (tn , yn ) − f (tn , y (tn )) + f (tn+1 , yn+1 ) − f (tn+1 , y (tn+1 ))] + ch3 . (163)
2
Since f is lipschitz continuous w.r.t. the second variable, then
Therefore,
1
en+1
≤ ken k + hλ(ken k +
en+1
) + ch3 . (166)
2
Page 83 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 84
So,
1
1
(1 − hλ)
en+1
≤ (1 + hλ) ken k + ch3 . (167)
2 2
Claim:[2]
n
c 2 1 + 21 hλ
ken k ≤ h 1
− 1 , n = 0, 1, · · · (168)
λ 1 − 2 hλ
Then, we can make h small enough to such that 0 < hλ < 2, then
∞ `
1 + 21 hλ
1 X 1 hλ hλ
= 1+ ≤ = exp .
1 − 12 hλ 1 − 21 hλ `! 1 − 12 hλ 1 − 12 hλ
` =0
Therefore,
n n
c 2 1 + 12 hλ c 2 1 + 21 hλ
c 2
nhλ
ken k ≤ h − 1 ≤ h ≤ h exp . (169)
λ 1 − 21 hλ λ 1 − 12 hλ λ 1 − 12 hλ
This bound of true for every negative integer n such that nh < T . Therefore,
c 2 nhλ c 2 T λ
.
ken k ≤ h exp ≤ h exp
(170)
λ 1 − 12 hλ λ 1 − 12 hλ
So, from the claim (170), we get ken k → 0, when h → 0. Therefore Trapezoidal Method is convergent .
Page 84 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 85
ξ1 = yn
ξ2 = yn + h [θf (tn + 0h, ξ1 ) + (1 − θ )(tn + 1h, ξ2 )]
yn+1 = yn + h[θf (tn + 0h, ξ1 ) + (1 − θ )f (tn + h, ξ2 )]
1 1
yn+1 = yn + hf tn + h, (yn + yn+1 ) . (172)
2 2
1 1
y (tn+1 ) = y (tn ) + hf tn + h, (y (tn ) + y (tn+1 )) . (173)
2 2
is of order 2 .
1 2 00
y (tn+1 ) = y (tn ) + hy 0 (tn ) + h y (tn ) + O (h3 ) (174)
2! !
∂ ∂
f (x0 + ∆x, y0 + ∆y ) = f (x0 , y0 ) + ∆x + ∆y f (x0 , y0 ) + O (h2 ). (175)
∂x ∂y
So,
1 1
y (tn+1 ) − y (tn ) + hf tn + h, (y (tn ) + y (tn+1 ))
2 2
1 2 00
= y (tn ) + hy (tn ) + h y (tn ) + O (h3 ) − y (tn )
0
2! !
1 ∂f (tn , yn ) 1 ∂f (tn , yn ) 2
− h f (tn , yn ) + (tn + h − tn ) + ( (y (tn ) + y (tn+1 )) − yn ) + O (h )
2 ∂t 2 ∂y
Page 85 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 86
1 2 00
= hy 0 (tn ) + h y (tn ) + O (h3 )
2! !
1 2 ∂f (tn , yn ) 1 2 ∂f (tn , yn ) 3
− hf (tn , yn ) + h + h + O (h )
2 ∂t 2 ∂y
!
1 ∂f (tn , yn ) ∂f (tn , yn ) 0
= hy 0 (tn ) + h2 + y (tn )
2! ∂t ∂y
!
0 1 2 ∂f (tn , yn ) 1 2 ∂f (tn , yn ) 3
− hy (tn ) + h + h + O (h )
2 ∂t 2 ∂y
= O (h3 ).
Theorem 6.8. (The convergence of Midpoint Rule Method) Midpoint Rule Method
1 1
y (tn+1 ) = y (tn ) + hf tn + h, (y (tn ) + y (tn+1 )) , (177)
2 2
is convergent.
1 1
y (tn+1 ) = y (tn ) + hf tn + h, (y (tn ) + y (tn+1 )) + O (h3 ), (178)
2 2
Subtracting (178) from (172), we get
1 1 1 1
en+1 = en + h f tn + h, (y (tn ) + y (tn+1 )) − f tn + h, (y (tn ) + y (tn+1 )) + ch3 . (179)
2 2 2 2
Since f is lipschitz continuous w.r.t. the second variable, then
f t + 1 h, 1 (y (t ) + y (t 1 1
n 2 2 n n+1 )) − f tn + h, ( y ( tn ) + y ( tn+1 ))
2 2
1
≤ λ|y − y (tn ) + yn+1 − y (tn+1 )|, λ > 0. (180)
2 n
Therefore,
1
en+1
≤ ken k + hλ(ken k +
en+1
) + ch3 . (181)
2
So,
1
1
(1 − hλ)
en+1
≤ (1 + hλ) ken k + ch3 . (182)
2 2
Claim:[2]
n
c 2 1 + 21 hλ
ken k ≤ h 1
− 1 , n = 0, 1, · · · (183)
λ 1 − 2 hλ
Page 86 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 87
Then, we can make h small enough to such that 0 < hλ < 2, then
∞ `
1 + 12 hλ
1 X 1 hλ hλ
= 1+ ≤ = exp .
1 − 12 hλ 1 − 21 hλ ` =0 `! 1 − 12 hλ 1 − 12 hλ
Therefore,
n n
c 2 1 + 12 hλ c 2 1 + 21 hλ
c 2
nhλ
ken k ≤ h − 1 ≤ h ≤ h exp . (184)
λ 1 − 21 hλ λ 1 − 12 hλ λ 1 − 12 hλ
This bound of true for every negative integer n such that nh < T . Therefore,
c 2 nhλ c 2 T λ
.
ken k ≤ h exp ≤ h exp (185)
λ 1 − 12 hλ λ 1 − 12 hλ
So, from the claim (185), we get ken k → 0, when h → 0. Therefore Midpoint Rule Method is convergent
.
6.5 Problems
Solution. J
Page 87 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 88
7 Multistep Methond
7.1 The Adams Method
where
Z tn + s Z h
−1 −1
bm = h pm (τ )dτ = h pm (tn+s−1 + τ )dτ n = 0, 1, 2, · · · .
tn+s−1 0
t − tn+l
pm ( t ) = Πs−1
l =0,l,m , Lagrange interpolation polynomials .
tn+m − tn+l
(1-step Adams-bashforth)
yn+1 = yn + hf (tn , yn ),
(2-step Adams-bashforth)
3 1
yn + 2 = yn + 1 + h f (tn+1 , yn+1 ) − f (tn , yn ) ,
2 2
(3-step Adams-bashforth)
23 4 5
yn+3 = yn+2 + h f (tn+2 , yn+2 ) − f (tn+1 , yn+1 ) + f (tn , yn ) .
12 3 12
Definition 7.2. (General s-step Method) The general s-step Method a can be written as
s
X s
X
am yn+m = h bm f(tn+m , yn+m ). (187)
m=0 m=0
Page 88 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 89
Theorem 7.1. (s-step method convergent order) The multistep method (187) is of order p ≥ 1 if and only
if there exists c , 0 s.t.
Where,
s
X s
X
ρ (w ) : = am wm and σ (w ) := bm w m . (189)
m=0 m=0
Theorem 7.2. (s-step method convergent order) The multistep method (187) is of order p ≥ 1 if and only
if
1. sm=0 am = 0, (i.e.ρ (1) = 0),
P
Where,
s
X s
X
ρ (w ) : = am wm and σ (w ) := bm w m . (190)
m=0 m=0
Lemma 7.1. (Root Condition) If the roots |λi | ≤ 1 for each i = 1, · · · , m and all roots with value 1 are
simple root then the difference method is said to satisfy the root condition.
Theorem 7.3. (The Dahlquist equivalence theorem) The multistep method (187) is convergent if and only
if
1. consistency: multistep method (187)is order of p ≥ 1 ,
2. stability: the polynomial ρ (w ) satisfies the root condition .
Theorem 7.5. (Dahlquist second barrier) The highest oder of an A-stable multistep method is 2 .
7.4 Problems
Solution. Since the quadrature formula (209) is order of p if it is exact for every f ∈ Pp−1 . we can chose
Page 89 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 90
the simplest basis (1, τ, τ 2 , τ 3 , · · · , τ p−1 ), and the order conditions read that
p
X Z b
bj cjm = τ m w (τ )dτ, m = 0, 1, · · · , p − 1. (191)
j =1 a
use Simpson’s rule to derive a 3-step method. Determine its order and whether it is convergent.
Page 90 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 91
So,
1 ξ2 ξ3
ρ (w ) − σ (w )ln(w ) = ξ 2 + 2ξ − (2 + 2ξ + ξ 2 )(ξ − + ···)
3 2 3
+2ξ +ξ 2
−2ξ +ξ 2 − 23 ξ 3
=
−2ξ 2 +ξ 3 − 32 ξ 4
− 13 ξ 3 + 61 ξ 4 − 19 ξ 5
1
= − ξ 4 + O (ξ 5 ).
2
Therefore, by the theorem
1
ρ (w ) − σ (w )ln(w ) = − ξ 4 + O (ξ 5 ).
2
Hence, this scheme is order of 3.
3. The stability Since,
s
X
ρ (w ) : = am wm = −1 + w2 = (w − 1)(w + 1). (202)
m=0
And w = ±1 are simple root which satisfy the root condition. Therefore, this scheme is stable.
Hence, it is of order 3 and convergent. convergent J
Problem 7.3. Restricting your attention to scalar autonomous y 0 = f (y ), prove that the ERK method with
tableau
0
1 1
2 2
1 1
2 0 2
1 0 0 1
1 1 1 1
6 3 3 6
is of order 4.
Page 91 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 92
Solution. J
where f : [t0 , t ∗ ] × R → R is continuous in its first variable and Lipschitz continuous in its second variable.
Prove that Euler’s method converges.
So,
y (tn+1 ) − y (tn ) − hf (tn , y (tn )) = y (tn ) + hy 0 (tn ) + O (h2 ) − y (tn ) − hf (tn , y (tn ))
= y (tn ) + hy 0 (tn ) + O (h2 ) − y (tn ) − hy 0 (tn ) (204)
= O (h2 ).
Therefore,
en+1
≤ ken k + hλ ken k + ch2
= (1 + hλ) ken k + ch2 .
Claim:[2]
c
ken k ≤ h[(1 + hλ)n − 1], n = 0, 1, · · ·
λ
Proof for Claim (221): The proof is by induction on n.
1. when n = 0, en = 0, hence ken k ≤ λc h[(1 + hλ)n − 1],
2. Induction assumption:
c
ken k ≤ h[(1 + hλ)n − 1]
λ
Page 92 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 93
3. Induction steps:
en+1
≤ (1 + hλ) ken k + ch2
c
≤ (1 + hλ) h[(1 + hλ)n − 1] + ch2
λ
c
= h[(1 + hλ)n+1 − 1].
λ
So, from the claim (221), we get ken k → 0, when h → 0. Therefore Forward Euler Method is convergent .
J
what’s the order of the scheme? Is it a convergent scheme? Is it A-stable? Justify your answers.
So,
ξ2 ξ3
ρ (w ) − σ (w )ln(w ) = ξ 2 + 3ξ − (3 + 3ξ + ξ 2 )(ξ − + ···)
2 3
+3ξ +ξ 2
−3ξ −3ξ 2 −ξ 3
=
+ 23 ξ 2 + 32 ξ 3 + 12 ξ 4
−ξ 3 −ξ 4 − 31 ξ 5
1
= − ξ 2 + O (ξ 3 ).
2
Therefore, by the theorem
1
ρ (w ) − σ (w )ln(w ) = − ξ 2 + O (ξ 3 ).
2
Hence, this scheme is order of 1. The stability Since,
s
X
ρ (w ) : = am wm = −2 + w + w2 = (w + 2)(w − 1). (208)
m=0
And w = −1 or w = −2 which does not satisfy the root condition. Therefore, this scheme is not stable.
Hence, it is also not A-stable.
J
Page 93 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 94
Solution. J
Page 94 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 95
8 Runge-Kutta Methods
8.1 Quadrature Formulas
Definition 8.1. (The Quadrature) The Quadrature is the procedure of replacing an integral with a finite
sum.
Definition 8.2. (The Quadrature Formula) Let w be a nonnegative function in (a,b) s.t.
Zb Z
b j
w (τ )dτ < ∞, τ w (τ )dτ < ∞, j = 1, 2, · · · .
0<
a a
Remark 8.1. The quadrature formula (209) is order of p if it is exact for every f ∈ Pp−1 .
Definition 8.3. (Explicit Runge-Kutta Formulas) Explicit Runge-Kutta is to integral from tn to tn+1 as
follows
Z tn + 1
y ( tn + 1 ) = y (tn ) + f (τ, y (τ ))dτ
tn
Z 1
= y (tn ) + h f (tn + hτ, y (tn + hτ ))dτ
0
Specifically, we have
ξ1 = yn ,
ξ2 = yn + ha21 f (tn , ξ1 )
ξ3 = yn + ha31 f (tn + c1 h, ξ1 ) + ha32 f (tn + c2 h, ξ2 )
..
.
ν−1
X
ξν = yn + h aνi f (tn + ci h, ξi ))
i =1
X ν
yn+1 = yn + h bj f (tn + cj h, ξj )).
j =1
Page 95 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 96
Definition 8.5. (Implicit Runge-Kutta Formulas) Implicit Runge-Kutta use the following scheme
ν
X
ξj = yn + h aj,i f (tn + ci h, ξi ), j = 1, 2, · · · , ν
i =1
X ν
yn + 1 = yn + h bj f (tn + cj h, ξj ).
j =1
Theorem 8.2. necessary & sufficient A necessary & sufficient condition for A-stable Runge-Kutta method
is
r (z ) < 1, z ∈ C− ,
where
r (z ) = 1 + zbT (I − zA)−1 1.
8.5 Problems
Page 96 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 97
where d is dimension.
Theorem 9.1. (Discrete maximum principle) Let A = tridiag{ai , bi , ci }ni=1 ∈ Rn×n be a tridiagional
matrix with the properties that
bi > 0, ai , ci ≤ 0, ai + bi + ci = 0.
Prove the following maximum principle: If u ∈ Rn is such that (Au )i =2,··· ,n−1 ≤ 0, then ui ≤ max{u1 , un }.
ak uk−1 + bk uk + ck uk +1 < 0.
This is contradiction to (Au )i =2,··· ,n−1 < 0. Therefore, If u ∈ Rn is such that (Au )i =2,··· ,n−1 < 0, then
ui < max{u1 , un }.
2. For (Au )i =2,··· ,n−1 = 0:
Since (Au )i =2,··· ,n−1 = 0, so
ak uk−1 + bk uk + ck uk +1 = 0.
Since ak + ck = −bk , so
And ak < 0, ck < 0, uk−1 − uk ≤ 0, uk +1 − uk ≤ 0, so uk−1 = uk = uk +1 , that is to say, uk−1 and uk +1 is also
the maximum points. Bu using the same argument again, we get uk−2 = uk−1 = uk = uk +1 = uk +2 .
Repeating the process, we get
u1 = u2 = · · · = un−1 = un .
Theorem 9.2. (Discrete Poincaré inequality) Let Ω = (0, 1) and Ωh be a uniform grid of size h. If Y ∈ Uh
is a mesh function on Ωh such that Y (0) = 0, then there is a constant C, independent of Y and h, for which
kY k2,h ≤ C
δ̄Y
2,h .
Page 97 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 98
Proof. I consider the following uniform partition (Figure. A1) of the interval (0, 1) with N points.
x1 = 0 x2 xN −1 xN = 1
Then,
X N
(Yi−1 − Yi ) = |YN |.
i =2
and
N N N 1/2 N 1/2
Yi−1 − Yi X 2 X Yi−1 − Yi 2
X X
|YN | ≤ |Yi−1 − Yi | = h ≤ h .
h h
i =2 i =2 i =2 i =2
Therefore
K K
X X Yi−1 − Yi 2
|YK |2 ≤ h2
h
i =2 i =2
K 2
Yi−1 − Yi .
X
2
= h (K − 1) h
i =2
1. When K = 2,
2
2 2 Y1 − Y2
|Y2 | ≤ h
h .
2. When K = 3,
Y1 − Y2 2 Y2 − Y3 2
!
2 2 +
|Y3 | ≤ 2h .
h h
Page 98 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 99
3. When K = N ,
Y1 − Y2 2 Y2 − Y3 2 YN −1 − YN 2
!
2 2
|YN | ≤ (N − 1)h + + · · · + .
h h h
Theorem 9.3. (Von Neumann stability analysis method) For the difference scheme
X
Ujn+1 = n
αp Uj−p ,
p∈N
Page 99 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 100
9.1 Problems
Problem 9.1. (Prelim Jan. 2011#7) Consider the Crank-Nicholson scheme applied to the diffusion equation
∂u ∂2 u
= 2
∂t ∂x
where t > 0, −∞ < x < ∞.
1. Show that the amplification factor in the Von Neumann analysis of the scheme us
1 + 21 z ∆t
g (ξ ) = 1
,z = 2 (cos (∆xξ ) − 1).
1− 2z
∆x2
∆t
Let µ = ∆x2
, then the scheme can be rewrote as
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1 ,
2
i.e.
µ n+1 µ µ n µ
− uj−1 + (1 + µ)ujn+1 − ujn++11 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
By using the Fourier transform, i.e.
ujn+1 = g (ξ )ujn , ujn = eij∆xξ ,
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = 1
,z = 2 (cos (∆xξ ) − 1).
1− 2z
∆x2
∆t
2. since z = 2 ∆x 2 (cos (∆xξ ) − 1), then z < 0, then we have
1 1
1 + z < 1 − z,
2 2
therefore g (ξ ) < 1. Since −1 < 1, then
1 1
z − 1 < z + 1.
2 2
Therefore,
1 + 21 z
g (ξ ) = > −1.
1 − 12 z
hence g (ξ ) < 1. So, the scheme is stable.
J
∆t 1 t∗
where b > 0, µ = (∆x )2
, ∆x = L+ 1 , and ∆t = N. Prove that, under suitable restrictions on µ and ∆x, the
n
error grid function e satisfy the estimate
ken k∞ ≤ t ∗ C ∆t + ∆x2 ,
Solution. Let ū be the exact solution and ūjn = ū (n∆t, j∆x ). Then from Taylor Expansion, we have
∂ n 1 ∂2
ūjn+1 = ūjn + ∆t ūj + (∆t )2 2 ū (ξ1 , j∆x ), tn ≤ ξ1 ≤ tn+1 ,
∂t 2 ∂t
∂ 1 ∂ 3 1 ∂4
n
ūj−1 = ūjn − ∆x ūjn − (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ2 ), xj−1 ≤ ξ2 ≤ xj ,
∂x 6 ∂x 24 ∂x
∂ 1 ∂ 3 1 ∂4
ūjn+1 = ūjn + ∆x ūjn + (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ3 ), xj ≤ ξ3 ≤ xj + 1 .
∂x 6 ∂x 24 ∂x
Then the truncation error T of this scheme is
ūjn+1 − ūjn n
ūj−1 − 2ūjn + ūjn+1 ūjn+1 − ūj−1
n
T = − +b
∆t ∆x2 ∆x
= O (∆t + (∆x )2 ).
Therefore
bµ∆x
ejn+1 = ejn + µ ej−1
n
− 2ejn + ejn+1 − ejn+1 − ej−1
n
+ c∆t (∆t + (∆x )2 ),
2
i.e.
! !
bµ∆x n bµ∆x n
ejn+1 = µ + ej−1 + (1 − 2µ)ejn + µ − ej +1 + c∆t (∆t + (∆x )2 ).
2 2
Then
en+1 ≤ µ + bµ∆x en + (1 − 2µ) en + µ − bµ∆x en + c∆t (∆t + (∆x )2 ).
j 2 j−1 j 2 j + 1
Therefore
en+1
≤ µ + bµ∆x
en
+ (1 − 2µ)
en
+ µ − bµ∆x
en
+ c∆t (∆t + (∆x )2 ).
j ∞ 2 j−1 ∞ j ∞ 2 j + 1 ∞
bµ∆x n bµ∆x n
en+1
≤ µ +
ke k∞ + (1 − 2µ) ken k∞ + µ − ke k∞ + c∆t (∆t + (∆x )2 ).
∞ 2 2
bµ∆x
If 1 − 2µ ≥ 0 and µ − 2 ≥ 0, i.e. µ ≤ 12 and 1 − 12 b∆x > 0, then
! !
bµ∆x bµ∆x
en+1
∞
≤ µ+ ken k∞ + ((1 − 2µ)) ken k∞ + µ − ken k∞ + c∆t (∆t + (∆x )2 )
2 2
= ken k∞ + c∆t (∆t + (∆x )2 ).
Then
ken k∞ ≤
en−1
+ c∆t (∆t + (∆x )2 )
∞
≤
en−2
+ c2∆t (∆t + (∆x )2 )
∞
..
≤ .
≤
e0
+ cn∆t (∆t + (∆x )2 )
∞
= ct ∗ (∆t + (∆x )2 ).
J
kAxk2 ≤ kxk2 ,
kAxk∞ ≤ kxk∞ ,
for any x ∈ Rm , provided µ ≤ 1.(In other words, the scheme may only be conditionally stable in the
max norm.)
n+1 n
u1 u1
n+1 u n
u2 2
un+1 = . and un = . .
.. ..
n+1 umn
um
So, the scheme may be written in the form un+1 = Aun , where A = C −1 B. By using the Fourier
transform, i.e.
ujn+1 = g (ξ )ujn , ujn = eij∆xξ ,
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = ,z = 2 (cos (∆xξ ) − 1).
1 − 12 z ∆x2
Moreover, g (ξ ) < 1, therefore, ρ (A) < 1.
kAxk2 ≤ kAk2 kxk2 = ρ (A) kxk2 ≤ kxk2 .
2. the scheme
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1
2
can be rewritten as
µ n+1 µ n+1 µ n µ
(1 + µ)ujn+1 = u + uj +1 + uj−1 + (1 − µ)ujn + ujn+1 .
2 j−1 2 2 2
then
µ µ µ
(1 − µ) u n + µ u n .
1 + µ ujn+1 ≤ uj−1
n+1 n+1
n
+ u +
j +1 j−1 u + j j +1
2 2 2 2
Therefore
µ
n+1
µ
n+1
µ
n
µ
n
(1 + µ)
ujn+1
≤
uj−1
n
+
uj +1
+
uj−1
+ (1 − µ)
uj
+
uj +1
.
∞ 2 ∞ 2 ∞ 2 ∞ ∞ 2 ∞
i.e.
µ
µ
µ µ
(1 + µ)
un+1
∞ ≤
un+1
∞ +
un+1
∞ + kun k∞ + (1 − µ) kun k∞ + kun k∞ .
2 2 2 2
if µ ≤ 1, then
un+1
≤ kun k∞ ,
∞
i.e.
kAun k∞ ≤ kun k∞ .
a2 (∆t )2 n a∆t
ujn+1 = ujn +
2
uj−1 − 2ujn + ujn+1 − ujn+1 − uj−1
n
,
2(∆x ) 2∆x
for the approximating the solution of the Cauchy problem for the advection equation
∂u ∂u
+a = 0, a > 0.
∂t ∂x
Use Von Neumann’s Method to show that the Lax-Wendroff scheme is stable provided the CFL condition
a∆t
≤ 1.
∆x
is enforced.
then we have
a2 (∆t )2 n a∆t
g (ξ )ujn = ujn + u − 2u n
+ u n
j +1 − u n
− u n
j−1 .
2(∆x )2 j−1 j 2∆x j +1
And then
a2 (∆t )2 i (j−1)∆xξ i (j +1)∆xξ
a∆t
i (j +1)∆xξ i (j−1)∆xξ
g (ξ )eij∆xξ = eij∆xξ + e − 2e ij∆xξ
+ e − e − e .
2(∆x )2 2∆x
Therefore
a2 (∆t )2 −i∆xξ a∆t
g (ξ ) = 1+ 2
e − 2 + ei∆xξ − ei∆xξ − e−i∆xξ
2(∆x ) 2∆x
a2 (∆t )2 a∆t
= 1+ (2 cos(∆xξ ) − 2) − (2i sin(∆xξ ))
2(∆x )2 2∆x
a2 (∆t )2 a∆t
= 1+ (cos(∆xξ ) − 1) − (i sin(∆xξ )) .
(∆x )2 ∆x
a∆t
Let µ = ∆x , then
g (ξ ) = 1 + µ2 (cos(∆xξ ) − 1) − µ (i sin(∆xξ )) .
If g (ξ ) < 1, then the scheme is stable, i,e
2
1 + µ2 (cos(∆xξ ) − 1) + (µ sin(∆xξ ))2 < 1.
i.e.
i.e.
i.e.
i.e
then we get µ < 1. The above process is invertible, therefore, we prove the result. J
Solution. J
Theorem 10.1. (1D Dirichlet-Poincaré inequality) Let a > 0, u ∈ C 1 ([−a, a]) and u (−a) = 0, then the
1D Dirichlet-Poincaré inequality is defined as follows
Z a Z a
u (x )2 dx ≤ 4a2 u 0 (x )2 dx.
−a −a
Therefore
Z x
u (x ) 0
≤
u (ξ )dξ
Z −a
x
≤ u 0 (ξ ) dξ
Z−aa
≤ u 0 (ξ ) dξ (x ≤ a)
−a
Z a !1/2 Z a !1/2
2
≤ 12 dξ u 0 (ξ ) dξ (Cauchy-Schwarz inequality)
−a −a
Z a !1/2
= (2a)1/2 u 0 (ξ )2 dξ .
−a
Therefore
2 Z a
u (x ) ≤ 2a u 0 (ξ )2 dξ.
−a
>a
Theorem 10.2. (1D Neumann-Poincaré inequality) Let a > 0, u ∈ C 1 ([−a, a]) and ū = −a
u (x )dx, then
the 1D Neumann-Poincaré inequality is defined as follows
Z a Z a
u (x ) − ū (x )2 dx ≤ 2a(a − c ) u 0 (x )2 dx.
−a −a
>a
Proof. Since, ū = −a
u (x )dx, then by intermediate value theorem, there exists a c ∈ [−a, a], s.t.
u (c ) = ū (x ).
then by calculus fact, we have
Z x
u (x ) − ū (x ) = u (x ) − u (c ) = u 0 (ξ )dξ.
c
Therefore
Z x
u (x ) − ū (x ) 0
≤
u (ξ )dξ
c
Z x
≤ u 0 (ξ ) dξ
Zc a
≤ u 0 (ξ ) dξ (x ≤ a)
c
Z a !1/2 Z a !1/2
2
2 0
≤ 1 dξ u (ξ ) dξ
(Cauchy-Schwarz inequality)
c c
Z a !1/2
= (a − c )1/2 u 0 (ξ )2 dξ .
−a
Therefore
2 Z a
u (x ) − ū (x ) ≤ (a − c ) u 0 (ξ )2 dξ.
−a
Integration on both sides of the above equation from −a to a w.r.t. x yields
Z a Za Z a
u (x ) − ū (x )2 dx ≤ (a − c ) u 0 (ξ )2 dξdx
−a
Z−aa −a
Za
= u 0 (ξ )2 dξ
(a − c )dx
−a −a
Z a
= 2a(a − c ) u 0 (ξ )2 dξ
Z−aa
= 2a(a − c ) u 0 (x )2 dx.
−a
Definition 10.1. (symmetric, continuous and coercive) We consider the bilinear form a : H × H → R on
a normed space H.
1. a(·, ·) is said to be symmetric provided that
Proof.
Theorem 10.3. (Lax-Milgram Theorem[1]) Given a Hilbert space H, a continuous, coercive bilinear form
a(·, ·) and a continuous functional F ∈ H 0 , there exists a unique u ∈ H s.t.
a(u, v ) = F (v ), ∀v ∈ H.
Theorem 10.4. (Céa Lemma[1]) Suppose V is subspace of H. a(·, ·) is continuous and coercive bilinear
form on V. Given F ∈ V 0 , u ∈ V , s.t.
a(u, v ) = F (v ), ∀v ∈ V .
a(uh , v ) = F (v ), ∀v ∈ Vh ,
we have
C
ku − uh kV ≤ min ku − vkV ,
α v∈Vh
Proof.
Theorem 10.5. (Convergence of 1d FEM) The linear basis FEM solution uh for
−u 00 (x ) = f (x ), x ∈ I = x ∈ [a, b ] ,
u (a) = u (b ) = 0.
Q1 u (x ) = u (xi ) + u 0 (x )(x − xi ).
Then we have
Z
u (x ) − Q1 u (x ) = x − y u 00 (y )dy.
I
This implies
Z
u (x ) − Q1 u (x )
C ( Ii )
= max x − y u 00 (y ) dy
x∈Ii I
Z i
≤ h u 00 (y ) dy
Ii
Z !1/2 Z !1/2
2
2 00
u (y ) dy
≤ h 1 dy
Ii Ii
Z !1/2
≤ h3/2 u 00 (y )2 dy
I
i
= h3/2
u 00 (x )
L2 (I ) .
i
And
Z
ku − uh k2L2 (I ) = (u − uh )2 dx
i
Ii
Z
≤ ku − uh k2C (I ) dx
i
Ii
2
≤ h ku − uh kC (I ) .
i
Therefore,
and
Therefore
ku − uh kL2 (Ii ) ≤ 2h2
u 00 (x )
L2 (I ) .
i
Hence
ku − uh kL2 (I ) ≤ 2h2
u 00 (x )
L2 (I ) .
2. For the linear basis we have the Element solution uh (xi ) = u (xi ) and uh (xi +1 ) = u (xi +1 ) on element
Ii = [xi , xi +1 ]. and
u h ( xi + 1 ) − u h ( xi ) u ( xi + 1 ) − u ( xi ) 1 xi + 1 0
Z Z
0 1
uh ( x ) = = = u (y )dy = u 0 (y )dy,
h h h xi h Ii
1 xi + 1 0
Z Z
h 1
u 0 (x ) = u 0 (x ) = u (x )dy = u 0 (x )dy.
h h xi h Ii
Therefore
Z
1
uh0 (x ) − u 0 (x ) = u 0 (y ) − u 0 (x )dy
h Ii
Z Z y
1
= u 00 (ξ )dξdy.
h Ii x
so
2 Z 2
u 0 − uh0
2 = uh0 (x ) − u 0 (x ) dx
L (Ii )
Ii
Z Z Z y !2
1
= u 00 (ξ )dξdy dx
h2 Ii Ii x
Z Z Z !2
1
≤ u 00 (ξ )dξdy dx
h2 Ii Ii Ii
Z Z !2 Z
1
= u 00 (ξ )dξdy dx
h2 Ii Ii Ii
Z Z !2
1 00
= dy u (ξ )dξ
h Ii Ii
Z !2
1 00
= h u (ξ )dξ
h Ii
Z !2
00
= h u (ξ )dξ
Ii
Z !1/2 Z !1/2 2
2
2 00
≤ h 1 dξ u (ξ ) dξ
Ii Ii
Z !
≤ h2 u 00 (ξ )2 dξ
Ii
hence
u 0 − uh0
≤ Ch
u 00
2 .
L2 ( I i ) L (I ) i
Therefore,
u 0 − uh0
≤ Ch
u 00
L2 (I ) .
L2 ( I )
10.2 Problems
Problem 10.1. (Prelim Jan. 2008#8) Let Ω ⊂ R2 be a bounded domain with a smooth boundary. Consider
a 2-D poisson-like equation
−∆u + 3u = x2 y 2 , in Ω,
= 0, on ∂Ω.
u
ku − uh kH 1 ≤ C inf ku − vh kH 1 ,
vh ∈Vh
with C independent of h, where Vh denotes a finite element subspace of H 1 (Ω) consisting of contin-
uous piecewise polynomials of degree k ≥ 1.
Solution. 1. For this pure Dirichlet Problem, the test functional space v ∈ H01 . Multiple the test func-
tion on the both sides of the original function and integral on Ω, we get
Z Z Z
− ∆uvdx + uvdx = xyvdx.
Ω Ω Ω
Let
Z Z Z
a(u, v ) = ∇u∇vdx + uvdx, f (v ) = xyvdx.
Ω Ω Ω
Then, the
(a) Ritz variational problem is: find uh ∈ H01 , such that
1
J (uh ) = min a(uh , uh ) − f (uh ).
2
(b) Galerkin variational problem is: find uh ∈ H01 , such that
a(uh , uh ) = f (uh ).
(b)
Z Z
2
a(u, u ) = (∇u ) dx + u 2 dx
Ω Ω
So,
Z Z
2
a(u, u ) = |∇u| dx + |u|2 dx
Ω Ω
= k∇uk2L2 (Ω) + kuk2L2 (Ω)
= kuk2H 1 (Ω) .
(c)
Z
f (v ) ≤ xyv dx
Ω
Z
≤ max |xy| |v| dx
Ω
Z !1/2 Z !1/2
2 2
≤ C 1 dx |v| dx
Ω Ω
≤ C kvkL2 (Ω) ≤ C kvkH 1 (Ω) .
by Lax-Milgram theorem, we get that e Galerkin method has a unique solution uh . Moreover,
a(vh , vh ) = f (vh ).
a(u, vh ) = f (vh ).
a(u − uh , vh ) = 0.
Then by coercivity
ku − uh k2H 1 (Ω) ≤ a(u − uh , u − uh )
= a(u − uh , u − vh ) + a(u − uh , vh − uh )
= a(u − uh , u − vh )
≤ ku − uh kH 1 (Ω) ku − vh kH 1 (Ω) .
Therefore,
ku − uh kH 1 ≤ C inf ku − vh kH 1 ,
vh ∈Vh
n o
Problem 10.2. (Prelim Aug. 2006#9) Let Ω := (x, y ) : x2 + y 2 < 1 , consider the poisson problem
−∆u + 2u = xy, in Ω,
= 0, on ∂Ω.
u
ku − uN kE ≤ C inf ku − vN kE ,
vN ∈VN
where
Z
kvkE := |∇v|2 + 2v 2 dxdy.
Ω
Solution. J
References
[1] S. C. Brenner and R. Scott, The mathematical theory of finite element methods, vol. 15, Springer, 2008.
108
[2] A. Iserles, A First Course in the Numerical Analysis of Differential Equations (Cambridge Texts in Applied
Mathematics), Cambridge University Press, 2008. 80, 84, 86, 92, 151
[3] Y. Saad, Iterative methods for sparse linear systems, Siam, 2003. 2, 39
[4] A. J. Salgado, Numerical math lecture notes: 571-572. UTK, 2013-14. 1
[5] S. M. Wise, Numerical math lecture notes: 571-572. UTK, 2012-13. 1
Appendices
and then λ = λ, so λ is real. Moreover, since ρ (A) ⊂ (0, ∞), then we have λ is positive.
Problem A.2. (Sample#2) Suppose dim(A) = n. If A has n distinct eigenvalues , then A is diagonalizable
.
Solution. (Sketch) Suppose n = 2, and let λ1 , λ2 be distinct eigenvalues of A with corresponding eigen-
vectors v1 , v2 . Now, we will use contradiction to show v1 , v2 are lineally independent. Suppose v1 , v2 are
lineally dependent, then
c1 v1 + c2 v2 = 0, (210)
c1 λ1 v1 + c2 λ1 v2 = 0. (212)
Since λ1 , λ2 and v2 , 0, then c2 = 0. Similarly, we can get c1 = 0. Hence, we get the contradiction.
A similar argument gives the result for n. Then we get A has n linearly independent eigenvectors . J
i.e.
1 + α + αv ∗ u = 0,
i.e.
1
α=− , 1 , −v ∗ u.
1 + v∗u
Aw = (In + uv ∗ )w = w + uv ∗ w = 0
Aw = (In + uv ∗ )βu = β (u + uv ∗ u ) = β (u + (v ∗ u )u ) = 0.
hence span(u ) ∈ w.
J
for any x ∈ Rn .
3. Let b ∈ Rn be given. Prove that x∗ ∈ Rn solves Ax = b if and only if x∗ minimizes the quadratic
function f : Rn → R defined by
1 T
f (x ) = x Ax − xT b.
2
√ √
Solution. √ 1. (a) Obviously, kxkA = xT Ax ≥ 0. When x = 0, then kxkA = xT Ax = 0; when kxkA =
xT Ax = 0, then we have (Ax, x ) = 0, since A is SPD, therefore, x ≡ 0.
√ √ √
(b) kλxkA = λxT Aλx = λ2 xT Ax = |λ| xT Ax = |λ| kxkA .
(c) Next we will show
x + y
A = kxkA +
y
A . First, we would like to show
y T Ax ≤ kxkA
y
.
A
Then
c.s.
y T Ax = y T RT Rx = (Ry )T Rx = (Rx, Ry ) ≤ kRxk2
Ry
= kxkA
y
.
2 A
And
2
x + y
= (x + y, x + y )A = (x, x )A + 2(x, y )A + (y, y )A
A
≤ kxkA + 2 y T Ax +
y
A
≤ kxkA + 2 kxkA
y
A +
y
A
2
= kxkA +
y
A .
therefore
x + y
= kxkA +
y
.
A A
i.e.
p p
λ1 kxk2 ≤ kRxk2 = kxkA ≤ λn kxk2 .
3. Since
∂ T ∂ T ∂
x Ax = x Ax + xT (Ax )
∂xi ∂xi ∂xi
0
..
.
0
T
= [0, · · · , 0, 1, 0, · · · , 0]Ax + x A 1 i
i
0
..
.
0
= (Ax )i + AT x = 2 (Ax )i .
i
and
∂ T ∂ T
x b = x b = [0, · · · , 0, 1, 0, · · · , 0]b = bi .
∂xi ∂xi i
Therefore,
1
∇f (x ) =
2Ax − b = Ax − b.
2
If Ax∗ = b, then ∇f (x∗ ) = Ax∗ − b = 0, therefore x∗ minimizes the quadratic function f. Conversely,
when x∗ minimizes the quadratic function f, then ∇f (x∗ ) = Ax∗ − b = 0, therefore Ax∗ = b.
J
Problem A.5. (Sample#9) Suppose that the spectrum of A ∈ Rn×n sym is denoted ρ (A) = {λ1 , λ2 , · · · , λn } ⊂ R.
Let S = {x1 , · · · , xn } be the orthonormal basis of eigenvectors of A, with Axk = λk xk , for k = 1, · · · , n. The
Rayleigh quotient of x ∈ Rn∗ is defined as
xT Ax
R(x) : = .
xT x
Prove the following facts:
1.
Pn 2
j = 1 λj α j
R(x) : = Pn 2
j =1 αj
where αj = xT xj .
2.
1. First, we need to show that x = nj=1 αj xj is the unique representation of x w.r.t. the or-
P
Solution.
thonormal basis S. Since S = {x1 , · · · , xn } is the orthonormal basis of eigenvectors of A, then nj=1 xT xj xj
P
2. Since,
Pn 2
xT Ax j =1 λj αj
R(x) := T = Pn 2
.
x x j = 1 αj
, then
Pn 2 Pn 2
j = 1 αj j = 1 αj
min λj Pn 2
≤ R(x) ≤ max λj Pn 2
,
j = 1 αj j =1 αj
j j
i.e.
Hence
Problem A.6. (Sample #31) Let A ∈ Rn×n be symmetric positive define (SPD). Let b ∈ Rn . Consider
solving Ax = b using the iterative method
Mxn+1 = N xn + b, n = 0, 1, 2, · · ·
xn+1 = M −1 N xn + M −1 b.
Let G = M −1 N = M −1 (M −A) = I −M −1 A, If we can prove that ρ (G ) < 1, then this method converges.
Let λ be any eigenvalue of G and x be the corresponding eigenvector, i.e.
Gx = λx.
then
(I − M −1 A)x = λx,
i.e.
(M − A)x = λMx,
i.e.
(1 − λ)Mx = Ax.
(1 − λ)x∗ Mx = x∗ Ax.
So we have
1 ∗
x∗ Mx = x Ax.
1−λ
taking conjugate transpose of which yields
1 1
x∗ M ∗ x = x∗ A∗ x = x∗ Ax.
1−λ 1−λ
Then, we have
!
1 1
x ∗ (M + M ∗ − A)x = + − 1 x∗ Ax
1−λ 1−λ
!
λ 1
= + x∗ Ax
1−λ 1−λ
1 − λ2 ∗
= x Ax.
|1 − λ|2
1 − λ2 > 0.
i.e.
|λ| < 1.
2.
T
where A = L + D + U . Since A is SPD, then MGS + MGS − A = D + L + D T + LT − A = D + LT − U is
SPD. Therefore, From the part 1, we get that the Gauss-Seidel Method converges.
J
Problem A.7. (Sample #32) Let A ∈ Rn×n be symmetric positive define (SPD). Let b ∈ Rn . Consider
solving Ax = b using the iterative method
Mxn+1 = N xn + b, n = 0, 1, 2, · · ·
Solution. Let ek = xk − x. And rewrite the scheme to the canonical form B = M, α = 1, then
x k +1 − x k
!
B + Axk = b = Ax.
α
so, we get
e k +1 − e k
!
B + Aek = 0.
α
Let v k +1 = ek +1 − ek , then
1 k +1
Bv + Aek = 0.
α
Taking the conjugate transport of the above equation, then we get
1 ∗ k +1 1
Bv + A ∗ ek = B∗ v k +1 + Aek = 0.
α α
therefore
1 B + B∗ k + 1
( )v + Aek = 0.
α 2
B + B∗
Let Bs = 2 . Then take the inner product of both sides with v k +1 ,
1
(B v k +1 , v k +1 ) + (Aek , v k +1 ) = 0.
α s
Since
1 k +1 1 1 1
ek = (e + e k ) − (e k +1 − e k ) = (e k +1 + e k ) − v k +1 .
2 2 2 2
Therefore,
1
0 = (B v k +1 , v k +1 ) + (Aek , v k +1 )
α s
1 1 1
= (B v k +1 , v k +1 ) + (A(ek +1 + ek ), v k +1 ) − (Av k +1 , v k +1 )
α s 2 2
1 α 1
= ((Bs − A)v k +1 , v k +1 ) + (A(ek +1 + ek ), v k +1 )
α 2 2
1 α 1
2
2
= ((Bs − A)v k +1 , v k +1 ) + (
ek +1
A −
ek
A )
α 2 2
M +M T −A
By assumption, Q = Bs − α2 A = 2 > 0, i.e. there exists m > 0, s.t.
2
(Qy, y ) ≥ m
y
2 .
Therefore,
m
k +1
2 1
k +1
2
k
2
v
+ (
e
−
e
) ≤ 0.
α 2 2 A A
i.e.
2m
k +1
2
k +1
2
k
2
v
+
e
≤
e
.
α 2 A A
Hence
2
2
ek +1
≤
ek
.
A A
and
2
ek +1
→ 0.
A
Problem A.8. (Sample #33) Consider a linear system Ax = b with A ∈ Rn×n . Richardson’s method is an
iteration method
Mxk +1 = N xk + b
xk +1 = (I − wA)xk + wb.
So the error transfer operator is T = I − wA. Then if λi is the eigenvalue of A, then 1 − wλi should be
the eigenvalue of T . The sufficient and necessary condition of convergence is ρ (T ) < 1, i.e.
|1 − wλi | < 1
2
w< .
λn
2
conversely, if w < λn , then ρ (T ) < 1, then the scheme converges.
2. The minimum is attachment at |1 − ωλn | = |1 − ωλ1 |(Figure.1), i.e.
ωλn − 1 = 1 − ωλ1 .
Therefore, we get
2
ωopt = .
λ1 + λn
J
Sn := I + A + A2 + · · · + An .
lim Sn = (I − A)−1 .
n→∞
Moreover,
Ak x
kAk
Ak−1 x
k
A
= sup ≤ sup ≤ · · · ≤ kAkk .
0,x∈Cn kxk 0,x∈C n kxk
From the properties of geometric series, Sn converges if only if |A| < 1. Therefore, we get if |A| < 1
then A is convergent. Conversely, if A is convergent, then |A| < 1. Hence Sn converges.
2.
(I − A)x
= kx − Axk ≥ kxk − kAxk ≥ kxk − kAk kxk = (1 − kAk) kxk .
If A is convergent, then kAk , 0. Therefore, if
(I − A)x
= 0, then kxk = 0, i.e. ker (I − A) = 0. Hence,
I − A is non-singular. From the definition of Sn , we get
n
X n
X +1
( I − A ) Sn = k
A − Ak = A0 − An+1 = I − An+1 .
k =0 k =1
Taking limit on both sides of the above equation with the fact |A| < 1, then we get
(I − A) lim Sn = I.
n→∞
lim Sn = (I − A)−1 .
n→∞
J
Problem A.10. (Sample #40) Show that if λ is an eigenvalue of A∗ A, where A ∈ Cn×n , then
0 ≤ λ ≤ kAk kA∗ k
Problem A.11. (Sample #41) Suppose A ∈ Cn×n and A is invertible. Prove that
r
λn
κ2 ≤ .
λ1
√
Solution. Since κ2 = kAk2
A−1
2 and kAk2 = max ρ (A) = λn . therefore
√
λ
κ2 = kAk2
A−1
2 = √ n .
λ1
J
Problem A.12. (Sample #34) Let A = [ai,j ] ∈ Cn×n be invertible and b ∈ Cn . Prove that the classical Jacobi
iteration method for approximating the solution to Ax = b is convergent, for any starting value x0 , if A is
strictly diagonally dominant, i.e.
X
ai,i < ai,k , ∀ i = 1, · · · , n.
k,i
So,
X X aij
kT k∞ = max |tij | = max | |.
i i aii
j i,j
Therefore,
X |aij |
1≥ .
|aii |
j,i
Hence, kT k∞ < 1 J
Problem A.13. (Sample #35) Let A = [ai,j ] ∈ Cn×n be invertible and b ∈ Cn . Prove that the classical
Gauss-Seidel iteration method for approximating the solution to Ax = b is convergent, for any starting
value x0 , if A is strictly diagonally dominant, i.e.
X
ai,i < ai,k , ∀ i = 1, · · · , n.
k,i
(D + L)(xk +1 − xk ) + Axk = b.
We want to show If A is diagonal dominant , then kTGS k < 1, then Jacobi Method convergences. From the
definition of T, we know that T for Gauss-Seidel iteration Method is as follows
which implies
P
j>i |aij |
( )
γ = maxi P | ≤ 1.
|aii | − j<i |aij
Moreover
X X
X
X
|((L + D )y )i0 | = | ai0 j yj + ai0 i0 yj | ≥ |ai0 i0 yj | − | ai0 j yj | = |ai0 i0 |
y
∞ − | ai0 j yj | ≥ |ai0 i0 |
y
∞ − |ai0 j |
y
∞ .
j<i0 j<i0 j<i0 j<i0
which implies
P
j>i0 |ai0 j |
y
≤
∞
P kxk∞ .
|ai0 i0 | − j<i0 |ai0 j |
So,
which implies
kTGS k∞ ≤ γ < 1.
Problem A.14. (Sample #38) Let A ∈ Cn×n be invertible and suppose b ∈ Cn∗ satisfies Ax = b. let the
perturbations δx, δb ∈ Cn satisfy Aδx = δb, so that A(x + δx ) = b + δb.
1. Prove the error (or perturbation) estimate
2. Show that for any invertible matrix A, the upper bound for kδbk
kbk
above can be attained for suitable
choice of b and δb. (In other words, the upper bound is sharp.)
Therefore
kδbk 1 1
≤ kδxk ,
≤ .
kAk −1
A
kbk kxk
Hence
1 kδbk kδxk
≤ .
κ (A) kbk kxk
Therefore
1 kAk
≤ .
kxk kbk
Hence,
kδxk kδbk
≤ κ (A) .
kxk kbk
So,
1 kδbk kδxk kδbk
≤ ≤ κ (A) .
κ (A) kbk kxk kbk
where δx = x̂ − x.
Solution. Since
A−1
kδAk < 1, then we have
A−1 δA
≤
A−1
kδAk < 1.
Therefore,
(I − A−1 δA)−1
≤ 1
.
1 − −1 δA
A
(I + A−1 δA)−1
≤ 1
.
1 −
A−1 δA
δx = x + δx − x
= (A + δA)−1 (b + δb ) − A−1 b
= (A + δA)−1 AA−1 (b + δb ) − A−1 b
−1
= (A + δA)−1 A−1 A−1 (b + δb ) − A−1 b
−1
= (A−1 A + A−1 δA) A−1 (b + δb ) − A−1 b
−1
= (I + A−1 δA) A−1 (b + δb ) − A−1 b
−1
= (I + A−1 δA) A−1 (b + δb ) − (I + A−1 δA)A−1 b
−1
= (I + A−1 δA) A−1 δb − A−1 δAA−1 b
Therefore,
1
kδxk ≤
A−1 δb
+
A−1 δAA−1 b
1 −
A δA
−1
1
≤
A−1
kδbk +
A−1
kδAk
A−1 b
1 −
A−1 δA
1
=
A−1
kδbk +
A−1
kδAk kxk
1 −
A δA
−1
!
κ (A) kδbk kδAk kxk
=
+
1 −
A−1 δA
kAk kAk
where δx = x̂ − x.
Solution. Since
A−1
kδAk < 1, then we have
A−1 δA
≤
A−1
kδAk < 1.
Therefore,
(I − A−1 δA)−1
≤ 1
.
1 −
A−1 δA
(I + A−1 δA)−1
≤ 1
.
1 −
A−1 δA
δx = x + δx − x
= (A + δA)−1 b − A−1 b
= (A + δA)−1 AA−1 b − A−1 b
−1
= (A + δA)−1 A−1 A−1 b − A−1 b
−1
= (A−1 A + A−1 δA) A−1 b − A−1 b
−1
= (I + A−1 δA) A−1 b − A−1 b
−1
= (I + A−1 δA) A−1 b − (I + A−1 δA)A−1 b
−1
= (I + A−1 δA) −A−1 δAA−1 b
Therefore,
1
kδxk ≤
A−1 δAA−1 b
1 −
A−1 δA
1
≤
A−1
kδAk
A−1 b
1 −
A−1 δA
1
=
A−1
kδAk kxk
1 −
A−1 δA
!
κ (A) kδAk kxk
=
1 −
A−1 δA
kAk
Problem A.17. (Sample #40) Show that if λ is an eigenvalue of A∗ A, where A ∈ Cn×n , then
0 ≤ λ ≤ kA∗ k kAk .
Problem A.19. (Sample #42) Suppose A ∈ Cn×n and A is invertible. Prove that
q
κ2 ≤ κ1 ( A ) κ∞ ( A ) .
Solution.
Claim A.2.
kAk22 ≤ kAk1 kAk∞ .
Proof.
kAk22 = ρ (A)2 = λ ≤ kAk1 kA∗ k1 = kAk1 kAk∞ .
where λ is the eigenvalue of A∗ A. J
Since κ2 = kAk2
A−1
2 , κ1 = kAk1
A−1
1 and κ∞ = kAk∞
A−1
∞ .
q
q q
−1
≤
p
A−1
A−1
≤
−1
−1
= κ1 ( A ) κ∞ ( A ) .
kAk2
A
2
kAk1 kAk∞ 1 ∞
kAk1
A
1
kAk∞
A
∞
J
Problem A.20. (Sample #44) Suppose A, B ∈ Rn×n and A is non-singular and B is singular. Prove that
1 kA − Bk
≤ ,
κ (A) kAk
where κ (A) = kAk ·
A−1
, and k·k is an reduced matrix norm.
Solution. Since B is singular, then there exists a vector x , 0, s.t. Bx = 0. Since A is non-singular, then
A−1 is also non-singular. Moreover, A−1 Bx = 0. Then, we have
x = x − A−1 Bx = (I − A−1 B)x.
So
kxk =
(I − A−1 B)x
≤
A−1 A − A−1 B
kxk ≤
A−1
kA − Bk kxk .
Since x , 0, so
1 ≤
A−1
kA − Bk .
1 kA − Bk
≤ ,
−1
A
kAk kAk
i.e.
1 kA − Bk
≤ .
κ (A) kAk
J
Problem A.21. (Sample #1) Let {xn } be a sequence generated by Newton’s method. Suppose that the initial
guess x0 is well chosen so that this sequence converges to the exact solution x∗ . Prove that if f (x∗ ) =
f 0 (x∗ ) = · · · = f m−1 (x∗ ) = 0, f m (x∗ ) , 0, xn converges linearly to x∗ with
e k +1 m−1
lim k
= .
k→∞ e m
f (x k )
x k +1 = x k − .
f 0 (x k )
Let ek = xk − x∗ , then
e k +1 = x k + 1 − x∗
f (x k )
= x k − 0 k − x∗
f (x )
f (x k )
= ek − .
f 0 (x k )
Therefore,
e k +1 f (x k )
= 1 − .
ek e k f 0 (x k )
Since x0 is well chosen so that this sequence converges to the exact solution x∗ , therefore we have the Taylor
expansion for f (xk ), f 0 (xk ) at x∗ , i.e.
Solution. J
Problem A.23. (Sample #3) Let a ∈ Rn and R > 0 be given. Suppose that f : B(a, R) → Rn , fi ∈
C 2 (B(a, R)), for each i = 1, · · · , n. Suppose that there is
a point ξ
∈ B(a, R), such that f(ξ ) = 0, and
that the Jacobian matrix Jf (x) is invertible, with estimate
[Jf (x)]−1
2 ≤ β, for any x ∈ B(a, R). Prove that
the sequence {xk } defined by Newton’s method,
Solution. J
Problem A.24. (Sample #6) Assume that f : R → R, f ∈ C 2 (R), f 0 (x ) > 0 for all x ∈ R, and f 00 (x ) > 0,
for all x ∈ R.
1. Suppose that a root ξ ∈ R exists. Prove that it is unique. Exhibit a function satisfying the assump-
tions above that has no root.
2. Prove that for any starting guess x0 ∈ R, Newton’s method converges, and the convergence rate is
quadratic.
Solution. 1. Let x1 and x2 are the two different roots. So, f (x1 ) = f (x2 ) = 0, then by Mean value
theorem, we have that there exists η ∈ [x1 , x2 ], such f (η ) = 0 which contradicts f 0 (x ) > 0.
2. example f (x ) = ex .
3. Let x∗ be the root of f (x ). From the Taylor expansion, we know
1
0 = f (x∗ ) = f (xk ) + f 0 (xk )(x∗ − xk ) + f 00 (θ )(x∗ − xk )2 ,
2
where θ is between x∗ and xk . Define ek = x∗ − xk , then
1
0 = f (x∗ ) = f (xk ) + f 0 (xk )(ek ) + f 00 (θ )(ek )2 .
2
so
h i−1 1h i−1
f 0 (xk ) f (xk ) = −(ek ) − f 0 (xk ) f 00 (θ )(ek )2 .
2
From the Newton’s scheme, we have
h i−1
xk +1 = xk − f 0 (xk ) f (xk )
x∗ = x∗
So,
h i−1 1h i−1
ek +1 = ek + f 0 (xk ) f (xk ) = − f 0 (xk ) f 00 (θ )(ek )2 ,
2
i.e.
f 00 (θ )
e k +1 = − h i (e k )2 ,
0
2 f (x )k
Therefore,
k + 1
f 00 (θ ) k 2 C1 k 2
e ≤ h i (e ) ≤ e .
2 f 0 (xk ) 2C2
This implies
2
xk +1 − x∗ ≤ C xk − x∗ .
f ( xk ) f (y )
y k = xk − 0
, xk + 1 = y k − 0 k
f ( xk ) f ( xk )
xk + 1 − x ∗ f 00 (xk )
lim = ,
k→∞ (yk − x∗ )(xk − x∗ ) f 0 ( xk )
xk +1 − x∗ 1 f 00 (xk )
!
lim = .
k→∞ (xk − x∗ )3 2 f 0 ( xk )
3. Would you say that this method is faster than Newton’s method given that its convergence is cubic?
Solution. 1. First, we will show that if xk ∈ [x − h, x + h], then yk ∈ [x − h, x + h]. By Taylor expansion
formula, we have
1 00
0 = f (x∗ ) = f (xk ) + f 0 (xk )(x∗ − xk ) + f (ξk )(x∗ − xk )2 ,
2!
where ξ is between x and xk . Therefore, we have
1 00
f (xk ) = −f 0 (xk )(x∗ − xk ) − f (ξk )(x∗ − xk )2 .
2!
Plugging the above equation to the first step of the Newton’s method, we have
1 f 00 (ξk ) ∗
y k = xk + ( x ∗ − xk ) + ( x − xk ) 2 .
2! f 0 (xk )
then
1 f 00 (ξk ) ∗
yk − x ∗ = ( x − xk ) 2 . (214)
2! f 0 (xk )
Therefore,
1 f 00 (ξk )
1 f 00 (ξk ) ∗
yk − x∗ = ∗ 2 (x − xk ) (x∗ − xk ) .
0
( x − x )
k ≤ 0
2! f (xk ) 2 f ( xk )
Since we can choose the initial value very close to x∗ , shah that
00
f (ξ ) (x∗ − x ) ≤ 1
f 0 (x ) k
k
By mean value theorem, we have there exists ηk between x∗ and xk , such that
(yk − x∗ )2 00
f ( yk ) = f (x∗ ) + f 0 (x∗ )(yk − x∗ ) + f ( γk )
2
(y − x∗ )2 00
= f 0 (x∗ )(yk − x∗ ) + k f ( γk ) ,
2
where γ is between yk and x∗ . Plugging the above two equations to the second step of the Newton’s
method, we get
(yk − x∗ )2 00
" #
1
xk + 1 − x ∗ = f 00
η (
k k x − x ∗
)( y k − x ∗
) − f 0 ∗
( x )( y k − x ∗
) − f ( γ k ) + ( y k − x ∗ 0 ∗
) f ( x )
f 0 ( xk ) 2
(yk − x∗ )2 00
" #
1 00 ∗ ∗
= f η k ( x k − x )( y k − x ) − f ( γk ) . (215)
f 0 ( xk ) 2
xk + 1 − x ∗ f 00 ηk (yk − x∗ )f 00 (γk )
= − .
(xk − x∗ )(yk − x∗ ) f 0 (xk ) 2(xk − x∗ )f 0 (xk )
xk + 1 − x ∗ f 00 ηk 1 f 00 (ξk ) ∗ f 00 (γk )
= − ( x − x k ) .
(xk − x∗ )(yk − x∗ ) f 0 (xk ) 4 f 0 (xk ) f 0 ( xk )
xk + 1 − x ∗ f 00 (x∗ )
lim = .
k→∞ (xk − x∗ )(yk − x∗ ) f 0 (x ∗ )
1 2 f 0 ( xk )
= .
yk − x∗ (x∗ − xk )2 f 00 (ξk )
Hence
!2
x − x∗ 1 f 00 (x∗ )
lim k +1 ∗ 3 = .
k→∞ (xk − x ) 2 f 0 (x ∗ )
J
Problem A.26. (Sample #1) Show that, if z is a non-zero complex number that iix on the boundary of the
linear stability domain of the two-step BDF method
4 1 2
yn+2 − yn+1 + yn = hf (xn+2 , yn+2 ),
3 3 3
then the real part of z must be positive. Thus deduce that this method is A-stable.
• coercive:
α kuk2V ≤ a(u, v ) , ∃α > 0, ∀u ∈ V ,
So, we have
λ
kuh kV ≤ .
α
a(u, v ) = L(v )
a(uh , v ) = L(v )
then we have the so called Galerkin Orthogonal a(u − uh , v ) = 0 for all v ∈ V . Then by coercivity
α ku − uh k2V ≤ a(u − uh , u − uh )
= a(u − uh , u − w ) + a(u − uh , w − uh )
= a(u − uh , u − w )
≤ γ ku − uh kV ku − wkV .
therefore
γ
ku − uh kV ≤ ku − wkV .
α
Hence
γ
ku − uh kV ≤ inf ku − wkV .
α w∈Sh
J
1 n µ as
ujn+1 =
uj−1 + ujn+1 − ujn+1 − uj−1
n
, µ= ,
2 2 h
for approximating solutions to the Cauchy problem for the advection equation
∂u ∂u
+ =0
∂t ∂x
where a > 0. Here h > 0 is the space step size, and s > 0 us the time step size.
1. Prove that, if s = C1 h, where C1 is fixed positive constant, then the local truncation error satisfies
the estimate
Tln ≤ C0 (s + h)
, where C0 > 0 us a constant independent of s and h.
2. Use the von Neumann analysis to show that the Lax-Friedrichs scheme is stable provided the CFL
as
condition
0 < µ = h ≤ 1 holds. In other words, compute the amplification factor, g (ξ ), and show
that g (ξ ) ≤ 1, for all values of ξ, provided µ ≤ 1.
Solution. 1. Then the Lax-Friedrichs method for solving the above partial differential equation is given
by:
ujn+1 − 12 (ujn+1 + uj−1
n
) n
ujn+1 − uj−1
+a =0
∆t 2 ∆x
Or, rewriting this to solve for the unknown ujn+1 ,
1 n ∆t
ujn+1 = n
(uj +1 + uj−1 )−a n
(u n − uj−1 )
2 2 ∆x j +1
i.e.
1 n µ as
ujn+1 =
uj−1 + ujn+1 − ujn+1 − uj−1
n
, µ= .
2 2 h
Let ū be the exact solution and ūjn = ū (n∆t, j∆x ). Then from Taylor Expansion, we have
∂ n 1 ∂2
ūjn+1 = ūjn + ∆t ūj + (∆t )2 2 ū (ξ1 , j∆x ), tn ≤ ξ1 ≤ tn+1 ,
∂t 2 ∂t
∂ 1 ∂2
n
ūj−1 = ūjn − ∆x ūjn + (∆x )2 2 ū (n∆t, ξ2 ), xj−1 ≤ ξ2 ≤ xj ,
∂x 2 ∂x
∂ 1 ∂2
ūjn+1 = ūjn + ∆x ūjn + (∆x )2 2 ū (n∆t, ξ3 ), xj ≤ ξ3 ≤ xj + 1 .
∂x 2 ∂x
Then the truncation error T n+1 of this scheme is
uin+1 − 21 (uin+1 + ui−1
n
) u n
− u n
T n+1 = + a i +1 i−1
∆t 2 ∆x
O (s )2 + O (h)2
= C
s
If s = C1 h, where C1 is fixed positive constant, then the local truncation error
T n+1 = C0 (s + h).
then we have
1 n µ
g (ξ )ujn = uj−1 + ujn+1 − ujn+1 − uj−1
n
,
2 2
and
1 i (j−1)∆xξ µ
g (ξ )eij∆xξ = e + ei (j +1)∆xξ − ei (j +1)∆xξ − ei (j−1)∆xξ .
2 2
Then we have
1 −i∆xξ µ
g (ξ ) = e + ei∆xξ − ei∆xξ − e−i∆xξ
2 2
= cos(∆xξ ) + iµ sin(∆xξ ).
From von Neumann analysis, we know that the Lax-Friedrichs scheme is stable if g (ξ ) ≤ 1, i.e.
i.e.
µ ≤ 1.
u n+1
≤ ku n k∞ ,
∞
Solution. This problem is similar to Sample #14. The scheme can be rewritten as
µ n+1 s n+1 µ n+1 µ n s µ
(1 + µ)ujn+1 = uj−1 − uj + uj +1 + uj−1 + (1 − µ − )ujn + ujn+1 .
2 2 2 2 2 2
Then we have
µ s µ µ s n µ n
1 + µ ujn+1 ≤ uj−1
n+1 n+1 n+1
n
+ uj + uj +1 + uj−1 + (1 − µ − ) uj + uj +1 .
2 2 2 2 2 2
Therefore
µ
s
µ
µ s µ
(1 + µ)
u n+1
∞ ≤
u n+1
∞ +
u n+1
∞ +
u n+1
∞ + ku n k∞ + (1 − µ − ) ku n k∞ + ku n k∞ .
2 2 2 2 2 2
if 0 < µ + 2s ≤ 1, then
µ
s
µ
µ s µ
(1 + µ)
u n+1
∞ ≤
u n+1
∞ +
u n+1
∞ +
u n+1
∞ + ku n k∞ + (1 − µ − ) ku n k∞ + ku n k∞ .
2 2 2 2 2 2
i.e.
s
s
(1 − )
u n+1
∞ ≤ (1 − ) ku n k∞
2 2
Hence
u n+1
≤ ku n k∞ .
∞
Problem A.30. (Sample #8) 1D Discrete Poincaré inequality: Let Ω = (0, 1) and Ωh be a uniform grid of
size h. If Y ∈ Uh is a mesh function on Ωh such that Y (0) = 0, then there is a constant C, independent of Y
and h, for which
kY k2,h ≤ C
δ̄Y
2,h .
Solution. I consider the following uniform partition (Figure. A1) of the interval (0, 1) with N points.
x1 = 0 x2 xN −1 xN = 1
Then,
X N
Yi−1 − Yi = |YN |.
i =2
and
N N X N
1/2 N 1/2
X X
Y − Y X Y − Y 2
h i−1 i i−1 i
|Yi−1 − Yi | = h2 .
|YN | ≤ ≤
h h
i =2 i =2 i =2 i =2
Therefore
K K
X X Yi−1 − Yi 2
2 2
|YK | ≤ h
h
i =2 i =2
K 2
Yi−1 − Yi .
X
= h2 (K − 1) h
i =2
1. When K = 2,
2
2 2 Y1 − Y2
|Y2 | ≤ h
h .
2. When K = 3,
Y − Y2 2 Y2 − Y3 2
!
|Y3 |2 ≤ 2h2 1 + .
h h
3. When K = N ,
Y − Y2 2 Y2 − Y3 2 − YN 2
!
Y
|YN |2 ≤ (N − 1)h2 1 + + · · · + N −1 .
h h h
Since Y1 = 0, so
N N
N (N − 1) 2 X Yi−1 − Yi 2
X
2 .
|Yi | ≤ h
2 h
i =1 i =2
And then
N N N
Yi−1 − Yi 2 Yi−1 − Yi 2
! X
1 X
2 N 2
X 1 1 2
|Yi | ≤ h = 2 + 2(N − 1) h .
(N − 1)2
2(N − 1) h h
i =1 i =2 i =2
1
Since h = N −1 , so
N N
Yi−1 − Yi 2
! X
2
X
2 1 1 2
h |Yi | ≤ + h
.
2 2(N − 1) h
i =1 i =2
then
N N
Yi−1 − Yi 2
! X
X
2 1 1
h |Yi | ≤ + h
.
2 2(N − 1) h
i =1 i =2
i.e,
!
1 1
2
kY k22,h ≤ +
δ̄Y
.
2,h
2 2(N − 1)
since N ≥ 2, so
2
kY k22,h ≤
δ̄Y
2,h .
Hence,
kY k2,h ≤ C
δ̄Y
2,h .
J
Problem A.31. (Sample #12) Discrete maximum principle: Let A = tridiag{ai , bi , ci }ni=1 ∈ Rn×n be a
tridiagional matrix with the properties that
bi > 0, ai , ci ≤ 0, ai + bi + ci = 0.
Prove the following maximum principle: If u ∈ Rn is such that (Au )i =2,··· ,n−1 ≤ 0, then ui ≤ max{u1 , un }.
ak uk−1 + bk uk + ck uk +1 < 0.
This is contradiction to (Au )i =2,··· ,n−1 < 0. Therefore, If u ∈ Rn is such that (Au )i =2,··· ,n−1 < 0, then
ui ≤ max{u1 , un }.
2. For (Au )i =2,··· ,n−1 = 0:
Since (Au )i =2,··· ,n−1 = 0, so
ak uk−1 + bk uk + ck uk +1 = 0.
Since ak + ck = −bk , so
And ak < 0, ck < 0, uk−1 − uk ≤ 0, uk +1 − uk ≤ 0, so uk−1 = uk = uk +1 , that is to say, uk−1 and uk +1 is also
the maximum points. Bu using the same argument again, we get uk−2 = uk−1 = uk = uk +1 = uk +2 .
Repeating the process, we get
u1 = u2 = · · · = un−1 = un .
kAxk2 ≤ kxk2 ,
kAxk∞ ≤ kxk∞ ,
for any x ∈ Rm , provided µ ≤ 1.(In other words, the scheme may only be conditionally stable in the
max norm.)
n+1 n
u1 u1
n+1 u n
u2 2
un+1 = . and un = . .
.. ..
n+1 umn
um
So, the scheme may be written in the form un+1 = Aun , where A = C −1 B. By using the Fourier
transform, i.e.
ujn+1 = g (ξ )ujn , ujn = eij∆xξ ,
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = ,z = 2 (cos (∆xξ ) − 1).
1 − 12 z ∆x2
Moreover, g (ξ ) < 1, therefore, ρ (A) < 1.
kAxk2 ≤ kAk2 kxk2 = ρ (A) kxk2 ≤ kxk2 .
2. the scheme
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1
2
can be rewritten as
µ n+1 µ n+1 µ n µ
(1 + µ)ujn+1 = u + uj +1 + uj−1 + (1 − µ)ujn + ujn+1 .
2 j−1 2 2 2
then
µ µ µ
(1 − µ) u n + µ u n .
1 + µ ujn+1 ≤ uj−1
n+1 n+1
n
+ u +
j +1 j−1 u + j j +1
2 2 2 2
Therefore
µ
n+1
µ
n+1
µ
n
µ
n
(1 + µ)
ujn+1
≤
uj−1
n
+
uj +1
+
uj−1
+ (1 − µ)
uj
+
uj +1
.
∞ 2 ∞ 2 ∞ 2 ∞ ∞ 2 ∞
i.e.
µ
µ
µ µ
(1 + µ)
un+1
∞ ≤
un+1
∞ +
un+1
∞ + kun k∞ + (1 − µ) kun k∞ + kun k∞ .
2 2 2 2
if µ ≤ 1, then
un+1
≤ kun k∞ ,
∞
i.e.
kAun k∞ ≤ kun k∞ .
a2 (∆t )2 n a∆t
ujn+1 = ujn +
2
uj−1 − 2ujn + ujn+1 − ujn+1 − uj−1
n
,
2(∆x ) 2∆x
for the approximating the solution of the Cauchy problem for the advection equation
∂u ∂u
+a = 0, a > 0.
∂t ∂x
Use Von Neumann’s Method to show that the Lax-Wendroff scheme is stable provided the CFL condition
a∆t
≤ 1.
∆x
is enforced.
then we have
a2 (∆t )2 n a∆t
g (ξ )ujn = ujn + u − 2u n
+ u n
j +1 − u n
− u n
j−1 .
2(∆x )2 j−1 j 2∆x j +1
And then
a2 (∆t )2 i (j−1)∆xξ a∆t
g (ξ )eij∆xξ = eij∆xξ + 2
e − 2eij∆xξ + ei (j +1)∆xξ − ei (j +1)∆xξ − ei (j−1)∆xξ .
2(∆x ) 2∆x
Therefore
a2 (∆t )2 −i∆xξ i∆xξ
a∆t
i∆xξ −i∆xξ
g (ξ ) = 1+ e − 2 + e − e − e
2(∆x )2 2∆x
a2 (∆t )2 a∆t
= 1+ 2
(2 cos(∆xξ ) − 2) − (2i sin(∆xξ ))
2(∆x ) 2∆x
a2 (∆t )2 a∆t
= 1+ (cos(∆xξ ) − 1) − (i sin(∆xξ )) .
(∆x )2 ∆x
a∆t
Let µ = ∆x , then
g (ξ ) = 1 + µ2 (cos(∆xξ ) − 1) − µ (i sin(∆xξ )) .
If g (ξ ) < 1, then the scheme is stable, i,e
2
1 + µ2 (cos(∆xξ ) − 1) + (µ sin(∆xξ ))2 < 1.
i.e.
i.e.
i.e.
i.e
Problem A.34. (Sample #16) Consider the Crank-Nicholson scheme applied to the diffusion equation
∂u ∂2 u
= 2
∂t ∂x
where t > 0, −∞ < x < ∞.
1. Show that the amplification factor in the Von Neumann analysis of the scheme us
1 + 21 z ∆t
g (ξ ) = ,z = 2 (cos (∆xξ ) − 1).
1 − 12 z ∆x2
∆t
Let µ = ∆x2
, then the scheme can be rewrote as
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1 ,
2
i.e.
µ n+1 µ µ n µ
− uj−1 + (1 + µ)ujn+1 − ujn++11 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
By using the Fourier transform, i.e.
ujn+1 = g (ξ )ujn , ujn = eij∆xξ ,
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = ,z = 2 (cos (∆xξ ) − 1).
1 − 12 z ∆x2
∆t
2. since z = 2 ∆x 2 (cos (∆xξ ) − 1), then z < 0, then we have
1 1
1 + z < 1 − z,
2 2
therefore g (ξ ) < 1. Since −1 < 1, then
1 1
z − 1 < z + 1.
2 2
Therefore,
1 + 21 z
g (ξ ) = > −1.
1 − 12 z
hence g (ξ ) < 1. So, the scheme is stable.
J
∆t 1 t∗
where b > 0, µ = (∆x )2
, ∆x = L+ 1 , and ∆t = N. Prove that, under suitable restrictions on µ and ∆x, the
n
error grid function e satisfy the estimate
ken k∞ ≤ t ∗ C ∆t + ∆x2 ,
Solution. Let ū be the exact solution and ūjn = ū (n∆t, j∆x ). Then from Taylor Expansion, we have
∂ n 1 ∂2
ūjn+1 = ūjn + ∆t
ūj + (∆t )2 2 ū (ξ1 , j∆x ), tn ≤ ξ1 ≤ tn+1 ,
∂t 2 ∂t
∂ 1 ∂ 3 1 ∂4
n
ūj−1 = ūjn − ∆x ūjn − (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ2 ), xj−1 ≤ ξ2 ≤ xj ,
∂x 6 ∂x 24 ∂x
∂ 1 ∂ 3 1 ∂4
ūjn+1 = ūjn + ∆x ūjn + (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ3 ), xj ≤ ξ3 ≤ xj + 1 .
∂x 6 ∂x 24 ∂x
Then the truncation error T of this scheme is
ūjn+1 − ūjn n
ūj−1 − 2ūjn + ūjn+1 ūjn+1 − ūj−1
n
T = − +b
∆t ∆x2 ∆x
= O (∆t + (∆x )2 ).
Therefore
bµ∆x
ejn+1 = ejn + µ ej−1
n
− 2ejn + ejn+1 − ejn+1 − ej−1
n
+ c∆t (∆t + (∆x )2 ),
2
i.e.
! !
bµ∆x n bµ∆x n
ejn+1 = µ+ n
ej−1 + (1 − 2µ)ej + µ − ej +1 + c∆t (∆t + (∆x )2 ).
2 2
Then
en+1 ≤ µ + bµ∆x en + (1 − 2µ) en + µ − bµ∆x en + c∆t (∆t + (∆x )2 ).
j 2 j−1 j 2 j +1
Therefore
en+1
≤ µ + bµ∆x
en
+ (1 − 2µ)
en
+ µ − bµ∆x
en
+ c∆t (∆t + (∆x )2 ).
j ∞ 2 j−1 ∞ j ∞ 2 j +1 ∞
bµ∆x n bµ∆x n
en+1
≤ µ +
ke k∞ + (1 − 2µ) ken k∞ + µ − ke k∞ + c∆t (∆t + (∆x )2 ).
∞ 2 2
bµ∆x 1
If 1 − 2µ ≥ 0 and µ − 2 ≥ 0, i.e. µ ≤ 2 and 1 − 12 b∆x > 0, then
! !
bµ∆x bµ∆x
en+1
∞
≤ µ+ ken k∞ + ((1 − 2µ)) ken k∞ + µ − ken k∞ + c∆t (∆t + (∆x )2 )
2 2
= ken k∞ + c∆t (∆t + (∆x )2 ).
Then
ken k∞ ≤
en−1
+ c∆t (∆t + (∆x )2 )
∞
≤
en−2
+ c2∆t (∆t + (∆x )2 )
∞
..
≤ .
≤
e0
+ cn∆t (∆t + (∆x )2 )
∞
= ct ∗ (∆t + (∆x )2 ).
Problem B.1. (Prelim Jan. 2011#1) Consider a linear system Ax = b with A ∈ Rn×n . Richardson’s method
is an iterative method
Mxk +1 = N xk + b
1
Solution. 1. Since M = w,N = M − A = w1 I − A, then we have
xk +1 = (I − wA)xk + bw.
So TR = I − wA, From the sufficient and & necessary condition for convergence, we should have
ρ (TR ) < 1. Since λi are the eigenvalue of A, then we have 1 − λi w are the eigenvalues of TR . Hence
Richardson’s method converges if and only if |1 − λi w| < 1, i.e
λn w − 1 = 1 − λ1 w,
i,e
2
w0 = .
λ1 + λn
J
|1 − λn | |1 − λ1 |
1
w
1 wopt 1
λn λ1
Problem B.2. (Prelim Jan. 2011#2) Let A ∈ Cm×n and b ∈ Cm . Prove that the vector x ∈ Cn is a least
squares solution of Ax = b if and only if r⊥ range(A), where r = b − Ax.
Problem B.3. (Prelim Jan. 2011#3) Suppose A, B ∈ Rn×n and A is non-singular and B is singular. Prove
that
1 kA − Bk
≤ ,
κ (A) kAk
where κ (A) = kAk ·
A−1
, and k·k is an reduced matrix norm.
Solution. Since B is singular, then there exists a vector x , 0, s.t. Bx = 0. Since A is non-singular, then
A−1 is also non-singular. Moreover, A−1 Bx = 0. Then, we have
x = x − A−1 Bx = (I − A−1 B)x.
So
kxk =
(I − A−1 B)x
≤
A−1 A − A−1 B
kxk ≤
A−1
kA − Bk kxk .
Since x , 0, so
1 ≤
A−1
kA − Bk .
1 kA − Bk
≤ ,
A−1
kAk kAk
i.e.
1 kA − Bk
≤ .
κ (A) kAk
J
Problem B.4. (Prelim Jan. 2011#4) Let f : Ω ⊂ Rn → Rn be twice continuously differentiable. Suppose
x∗ ∈ Ω is a solution of f (x ) = 0, and the Jacobian matrix of f, denoted Jf , is invertible at x∗ .
1. Prove that if x0 ∈ Ω is sufficiently close to x∗ , then the following iteration converges to x∗ :
Solution. Let x∗ be the root of f(x ) i.e. f(x∗ )=0. From the Newton’s scheme, we have
h i−1
xk +1 = xk − J (x0 ) f(xk )
x∗ = x∗
Therefore, we have
h i−1
x∗ − xk +1 = x∗ − xk + J (x0 ) (f(xk ) − f(x∗ ))
h i−1
= x∗ − xk − J(x0 ) J(ξ )(x∗ − xk ).
therefore
J(ξ ) ∗
x∗ − xk +1 ≤ 1 −
0
x − xk
J(x )
From theorem
Theorem B.1. Suppose J : Rm → Rn×n is a continuous matrix-valued function. If J(x*) is nonsingular, then
there exists δ > 0 such that, for all x ∈ Rm with kx − x∗ k < δ, J(x) is nonsingular and
J (x )−1
< 2
J (x∗ )−1
.
we get
1
x∗ − xk +1 ≤ x∗ − xk .
2
Which also shows the convergence is typically only linear.
J
where f : [t0 , t ∗ ] × R → R is continuous in its first variable and Lipschitz continuous in its second variable.
Prove that Euler’s method converges.
So,
y (tn+1 ) − y (tn ) − hf (tn , y (tn )) = y (tn ) + hy 0 (tn ) + O (h2 ) − y (tn ) − hf (tn , y (tn ))
= y (tn ) + hy 0 (tn ) + O (h2 ) − y (tn ) − hy 0 (tn ) (219)
= O (h2 ).
Therefore,
en+1
≤ ken k + hλ ken k + ch2
= (1 + hλ) ken k + ch2 .
Claim:[2]
c
ken k ≤ h[(1 + hλ)n − 1], n = 0, 1, · · ·
λ
Proof for Claim (221): The proof is by induction on n.
1. when n = 0, en = 0, hence ken k ≤ λc h[(1 + hλ)n − 1],
2. Induction assumption:
c
ken k ≤ h[(1 + hλ)n − 1]
λ
3. Induction steps:
en+1
≤(1 + hλ) ken k + ch2
c
≤ (1 + hλ) h[(1 + hλ)n − 1] + ch2
λ
c
= h[(1 + hλ)n+1 − 1].
λ
So, from the claim (221), we get ken k → 0, when h → 0. Therefore Forward Euler Method is convergent
. J
what’s the order of the scheme? Is it a convergent scheme? Is it A-stable? Justify your answers.
So,
ξ2 ξ3
ρ (w ) − σ (w )ln(w ) = ξ 2 + 3ξ − (3 + 3ξ + ξ 2 )(ξ − + ···)
2 3
+3ξ +ξ 2
−3ξ −3ξ 2 −ξ 3
=
+ 23 ξ 2 + 32 ξ 3 + 12 ξ 4
−ξ 3 −ξ 4 − 31 ξ 5
1
= − ξ 2 + O (ξ 3 ).
2
Therefore, by the theorem
1
ρ (w ) − σ (w )ln(w ) = − ξ 2 + O (ξ 3 ).
2
Hence, this scheme is order of 1. Since,
s
X
ρ (w ) : = am wm = −2 + w + w2 = (w + 2)(w − 1). (223)
m=0
And w = −1 or w = −2 which does not satisfy the root condition. Therefore, this scheme is not stable.
Hence, it is also not A-stable. J
Problem B.7. (Prelim Jan. 2011#7) Consider the Crank-Nicholson scheme applied to the diffusion equa-
tion
∂u ∂2 u
= 2
∂t ∂x
where t > 0, −∞ < x < ∞.
1. Show that the amplification factor in the Von Neumann analysis of the scheme us
1 + 21 z ∆t
g (ξ ) = 1
,z = 2 (cos (∆xξ ) − 1).
1− 2z
∆x2
ujn+1 − ujn
n+1 n+1 n+1 n
− 2ujn + ujn+1
1 uj−1 − 2uj + uj +1 uj−1
= +
∆t 2 ∆x2 ∆x2
∆t
Let µ = ∆x2
, then the scheme can be rewrote as
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1 ,
2
i.e.
µ n+1 µ µ n µ
− uj−1 + (1 + µ)ujn+1 − ujn++11 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
By using the Fourier transform, i.e.
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = 1
,z = 2 (cos (∆xξ ) − 1).
1− 2z
∆x2
∆t
2. since z = 2 ∆x 2 (cos (∆xξ ) − 1), then z < 0, then we have
1 1
1 + z < 1 − z,
2 2
therefore g (ξ ) < 1. Since −1 < 1, then
1 1
z − 1 < z + 1.
2 2
Therefore,
1 + 21 z
g (ξ ) = > −1.
1 − 12 z
hence g (ξ ) < 1. So, the scheme is stable.
J
∆t 1 t∗
where b > 0, µ = (∆x )2
, ∆x = L+ 1 , and ∆t = N. Prove that, under suitable restrictions on µ and ∆x, the
n
error grid function e satisfy the estimate
ken k∞ ≤ t ∗ C ∆t + ∆x2 ,
Solution. Let ū be the exact solution and ūjn = ū (n∆t, j∆x ). Then from Taylor Expansion, we have
∂ n 1 ∂2
ūjn+1 = ūjn + ∆t
ūj + (∆t )2 2 ū (ξ1 , j∆x ), tn ≤ ξ1 ≤ tn+1 ,
∂t 2 ∂t
∂ 1 ∂ 2 1 ∂ 3 1 ∂4
n
ūj−1 = ūjn − ∆x ūjn + (∆x )2 2 ūjn − (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ2 ), xj−1 ≤ ξ2 ≤ xj ,
∂x 2 ∂x 6 ∂x 24 ∂x
∂ 1 ∂ 2 1 ∂ 3 1 ∂4
ūjn+1 = ūjn + ∆x ūjn + (∆x )2 2 ūjn + (∆x )3 3 ūjn + (∆x )4 4 ū (n∆t, ξ3 ), xj ≤ ξ3 ≤ xj + 1 .
∂x 2 ∂x 6 ∂x 24 ∂x
Then the truncation error T of this scheme is
ūjn+1 − ūjn n
ūj−1 − 2ūjn + ūjn+1 ūjn+1 − ūj−1
n
T = − +b
∆t ∆x2 ∆x
= O (∆t + (∆x )2 ).
Therefore
bµ∆x
ejn+1 = ejn + µ ej−1
n
− 2ejn + ejn+1 − ejn+1 − ej−1
n
+ c∆t (∆t + (∆x )2 ),
2
i.e.
! !
bµ∆x n bµ∆x n
ejn+1 = µ+ n
ej−1 + (1 − 2µ)ej + µ − ej +1 + c∆t (∆t + (∆x )2 ).
2 2
Then
en+1 ≤ µ + bµ∆x en + (1 − 2µ) en + µ − bµ∆x en + c∆t (∆t + (∆x )2 ).
j 2 j−1 j 2 j +1
Therefore
en+1
≤ µ + bµ∆x
en
+ (1 − 2µ)
en
+ µ − bµ∆x
en
+ c∆t (∆t + (∆x )2 ).
j ∞ 2 j−1 ∞ j ∞ 2 j +1 ∞
bµ∆x n bµ∆x n
en+1
≤ µ +
ke k∞ + (1 − 2µ) ken k∞ + µ − ke k∞ + c∆t (∆t + (∆x )2 ).
∞ 2 2
bµ∆x
If 1 − 2µ ≥ 0 and µ − 2 ≥ 0, i.e. µ ≤ 12 and 1 − 12 b∆x > 0, then
! !
bµ∆x bµ∆x
en+1
∞
≤ µ+ ken k∞ + ((1 − 2µ)) ken k∞ + µ − ken k∞ + c∆t (∆t + (∆x )2 )
2 2
= ken k∞ + c∆t (∆t + (∆x )2 ).
Then
ken k∞ ≤
en−1
+ c∆t (∆t + (∆x )2 )
∞
≤
e
+ c2∆t (∆t + (∆x )2 )
n−2
∞
..
≤ .
≤
e0
+ cn∆t (∆t + (∆x )2 )
∞
= ct ∗ (∆t + (∆x )2 ).
J
Problem B.9. (Prelim Aug. 2010#1) Prove that A ∈ Cm×n (m > n) and let A = Q̂R̂ be a reduced QR
factorization.
1. Prove that A has rank n if and only if all the diagonal entries of R̂ are non-zero.
2. Suppose rank(A) = n, and define P = Q̂Q̂∗ . Prove that range(P ) = range(A).
3. What type of matrix is P?
Solution. 1. From the properties of reduced QR factorization, we knowQthat Q̂ has orthonormal columns,
therefore det(Q̂ ) = 1 and R̂ is upper triangular matrix, so det(R̂) = ni=1 rii . Then
n
Y
det(A) = det(Q̂R̂) = det(Q̂ ) det(R̂) = rii .
i =1
Therefore, A has rank n if and only if all the diagonal entries of R̂ are non-zero.
2. (a) range(A) ⊆ range(P ): Let y ∈ range(A), that is to say there exists a x ∈ Cn s.t. Ax = y. Then by
reduced QR factorization we have y = Q̂R̂x. then
P y = P Q̂R̂x = Q̂Q̂∗ Q̂R̂x = Q̂R̂x = Ax = y.
therefore y ∈ range(P ).
(b) range(P ) ⊆ range(A): Let v ∈ range(P ), that is to say there exists a v ∈ Cn , s.t. v = P v = Q̂Q̂∗ v.
Claim B.1.
Q̂Q̂∗ = A (A∗ A)−1 A∗ .
Proof.
−1
A (A∗ A)−1 A∗
= Q̂R̂ R̂∗ Q̂∗ Q̂R̂ R̂∗ Q̂∗
−1
= Q̂R̂ R̂∗ R̂ R̂∗ Q̂∗
−1
= Q̂R̂R̂−1 R̂∗ R̂∗ Q̂∗
= Q̂Q̂∗ .
J
Therefore by the claim, we have
Problem B.10. (Prelim Aug. 2010#4) Prove that A ∈ Rn×n is SPD if and only if it has a Cholesky factor-
ization.
A = LU = U T U .
xT Ax = xT U T U x = (U x )T U x.
with equality only when y = 0, i.e. x=0 (since U is non-singular). Hence A is SPD.
J
kAxk2 ≤ kxk2 ,
kAxk∞ ≤ kxk∞ ,
for any x ∈ Rm , provided µ ≤ 1.(In other words, the scheme may only be conditionally stable in the
max norm.)
n+1 n
u1 u1
n+1 u n
u2 2
un+1 = . and un = . .
.. ..
n+1 umn
um
So, the scheme may be written in the form un+1 = Aun , where A = C −1 B. By using the Fourier
transform, i.e.
ujn+1 = g (ξ )ujn , ujn = eij∆xξ ,
then we have
µ n µ µ n µ
− g (ξ )uj−1 + (1 + µ)g (ξ )ujn − g (ξ )ujn+1 = uj−1 + (1 − µ)ujn + ujn+1 .
2 2 2 2
And then
µ µ µ µ
− g (ξ )ei (j−1)∆xξ + (1 + µ)g (ξ )eij∆xξ − g (ξ )ei (j +1)∆xξ = ei (j−1)∆xξ + (1 − µ)eij∆xξ + ei (j +1)∆xξ ,
2 2 2 2
i.e.
µ µ µ µ
g (ξ ) − e−i∆xξ + (1 + µ) − ei∆xξ ej∆xξ = e−i∆xξ + (1 − µ) + ei∆xξ ej∆xξ ,
2 2 2 2
i.e.
g (ξ ) (1 + µ − µ cos(∆xξ )) = 1 − µ + µ cos(∆xξ ).
therefore,
1 − µ + µ cos(∆xξ )
g (ξ ) = .
1 + µ − µ cos(∆xξ )
hence
1 + 12 z ∆t
g (ξ ) = ,z = 2 (cos (∆xξ ) − 1).
1 − 12 z ∆x2
Moreover, g (ξ ) < 1, therefore, ρ (A) < 1.
kAxk2 ≤ kAk2 kxk2 = ρ (A) kxk2 ≤ kxk2 .
2. the scheme
µ n+1
ujn+1 = ujn + uj−1 − 2ujn+1 + ujn++11 + uj−1
n
− 2ujn + ujn+1
2
can be rewritten as
µ n+1 µ n+1 µ n µ
(1 + µ)ujn+1 = u + uj +1 + uj−1 + (1 − µ)ujn + ujn+1 .
2 j−1 2 2 2
then
µ µ µ
(1 − µ) u n + µ u n .
1 + µ ujn+1 ≤ uj−1
n+1 n+1
n
+ u +
j +1 j−1 u + j j +1
2 2 2 2
Therefore
µ
n+1
µ
n+1
µ
n
µ
n
(1 + µ)
ujn+1
≤
uj−1
n
+
uj +1
+
uj−1
+ (1 − µ)
uj
+
uj +1
.
∞ 2 ∞ 2 ∞ 2 ∞ ∞ 2 ∞
i.e.
µ
µ
µ µ
(1 + µ)
un+1
∞ ≤
un+1
∞ +
un+1
∞ + kun k∞ + (1 − µ) kun k∞ + kun k∞ .
2 2 2 2
if µ ≤ 1, then
un+1
≤ kun k∞ ,
∞
i.e.
kAun k∞ ≤ kun k∞ .
a2 (∆t )2 n a∆t
ujn+1 = ujn +
2
uj−1 − 2ujn + ujn+1 − ujn+1 − uj−1
n
,
2(∆x ) 2∆x
for the approximating the solution of the Cauchy problem for the advection equation
∂u ∂u
+a = 0, a > 0.
∂t ∂x
Use Von Neumann’s Method to show that the Lax-Wendroff scheme is stable provided the CFL condition
a∆t
≤ 1.
∆x
is enforced.
then we have
a2 (∆t )2 n a∆t
g (ξ )ujn = ujn + u − 2u n
+ u n
j +1 − u n
− u n
j−1 .
2(∆x )2 j−1 j 2∆x j +1
And then
a2 (∆t )2 i (j−1)∆xξ i (j +1)∆xξ
a∆t
i (j +1)∆xξ i (j−1)∆xξ
g (ξ )eij∆xξ = eij∆xξ + e − 2e ij∆xξ
+ e − e − e .
2(∆x )2 2∆x
Therefore
a2 (∆t )2 −i∆xξ a∆t
g (ξ ) = 1+ 2
e − 2 + ei∆xξ − ei∆xξ − e−i∆xξ
2(∆x ) 2∆x
a2 (∆t )2 a∆t
= 1+ (2 cos(∆xξ ) − 2) − (2i sin(∆xξ ))
2(∆x )2 2∆x
a2 (∆t )2 a∆t
= 1+ (cos(∆xξ ) − 1) − (i sin(∆xξ )) .
(∆x )2 ∆x
a∆t
Let µ = ∆x , then
g (ξ ) = 1 + µ2 (cos(∆xξ ) − 1) − µ (i sin(∆xξ )) .
If g (ξ ) < 1, then the scheme is stable, i,e
2
1 + µ2 (cos(∆xξ ) − 1) + (µ sin(∆xξ ))2 < 1.
i.e.
i.e.
i.e.
i.e
then we get µ < 1. The above process is invertible, therefore, we prove the result. J
Problem B.13. (Prelim Jan. 2008#8) Let Ω ⊂ R2 be a bounded domain with a smooth boundary. Consider
a 2-D poisson-like equation
−∆u + 3u = x2 y 2 , in Ω,
= 0, on ∂Ω.
u
ku − uh kH 1 ≤ C inf ku − vh kH 1 ,
vh ∈Vh
with C independent of h, where Vh denotes a finite element subspace of H 1 (Ω) consisting of contin-
uous piecewise polynomials of degree k ≥ 1.
Solution. 1. For this pure Dirichlet Problem, the test functional space v ∈ H01 . Multiple the test func-
tion on the both sides of the original function and integral on Ω, we get
Z Z Z
− ∆uvdx + uvdx = xyvdx.
Ω Ω Ω
Integration by part yields
Z Z Z
∇u∇vdx + uvdx = xyvdx.
Ω Ω Ω
Let
Z Z Z
a(u, v ) = ∇u∇vdx + uvdx, f (v ) = xyvdx.
Ω Ω Ω
Then, the
(a) Ritz variational problem is: find uh ∈ H01 , such that
1
J (uh ) = min a(uh , uh ) − f (uh ).
2
(b) Galerkin variational problem is: find uh ∈ H01 , such that
a(uh , uh ) = f (uh ).
2. Next, we will use Lax-Milgram to prove the uniqueness.
(a)
Z Z
a(u, v ) ≤ |∇u∇v| dx + |uv| dx
Ω Ω
≤ k∇ukL2 (Ω) k∇vkL2 (Ω) + kukL2 (Ω) kvkL2 (Ω)
≤ k∇ukL2 (Ω) k∇vkL2 (Ω) + C k∇ukL2 (Ω) k∇vkL2 (Ω)
≤ C k∇ukL2 (Ω) k∇vkL2 (Ω)
≤ C kukH 1 (Ω) kvkH 1 (Ω)
(b)
Z Z
2
a(u, u ) = (∇u ) dx + u 2 dx
Ω Ω
So,
Z Z
a(u, u ) = |∇u|2 dx + |u|2 dx
Ω Ω
= k∇uk2L2 (Ω) + kuk2L2 (Ω)
= kuk2H 1 (Ω) .
(c)
Z
f (v ) ≤ xyv dx
Ω
Z
≤ max |xy| |v| dx
Ω
Z !1/2 Z !1/2
≤ C 12 dx |v|2 dx
Ω Ω
≤ C kvkL2 (Ω) ≤ C kvkH 1 (Ω) .
by Lax-Milgram theorem, we get that e Galerkin method has a unique solution uh . Moreover,
a(vh , vh ) = f (vh ).
a(u, vh ) = f (vh ).
a(u − uh , vh ) = 0.
Then by coercivity
ku − uh k2H 1 (Ω) ≤ a(u − uh , u − uh )
= a(u − uh , u − vh ) + a(u − uh , vh − uh )
= a(u − uh , u − vh )
≤ ku − uh kH 1 (Ω) ku − vh kH 1 (Ω) .
Therefore,
ku − uh kH 1 ≤ C inf ku − vh kH 1 ,
vh ∈Vh
C Project 1 MATH571
COMPUTATIONAL ASSIGNMENT # 1
MATH 571
1. Instability of Gram–Schmidt
The purpose of the first part of your assignment is to investigate the instability of the
classical Gram–Schmidt orthogonalization process. Lecture 9 in BT is somewhat related to
this and could be a good source of inspiration.
1.– Write a piece of code that implements the classical Gram–Schmidt process, see Algo-
rithm 7.1 in BT. Ideally, this should be implemented in the form of a QR factorization,
that is, given a matrix A ∈ Rm×n your method should return two matrices Q ∈ Rm×n
and R ∈ Rn×n , where the matrix Q has (or at least should have) orthonormal columns
and A = QR.
2.– With the help of the developed piece of code, test the algorithm on a matrix A = R20×10
with:
• entries uniformly distributed over the interval [0, 1].
• entries given by
j−1
2i − 21
ai,j = .
19
• entries given by
1
ai,j = ,
i+j−1
this is the so-called Hilbert matrix.
3.– For each one of these cases compute Q? Q. Since Q, in theory, has orthonormal columns
what should you get? What do you actually get?
4.– Implement the modified Gram-Schmidt process (Algorithm 8.1 in BT) and repeat steps
1.—3. What do you observe?
TTH 12:40pm
Wenqiang Feng
Contents
Problem 1 3
Problem 2 5
Problem 1
1. See Listing 3.
2. See Listing 2.
3. We should get the Identity square matrices. But we did not get the actual Identity square matrices
through the Gram-Schmidt Algorithm. For case 1-2, we only get the matrices which diag(Q∗ Q) =
n
z }| {
{1, · · · , 1} and the other elements approximate to 0 in the sense of C × 10−16 ∼ 10−17 . For case 3,
Classical Gram-Schmidt Algorithm is not stable for case 3, since some elements of matrix Q∗ Q do not
approximate to 0, then the matrix Q∗ Q is not diagonal any more.
4. For case 1-2, we also did not get the actual Identity square matrices by using the Modified Gram-
n
z }| {
∗
Schmidt Algorithm. We only get the matrices which diag(Q Q) = {1, · · · , 1} and the other elements
approximate to 0 in the sense of C × 10−17 ∼ 10−18 . For case 3, the Modified Gram-Schmidt Algo-
n
z }| {
∗
rithm works well for case 3, we get the matrix which diag(Q Q) = {1, · · · , 1} and the other elements
approximate to 0 in the sense of C × 10−8 ∼ 10−13 . So, Modified Gram-Schmidt Algorithm is more
stable than the Classical one.
[m,n]=size(V);
Q=zeros(m,n);
R=zeros(n);
R(1,1)=norm(V(:,1));
25 Q(:,1)=V(:,1)/R(1,1);
for k=2:n
R(1:k-1,k)=Q(:,1:k-1)’*V(:,k);
Q(:,k)=V(:,k)-Q(:,1:k-1)*R(1:k-1,k);
R(k,k)=norm(Q(:,k));
30 if R(k,k) == 0
break;
end
Q(:,k)=Q(:,k)/R(k,k);
end
[m,n]=size(V);
Q=zeros(m,n);
R=zeros(n);
25 for k=1:n
R(k,k)=norm(V(:,k));
i f R(k,k) == 0
break;
end
30 Q(:,k)=V(:,k)/R(k,k);
for j=k+1:n
R(k,j)=Q(:, k)’* V(:,j);
V(:, j) = V(:, j)-R(k, j) * Q(:,k);
end
35 end
Problem 2
1. I Chose N = 10 and I got the polynomial p10 is as follow:
2. See Figure 1.
3. See Listing 6.
4. Since A is Vandermonde Matrix and all the points xi are different, then det(A) 6= 0. Therefore A has
full rank. Page 168 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 169
Wenqiang Feng MATH 571 ( TTH 12:40pm): Coding Assignment #1 Problem 2 (continued)
5. I varied N from 3 to 15. For every fixed N , I varied n form 1 to N . Then I got the following table
√
(Table.1). From table (Table.1), we can get that n ≈ 2 N + 1, where the N is the number of the
partition.
N \h 1 2 3 4 5 6 7 8 ··· 1
−17 −17
3 0.23 3.96 · 10 5.55 · 10
4 0.82 0.56 0.56 5.10 · 10−17
5 0.50 0.28 0.28 9.04 · 10−16 9.32 · 10−16
6 0.84 0.62 0.62 0.43 0.43 8.02 · 10−15
7 0.71 0.46 0.46 0.25 0.25 3.32 · 10−15 3.96 · 10−15
8 0.89 0.64 0.64 0.45 0.45 0.30 0.30 1.39 · 10−14
..
.
10 x1=-1:2/(2*N):1;
a = polyfit(x,y,n);
p = polyval(a,x1)
plot(x,y,’o’,x1,p,’-’)
15 for m=1:10
least_squares(x, y, m)
end
x=-1:2/N:1;
b=fun(x);
10
A=MatrixGen(x,n);
cof=GSsolver(A,b);
q=0;
for i=1:n+1
15 q=q+cof(i)*(x.ˆ(i-1));
end
error(j)=norm(q-b);
j=j+1;
error
20 end
end
function A=MatrixGen(x,n)
25 m=size(x,2);
A=zeros(m,n+1);
for i=1:m
for j=1:n+1
A(i,j)=x(i).ˆ(j-1);
30 end
end
function x=GSsolver(A,b)
35 [Q,R]=mgschmidt(A);
x= R\(Q’*b’);
D Project 2 MATH571
COMPUTATIONAL ASSIGNMENT # 2
MATH 571
kxk+1 − xk k < ,
or, since we are just trying to learn,
k+1
kx
Page − xk <
178 of,236
where x is the exact solution.
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 179
3. Set A = Λ + Σ.
4. Are the matrices Λ and Σ SPD? Do they commute?
5. Choose x ∈ RN and set f = Ax. Run ADI(A, f, , Λ, Σ, τ ) for different values of τ . Which one seems to
be the optimal one?
The following are not obligatory but can be used as extra credit:
6. Write an expression for xk+1 in terms of xk only. Hint: Try adding and subtracting (1) and (2). What
do you get?
7. From this expression find the equation that controls the error ek = x − xk .
8. Assume that (A1 A2 x, x) ≥ 0, show that in this case [x, y] = (A1 A2 x, y) is an inner product. If that is the
case we will denote kxkB = [x, x]1/2 .
9. Under this assumption we will show convergence of the ADI scheme. To do so:
• Take the inner product of the equation that controls the error with ek+1 + ek .
• Add over k = 1, K. We should obtain
K
X
keK+1 k22 + τ kek+1 + ek k2A + 2τ 2 keK+1 k2B = ke0 k22 + 2τ 2 ke0 k2B .
k=1
TTH 12:40pm
Wenqiang Feng
Contents
Problem 1 3
Problem 2 8
Let N dim to be the Dimension of the matrix and N iter to be the iterative numbers. In the whole report, b
was generated by Ax, where x is a corresponding vector and x’s entries are random numbers between 0 and
10. The initial iterative values of x are given by ~0.
Problem 1
1. Listing 1 shows the implement of Jacobi Method.
(a) From the records of the iterative number, I got the following results:
For case (2), the Jacobi Method is not convergence, because it has a big Condition Number. For
case (1) and case (3), if N dim is small, roughly speaking, N dim ≤ 10 − 20, then the N dim and
N iter have the roughly relationship N iter = log(N dim + C), when N dim is large, the N iter is
not depent on the N dim (see Figure (1)).
(b) When ω = 1, the SOR Method degenerates to the Gauss-seidel Method. For Gauss-seidel Method,
I get the similar results as Jacobi Method (see Figure (2)). But, the Gauss-Seidel Method is more
stable than Jacobi Method and case (3) is more stable than case (1) (see Figure (1) and Figure
(2)).
Jacobi iteration
45 GS iteration with
40
The iterative steps
35
30
25
10 20 30 40 50 60 70 80 90
The value of N
100
90
Jacobi iteration
GS iteration with
80
The iterative steps
70
60
50
40
30
10 20 30 40 50 60 70 80 90
The value of N
Figure 2: The relationship between N dim and N iter for case (3)
45
40
The iterative steps
Jacobi iteration
35 GS iteration with
SOR iteration with w=0.999
30
25
10 20 30 40 50 60 70 80 90
The value of N
ii. For case (2), In general, the SOR Method is not convergence, but SOR is convergence for
some small N dim ; Page 183 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 184
Wenqiang Feng MATH 571 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
Jacobi iteration
50 GS iteration with
SOR iteration with w=1.001
45
The iterative steps
40
35
30
25
10 20 30 40 50 60 70 80 90
The value of N
iii. For case(3), the optimal w is around 1.14; This numerical result is same as the theoretical
result. Let D = diag(diag(A)); E = A − D; T = D\E,
2
wopt = p ≈ 1.14.
1 − ρ(T )2
100
90
80 Jacobi iteration
GS iteration with
SOR iteration with w=1.14
70
The iterative steps
60
50
40
30
10 20 30 40 50 60 70 80 90
The value of N
35 while (error>tol&&iter<max_iter)
x1=x;
x= D\(E*x+b);
error=norm(x-x1);
iter=iter+1;
40 end
Page
Listing 185 of
2: SOR 236
Method
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 186
Wenqiang Feng MATH 571 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
function [x iter]=sor(A,b,w,x,tol,max_iter)
% jacobi: Solve the linear system with SOR iterative algorithm
%
% USAGE
5 % jacobi(A,b,epsilon,x0,tol,max_iter)
%
% INPUT
% A: N by N LHS coefficients matrix
% b: N by 1 RHS vector
10 % w: Relaxation parameter
% x: Initial guess
% tol: The stop tolerance
% max_iter: maxmum iterative steps
%
15 % OUTPUT
% x: The solutions
% iter: iterative steps
%
% AUTHOR
20 % Wenqiang Feng
% Department of Mathematics
% University of Tennessee at Knoxville
% E-mail: wfeng@math.utk.edu
% Date: 11/13/2013
25
n=size(A,1);
% Set default parameters
i f (nargin<4), x=zeros(n,1);tol=1e-16;max_iter=500;end;
%Initial some parameters
30 error=norm(b - A*x)/norm( b );
iter=0 ;
% s p l i t the matrix for Jacobi interative method
D=diag(diag( A ));
b = w * b;
35 M = w * tril( A, -1 ) + D;
N = -w * triu( A, 1 ) + ( 1.0 - w ) * D;
while (error>tol&&iter<max_iter)
x1=x;
40 x= M\(N*x+b);
error=norm(x-x1)/norm( x );
iter=iter+1;
end
Problem 2
1. Listing 3 shows the implement of ADI Method.
2. Yes, The Σ and Λ are the SPD matrices. Moreover, they are commute, since ΣΛ = ΛΣ.
6. Now, I will show [x, y] = (A1 A2 x, y) is an inner product, i.e, I will show the ||x||2B = [x, x] satisfies
parallelogram law:
It’s easy to show that the B-norm ||x||2B = [x, x] satisfies the parallelogram law,
So, The norm space can induce a inner product, so [x, y] = (A1 A2 x, y) is a inner product.
Therefore,
Page 187 of 236
||ek+1 ||22 + τ ||ek+1 + ek ||2A + τ 2 ||ek+1 ||2B = ||ek ||22 + τ 2 ||ek ||2B . (10)
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 188
Wenqiang Feng MATH 571 ( TTH 12:40pm): Computational Assignment #2 Problem 2 (continued)
Therefore, from (11), we get ||ek+1 + ek ||2A → 0 ∀τ > 0. So 21 (xk+1 + xk ) → x with respect to || · ||A .
while (error>tol&&iter<max_iter)
x1=x;
x=(tau*I+A1)\((tau*I-A2)*x+b); % the first half step
x=(tau*I+A2)\((tau*I-A1)*x+b); % the second half step
40 error=norm(x-x1);
iter=iter+1;
end
Page 188 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 189
TTH 12:40pm
Wenqiang Feng
Contents
Problem 1 3
Problem 2 4
Problem 3 4
Problem 1
Given the equation
(
−u00 + u = f, in Ω
0 0
(1)
−u (0) = u (1) = 0, on ∂Ω
devise a finite difference scheme for this problem that results in a tridiagonal matrix. The scheme must be
consistent of order O(2) in the C(Ω̄h ) norm and you should prove this.
Proof: I consider the following uniform partition (Figure. 1) of the interval (0, 1) with N + 1 points. For the
Neumann Boundary, we introduce two ghost point x−1 and xN +1 .
x−1 x0 = 0 x1 xN −1 xN = 1 xN +1
From the homogeneous Neumann boundary condition, we know that U1 = U−1 and UN +1 = UN −1 . Therefore
1. for i = 0, from the scheme,
1 2 1 2 2
− 2
U−1 + 2 U0 − 2 U1 + U0 = (1 + 2 )U0 − 2 U1 = F0
h h h h h
2. for i = 1, · · · , N − 1, we get
1 2 1 1 2 1
− 2
Ui−1 + 2 Ui − 2 Ui+1 + Ui = − 2 Ui−1 + (1 + 2 )Ui − 2 Ui+1 = Fi .
h h h h h h
3. for i = N
1 2 1 2 2
− UN −1 + 2 UN − 2 UN +1 + UN = (1 + 2 )UN − 2 UN −1 = FN .
h2 h h h h
So the algebraic system is
AU = F,
where
1 + h22 − h22 U0 F0
−1 1 + h22 − h12 U1 F1
h2
.. .. .. .. ..
A= . . . ,U = F = .
. .
1
− h2 1 + h22 − h12 UN −1 FN −1
− h22 1 + h22 UN FN
Next, I will show this scheme is of order O(2). From the Taylor expansion, we know
h2 00 h3
Ui+1 = u(xi+1 ) = u(xi ) + hu0 (xi ) + u (xi ) + u(3) (xi ) + O(h4 )
2 2
Page 2 of 236
192 3
h h
Ui−1 = u(xi−1 ) = u(xi ) − hu0 (xi ) + u00 (xi ) − u(3) (xi ) + O(h4 ).
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 193
Wenqiang Feng MATH 572 ( TTH 12:40pm): Exam problem 4-5 Problem 1
Therefore,
Ui+1 − 2Ui + Ui−1 u(xi+1 ) − 2u(xi ) + u(xi−1 )
− =− = −u00 (xi ) + O(h2 ).
h2 h2
Therefore, the scheme (2) is of order O(h2 ).
Problem 2
Let A = tridiag{ai , bi , ci }ni=1 ∈ Rn×n be a tridiagional matrix with the properties that
bi > 0, ai , ci ≤ 0, ai + bi + ci = 0.
Prove the following maximum principle: If u ∈ Rn is such that (Au)i=2,··· ,n−1 ≤ 0, then ui ≤ max{u1 , un }.
This is contradiction to (Au)i=2,··· ,n−1 < 0. Therefore, If u ∈ Rn is such that (Au)i=2,··· ,n−1 < 0, then
ui ≤ max{u1 , un }.
ak uk−1 + bk uk + ck uk+1 = 0.
Since ak + ck = −bk , so
And ak < 0, ck < 0, uk−1 − uk ≤ 0, uk+1 − uk ≤ 0, so uk−1 = uk = uk+1 , that is to say, uk−1 and uk+1
is also the maximum points. Bu using the same argument again, we get uk−2 = uk−1 = uk = uk+1 =
uk+2 . Repeating the process, we get
u1 = u2 = · · · = un−1 = un .
Problem 3
Prove the following discrete Poincaré inequality: Let Ω = (0, 1) and Ωh be a uniform grid of size h. If Y ∈ Uh
is a mesh function on Ωh such that Y (0) = 0, then there is a constant C, independent of Y and h, for which
kY k2,h ≤ C
δ̄ Y
2,h .
x1 = 0 x2 xN −1 xN = 1
Then,
XN
Yi−1 − Yi = |YN |.
i=2
and
N N N
!1/2 N !1/2
X X Yi−1 − Yi X X Yi−1 − Yi 2
|YN | ≤ |Yi−1 − Yi | = h ≤
h 2
.
i=2 i=2
h i=2 i=2
h
Therefore
K
! K !
X X Yi−1 − Yi 2
|YK | 2
≤ 2
h
h
i=2 i=2
K
X Yi−1 − Yi 2
2
= h (K − 1) .
h
i=2
1. When K = 2,
Y1 − Y2 2
|Y2 | 2
≤ 2
h .
h
2. When K = 3,
!
Y1 − Y2 2 Y2 − Y3 2
|Y3 | 2
≤ 2h 2
h + h .
3. When K = N ,
!
Y1 − Y2 2 Y2 − Y3 2 YN −1 − YN 2
|YN | 2
≤ (N − 1)h 2 + + ··· + .
Page
h 194 of 236
h h
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 195
Wenqiang Feng MATH 572 ( TTH 12:40pm): Exam problem 4-5 Problem 3
Since Y1 = 0, so
N
X N 2
N (N − 1) 2 X Yi−1 − Yi
|Yi |2 ≤ h .
i=1
2 i=2
h
And then
N N X N
1 X N X Yi−1 − Yi 2 1 1 Yi−1 − Yi 2
|Y |2
≤ h2 = + h2 .
2(N − 1) i=2
i
(N − 1)2 i=1 h 2 2(N − 1) i=2
h
1
Since h = N −1 , so
N N
X 1 1 X Yi−1 − Yi 2
h2 2
|Yi | ≤ + h2 .
2 2(N − 1) h
i=1 i=2
then
N XN
X 1 1 Yi−1 − Yi 2
h 2
|Yi | ≤ + h .
2 2(N − 1) h
i=1 i=2
i.e,
2 1 1
2
kY k2,h ≤ +
δ̄ Y
.
2 2(N − 1) 2,h
since N ≥ 2, so
2
2
kY k2,h ≤
δ̄ Y
2,h .
Hence,
kY k2,h ≤ C
δ̄ Y
2,h .
F Project 1 MATH572
COMPUTATIONAL ASSIGNMENT # 1
MATH 572
While this is good and we should not use methods that do not satisfy this condition, this type of result is
of little help in practice. In other words, we usually compute with a fixed h and, even if we know y(tn ), we
do not know the exact solution at the next time step and, thus, cannot assess how small the local error
en+1 = y(tn+1 ) − yn+1
is. Here we will study two strategies to estimate this quantity. Your assignment will consist in implementing
these two strategies and use them for the solution of a Cauchy problem
y 0 = f (t, y) t ∈ (t0 , T ), y(t0 ) = y0 ,
where
1. f = y − t, (t0 , T ) = (0, 10), y0 = 1 + δ, with δ ∈ {0, 10−3 }.
2. f = λy + sin t − λ cos t, (t0 , T ) = (0, 5), y0 = 0, λ ∈ {0, ±5, ±10}.
3. f = 1 − yt , (t0 , T ) = (2, 20), y0 = 2.
4. The Fresnel integral is given by
Z t
φ(t) = sin(s2 )ds.
0
Set it as a Cauchy problem and generate a table of values on [0, 10]. If possible obtain a plot of the
function.
5. The dilogarithm function Z x
ln(1 − t)
f (x) = − dt
0 t
on the interval [−2, 0].
Step Bisection. The local error analysis that is usually carried out with the help of Taylor expansions
yields, for a method of order s, that
ken+1 k ≤ Chs+1 .
The constant C here is independent of h but it might depend on the exact solution y and the current step
tn . To control the local error we will assume that C does not change as n changes. Let v denote the value of
the approximate solution at tn+1 obtained by doing one step of length h from tn . Let u be the approximate
solution at tn+1 obtained by taking two steps of size h/2 from tn . The important thing here is that both u
and v are computable. By the assumption on C we have
y(tn+1 ) = v + Chs+1 ,
y(tn+1 ) = u + 2C(h/2)s+1 ,
which implies
ku − vk
ken+1 k ≤ Chs+1 = .
1 − 2−s
Notice that the quantity on the right of this expression is completely computable. In a practical realization
one can then monitor ku − vk to make sure that it is below a prescribed tolerance. If it is not, the time step
Page 197 of 236
Date: Due March 13, 2014.
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 198
can be reduced (halved) to improve the local truncation error. On the other hand, if this quantity is well
below the prescribed tolerance, the time step can be doubled.
Implement this strategy for the fourth order ERK scheme
0
1 1
2 2
1 1
2 0 2
1 0 0 1
1 1 1 1
6 3 3 6
TTH 12:40pm
Wenqiang Feng
Contents
Adaptive Runge-Kutta Methods Formulas 3
Problem 1 3
Problem 2 4
Problem 3 7
Problem 4 8
Problem 5 8
The formula for the fourth order Runge-Kutta (4th RK) method can be read as following
y(t0 ) = y0 ,
K1 = hf (ti , yi )
K = hf (t + h , y + K1 )
2 i 2 i 2
(2)
K = hf (t + h
, y + K2
3 i 2 i 2 )
K4 = hf (ti + h, yi + K3 )
y 1
i+1 = yi + 6 (K1 + K2 + K3 + K4 )
y(t0 ) = y0 ,
K1 = hf (ti , yi )
K
h
K2 = hf (ti + 4 , yi + 41 )
3h 3 9
K3 = hf (ti + 8 , yi + 32 K1 + 32 K2 )
K4 = hf (ti + 12h 1932 7200 7296
13 , yi + 2197 K1 − 2197 K2 + 2197 K3 )
(3)
K5 = hf (ti + h, yi + 439 3680 845
216 K1 − 8K2 + 513 K3 − 4104 K4 )
K = hf (t + h , y − 8 K + 2K − 3544 K + 1859 K − 11 )
6 i i
2 27 1 2 2565 3 4104 4 40
y = y + 16
K + 6656
K + 28561
K − 9
K + 2
K
i+1 i 135 1 12825 3 56430 4 50 5 55 6
25 1408
ỹi+1 = yi + 216 K1 + 2656 K3 + 2197
4104 K 4 − 1
5 K 5 .
The error
1
E= |yi+1 − ỹi+1 | (4)
h
will be used as an estimator. If E ≤ T ol, y will be kept as the current step solution and then move to the
next step with time step size δh. If E > T ol, recalculate the current step with time step size δh, where
1/4
T ol
δ = 0.84 .
E
Problem 1
1. The 4th RK method and RKF method for Problem 1.1
(a) Results for Problem 1.1. From the figure (Fig.1) we can see that the 4th RK method and
RKF method are both convergent for Problem 1.1. The 4th RK method is convergent with 4
steps and RKF method with 2 steps and reached error 4.26 × 10−14 .
Page 201 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 202
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
10
10
9
8
8
6 6
y
y
5
4
4
3
2
1 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
x x
Figure 1: The 4th RK method and RKF method for Problem 1.1
(a) Results for Problem 1.2. From the figure (Fig.2) we can see that the 4th RK method and
RKF method are both convergent for Problem 1.2. The 4th RK method is convergent with 404
steps and reached error 9.9 × 10−6 . RKF method with 29 steps and reached error 2.3 × 10−9 .
(b) Figures (Fig.2)
30 30
25 25
20 20
y
15 15
10 10
5 5
0 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
x x
Figure 2: The 4th RK method and RKF method for Problem 1.2
Problem 2
Page 202 of 236
1. The 4th RK method and RKF method for Problem 2.1
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 203
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 2 (continued)
(a) Results for Problem 2.1. From the figure (Fig.3) we can see that the 4th RK method and
RKF method are both convergent for Problem 2.1. The 4th RK method is convergent with 24
steps and reached error 7.1 × 10−6 . RKF method with 8 steps and reached error 9.4 × 10−10 .
(b) Figures (Fig.3)
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
1 1
y
y
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x x
Figure 3: The 4th RK method and RKF method for Problem 2.1
(a) Results for Problem 2.2. From the figure (Fig.4) we can see that the 4th RK method and
RKF method are both divergent for Problem 2.2.
(b) Figures (Fig.4)
x 10
10 Problem.2.2,With steps =10002,error=8.087051e+00 x 10
6 Problem.2.2,with step=1001,error=3.725290e−09
0 0
Runge−Kutta−4th Runge−Kutta−Fehlberg
−1 −0.5
−2 −1
−3 −1.5
y
−4 −2
−5 −2.5
−6 −3
−7 −3.5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5
x x
Figure 4: The 4th RK method and RKF method for Problem 2.2
Pagefor
3. The 4th RK method and RKF method 203Problem
of 236 2.3
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 204
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 2 (continued)
(a) Results for Problem 2.3. From the figure (Fig.5) we can see that the 4th RK method and
RKF method are both convergent for Problem 2.3. The 4th RK method is convergent with 96
steps and reached error 9.98 × 10−6 . RKF method with 69 steps and reached error 1.3 × 10−11 .
(b) Figures (Fig.5)
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
y
y
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1 −1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x x
Figure 5: The 4th RK method and RKF method for Problem 2.3
(c) The 4th RK method and RKF method for Problem 2.4
i. Results for Problem 2.4. From the figure (Fig.6) we can see that the 4th RK method and
RKF method are both divergent for Problem 2.4.
ii. Figures (Fig.6)
x 10
21 Problem.2.4,With steps =10002,error=1.967067e+13 x 10
5 Problem.2.4,with step=1001,error=1.396984e−09
0 0
Runge−Kutta−4th Runge−Kutta−Fehlberg
−2
−1
−4
−2 −6
−8
−3
y
−10
−4 −12
−14
−5
−16
−6 −18
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5
x x
Figure 6: The 4th RK method and RKF method for Problem 2.4
(d) The 4th RK method and RKF method for Problem 2.5
i. Results for Problem 2.5. From
Pagethe
204figure (Fig.7) we can see that the 4th RK method
of 236
and RKF method are both convergent for Problem 2.5. The 4th RK method is convergent
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 205
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 2 (continued)
with 88 steps and reached error 8.77 × 10−6 . RKF method with 114 steps and reached error
2.57 × 10−10 .
ii. Figures (Fig.7)
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
y
y
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1 −1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x x
Figure 7: The 4th RK method and RKF method for Problem 2.5
Problem 3
1. The 4th RK method and RKF method for Problem 3
(a) Results for Problem 3. From the figure (Fig.8) we can see that the 4th RK method and RKF
method are both convergent for Problem 3. The 4th RK method is convergent with 4 steps and
reached error 1.77 × 10−15 . RKF method with 2 steps and reached error 2 × 10−15 .
(b) Figures (Fig.8)
10 10
9 9
8 8
7 7
y
6 6
5 5
4 4
3 3
2 2
2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20
x x
Page 205
Figure 8: The 4th RK method and of
RKF236method for Problem 3
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 206
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 3 (continued)
Problem 4
1. The 4th RK method and RKF method for Problem 4
(a) Results for Problem 4. From the figure (Fig.9) we can see that the 4th RK method and RKF
method are both convergent for Problem 4. The 4th RK method is convergent with 438 steps and
reached error 9.9 × 10−6 . RKF method with 134 steps and reached error 3.68 × 10−14 .
(b) Figures (Fig.9)
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
y
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
x x
Problem 5
1. The 4th RK method and RKF method for Problem 5
(a) Results for Problem 5. Since, x = 0 is the singular point for the problems and y0 = limx→0− =
1. So, the schemes do not work for the interval [−2, 0]. But schemes works for the interval [−2, 0−δ]
and δ > 1 × 1016 . I changed the problem to the following
ln(1+x)
f 0 (x) = x ,x ∈ [δ, 2]
f (δ) = 0.
The (Fig.8) gives the result for the interval [δ, 2] and δ = 1 × 1010 .
(b) Figures (Fig.10)
1.5 1.5
y
y
1 1
0.5 0.5
−2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0
x x
Figure 10: The 4th RK method and RKF method for Problem 5
Listing 2: Main function for problem1-5 with 4-th oder Runge-Kutta Method
% Script file: main1.m
% The RHS of the differential equation is defined as
% a handle function
% author:Wenqiang Feng
5 % Email: wfeng1@utk.edu
% date: Mar 8, 2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% common parameters
clc
10 clear all
n=1;
tol=1e-5;
choice=5; % The choice of the problem number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 %% the parameters for each problems
switch choice
case 1.1
% problem 11
f=@(x,y) y-x; %The right hand term
20 xinit=0;
xfinal=10;
yinit=1;%+1e-3; %The initial condition
case 1.2
Page 208 of 236
% problem 12
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 209
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 5
% problem 5
f=@(x,y) log(1+x)/x; %The right hand term
80 xinit=1e-10;
xfinal=2;
yinit=0; %The initial condition
end
85
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% computing the numberical solutions
y0=100*ones(1,n+1);
90 [x1,y1]=Runge_Kutta_4(f,xinit,yinit,xfinal,n);
k2 = h*f(t+h/4, y+k1/4);
k3 = h*f(t+3*h/8, y+3*k1/32+9*k2/32);
k4 = h*f(t+12*h/13, y+1932*k1/2197-7200*k2/2197+7296*k3/2197);
k5 = h*f(t+h, y+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
15 k6 = h*f(t+h/2, y-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
y1 = y + 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55;
y2 = y + 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
E=abs(y1-y2);
R = E/h;
20 delta = 0.84*(tol/R)ˆ(1/4);
i f E<=tol
t = t+h;
y = y1;
i = i+1;
25 fprintf(’Step %d: t = %6.4f, y = %18.15f\n’, i, t, y);
u(i)=y;
time(i)=t;
h = delta*h;
else
30 h = delta*h;
end
i f (i>1000)
disp(’the partitions excess 1000’)
break;
35 end
end
time=[t0,time];
u=[u0,u];
case 2.1
% problem 21
lambda=0;
25 f=@(x,y) lambda*y+sin(x)-lambda* cos(x); %The right hand term
xinit=0;
xfinal=5;
yinit=0; %The initial condition
case 2.2
30 % problem 22
lambda=5;
f=@(x,y) lambda*y+sin(x)-lambda* cos(x); %The right hand term
xinit=0;
xfinal=5;
35 yinit=0; %The initial condition
case 2.3
% problem 23
lambda=-5;
f=@(x,y) lambda*y+sin(x)-lambda* cos(x); %The right hand term
40 xinit=0;
xfinal=5;
yinit=0; %The initial condition
case 2.4
% problem 24
45 lambda=10;
f=@(x,y) lambda*y+sin(x)-lambda* cos(x); %The right hand term
xinit=0;
xfinal=5;
yinit=0; %The initial condition
50 case 2.5
% problem 25
lambda=-10;
f=@(x,y) lambda*y+sin(x)-lambda* cos(x); %The right hand term
xinit=0;
55 xfinal=5;
yinit=0; %The initial condition
case 3
% problem 3
f=@(x,y) 1-y/x; %The right hand term
60 xinit=2;
xfinal=20;
yinit=2; %The initial condition
case 4
% problem 4
65 f=@(x,y) sin(xˆ2); %The right hand term
xinit=0;
xfinal=10;
yinit=0; %The initial condition
70 case 5
% problem 5
f=@(x,y) log(1+x)/x; %The right hand term
xinit=1e-10;
xfinal=2; Page 212 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 213
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot
figure
90 plot(time,u,’-.’)
xlabel(’x’)
ylabel(’y’)
legend(’Runge-Kutta-Fehlberg’)
title(sprintf(’Problem.%1.1f,with step=%d,error=%1e’,choice,step,error),...
95 ’FontSize’, 14)
G Project 2 MATH572
COMPUTATIONAL ASSIGNMENT #2
MATH 572
The purpose of this assignment is to explore techniques for the numerical solution of boundary value and
initial boundary value problems and to introduce some ideas that we did not discuss in class but, nevertheless,
are quite important. You should submit the solution of at least two (2) of the following problems. Submitting
the solution to the third can be used for extra credit.
where αi , βi , γi are defined in terms of and h. Set ∈ {1, 10−1 , 10−3 , 10−6 } and compute the solution for
different values of h. What do you observe for h > ? For h ≈ ? For h < ?
• Show that
2 i
h +1
Ûi = 1, Ǔi = 2 , i = 0, N
h −1
are two linearly independent solutions of the difference equation. Find the dicrete solution U of the
problem in terms of Û and Ǔ . Using this representation, determine the relation between and h that
ensures that there are no oscillations in U . Does this coincide with your observations of the previous item?
2
Hint: Consider the sign of h +1 .
2
h −1
• Replace the centered difference approximation of the first derivative u0 by the up-wind difference u0 (xi ) ≈
h−1 (u(xi ) − u(xi−1 )). Repeat the previous two items and draw conclusions.
• Show that, using an up-wind approximation the arising matrix satisfies a discrete maximum principle.
Write a piece of code that, for a given a and f computes the finite element solution to this problem over a
mesh Th = {Ij }Nj=1 , Ij = [xj−1 , xj ] with hj = xj − xj−1 not necessarily uniform.
• Set a = 1 and choose f so that u = x3 . Compute the finite element solution on a sequence of uniform
meshes of size h = 1/N and verify the estimate
(3) ku − U kH 1 (0,1) ≤ Ch = CN −1 .
3
• Set a = 1 and f = − 4√ x
/ L2 (0, 1). This problem, however, is still well posed. Show
and notice that f ∈
this. For this case repeat the previous item. What do you observe?
• Set a(x) = 1 if 0 ≤ x < 1/π and a(x) = 2 otherwise. Choose f ≡ 1 and compute the exact solution.
Repeat the first item. What do you observe? Recall that to compute the exact solution we must include
the interface conditions: u and au0 are continuous.
The last two items show that in the case when either the right hand side or the coefficient in the equation
are not smooth, the solution does not satisfy u00 ∈ L2 (0, 1) and so the error estimate (3) cannot be obtained
with uniform meshes. Notice, also, that in both cases the solution is smooth except perhaps at very few
points, so that if we were able to handle these, problematic, points we should be able to recover (3). The
purpose of a posteriori error estimates is exactly this.
Let us recall the weak formulation of (2). Define:
Z 1 Z 1
A(v, w) = av 0 w0 , L(v) = f v,
0 0
If U is the finite element solution to (2) and v ∈ H01 (0, 1), we have
Z 1 Z 1 N Z
X
0 0
A(u − U, v) = A(u, v) − A(U, v) = L(v) − A(U, v) = fv − aU v = f v − aU 0 v 0 .
0 0 j=1 Ij
where
j(w(x)) = w(x + 0) − w(x − 0)
is the so-called jump. Let us now set v = w − Ih w, where Ih is the Lagrange interpolation operator. In this
case then v(xj ) = 0 (why?) and
Page 216 of 236
kvk 2 = kw − I wk 2 ≤ ch kw0 k 2 .
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 217
Consequently,
X
A(u − U, w − Ih w) ≤ C hj kf + (aU 0 )0 kL2 (Ij ) kw0 kL2 (Ij )
Ij ∈Th
1/2 1/2
X X
≤C h2j kf + (aU 0 )0 k2L2 (Ij ) kw0 kL2 (Ij )
Ij ∈Th Ij ∈Th
1/2
X
=C h2j kf + (aU 0 )0 k2L2 (Ij ) kw0 kL2 (0,1)
Ij ∈Th
What is the use of all this? Define rj = hj kf + (aU 0 )0 kL2 (Ij ) then, using Galerkin orthogonality we obtain
1/2
XN
1 1
ku − U k2H 1 (0,1) ≤ A(u − U, u − U ) = A(u − U, u − U − Ih (u − U )) ≤ C r2j ku − U kH 1 (0,1) .
c1 c1 j=1
In other words, we bounded the error in terms of computable and local quantities rj . This allows us to
devise an adaptive method:
• (Solve) Given Th find U .
• (Estimate) Compute the rj ’s.
• (Mark) Choose ` for which r` is maximal.
• (Refine) Construct a new mesh by bisecting I` and leaving all the other elements unchanged.
Implement this method and show that (3) is recovered.
P PN
You might also want to try choosing a set of minimal cardinality M so that j∈M r2j ≥ 12 j=1 r2j and
bisecting the cells Ij with j ∈ M.
TTH 12:40pm
Wenqiang Feng
Contents
The Convection Diffusion Equation 3
Problem 1 3
Problem 2 8
Heat Equation 13
Problem 3 13
Problem 1
1. The exact solution
From the problem, we know that the characteristic function is
−λ2 + λ = 0.
So, λ = 0, 1 . Therefore, the general solution is
1 1
u = c1 e0x + c2 e x = c1 + c2 e x .
By using the boundary conditions, we get the solution is
1 1 1
u(x) = 1 − 1 + 1 e x.
1 − e 1 − e
And u(x) is monotone.
2. Central Finite difference scheme
I consider the following partition for finite difference method:
0 1
x0 x1 xN −1 xN
3. when i = N − 1
UN −2 − 2UN −1 + UN UN − UN −2
− + = 0,
h2 2h
i.e.
1 2 1
− + UN −2 + 2 UN −1 + − UN = 0.
h2 2h h 2h h2
Since UN = 0, then,
1 2
− + UN −2 + UN −1 = 0. (5)
h2 2h h2
AU = F,
where
2 1
−
−
h2
1
2h
2
h2
1
h2 + 2h h2 − h2
2h
.. .. ..
A=
. . . ,
− h2 + 2h1 2 1
−
h2
1
2h
2
h2
− h2 + 2h h2
1
U1 h2 + 2h
..
U2 .
U = .. ,F =
0
.
.
..
UN −2 .
UN −1
Table 1: l∞ norms for the Central Difference Method with = {1, 10−1 , 10−3 , 10−6 }
(a) when h < the scheme is convergent with optimal convergence order (Figure.2), i.e.
Page 221 of 236
ku − uh kl∞ ≈ 0.01h1.9992 ,
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 222
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
(b) when h ≈ the scheme is convergent with optimal convergence order (Figure.2), i.e.
ku − uh kl∞ ≈ 3.201h2.0072 ,
(c) when h > the scheme is not stable and the solution has oscillation.
l i n e ar r e gr e ssi on f or ǫ = 1
l i n e ar r e gr e ssi on f or ǫ = 0. 1
log(error)
log(h)
u = c1 Û + c2 Ǔ
is also solution to 1. We also need this solution to satisfy the boundary conditions, so
2
u = c1 + c2 2 h +1
=1
h −1
2 N
u = c + c 2 h +1
= 0.
1 2 −1 h
so
(2 + h)N (2 − h)N
c1 = − , c2 = .
(2 + h)(2 − h)N −1 − (2 + h)N (2 + h)(2 − h)N −1 − (2 + h)N
So
U0 − 2U1 + U2 U1 − U0
− 2
+ = 0,
h h
i.e.
1 2 1
− + U0 + + U1 − U2 = 0.
h2 h h2 h h2
Since, U0 = 1, so we get
2 1 1
2
+ U1 − U2 = + . (6)
h h h2 h 2 h
(c) when i = N − 1
UN −2 − 2UN −1 + UN UN −1 − UN −2
− 2
+ = 0,
h h
i.e.
1 2 of 1236
− + UNPage
−2 + 223 + UN −1 − 2 UN = 0.
h2 h h2 h h
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 224
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
Since UN = 0, then,
1 2 1
− 2
+ UN −2 + 2
+ UN −1 = 0. (8)
h h h h
AU = F,
where
2 1
h2 + h − h2
− 1 2 1
h2 + h h2 + h − h2
.. .. ..
A=
. . . ,
− h2 + h1
2 1
h2 + h − h2
1 2
− h2 + h h2+ h1
1
U1 h2 + h
..
U2 .
U = .. ,F =
0
.
.
..
UN −2 .
UN −1
Table 2: l∞ norms for the Up-wind Difference Method with = {1, 10−1 , 10−3 , 10−6 }
(a) when h < the scheme is convergent with optimal convergence order (Figure.2), i.e.
ku − uh kl∞ ≈ 0.0471h0.946 ,
(b) when h ≈ the scheme is convergent, but the convergence order is not optimal (Figure.2), i.e.
ku − uh kl∞ ≈ 0.4398h0.6852 ,
(c) when h > the scheme is convergent, and the solution has no oscillation.
Page 224 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 225
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 1 (continued)
l i n e ar r e gr e ssi on f or ǫ = 1
l i n e ar r e gr e ssi on f or ǫ = 0. 1
log(error)
log(h)
Lemma 0.1 Let A = tridiag{ai , bi , ci }ni=1 ∈ Rn×n be a tridiagional matrix with the properties that
bi > 0, ai , ci ≤ 0, ai + bi + ci = 0.
Then the following maximum principle holds: If u ∈ Rn is such that (Au)i=2,··· ,n−1 ≤ 0, then ui ≤
max{u1 , un }.
From the Up-wind Difference scheme, we get that a1 = 0, ai = − h2 + h1 , i = 2, · · · , n, bi =
2 1
h2 + h , i = 1, · · · , n and ci = − h2 , i = 1, · · · , n − 1, moreover (Au)i=2,··· ,n−1 = 0. Therefore,
bi > 0, ai , ci ≤ 0, ai + bi + ci = 0.
Since (Au)i=2,··· ,n−1 = 0, so the corresponding matrix arising from the up-wind scheme satisfies the
discrete maximum principle(Lemma 0.1).
Problem 2
1. Partition
I consider the following partition for finite element method:
0 1
x1 x2 xN −1 xN
I will use the linear basis function, i.e. for each element I = [xi , xi+1 ]
( xi+1 −x
φ1 (x) = xi+1 −xi
φI (x) = x−xi
φ2 (x) = xi+1 −xi .
3. Weak Formula
Multiplying the testing function v ∈ H01 to both side of the problem, then integrating by part we get
the following weak formula
Z 1 Z 1
0 0
a(x)u v dx = f vdx.
0 0
4. Approximate Problem
The approximate problem is to find uh ∈ H 1 , s.t
where
Z 1 Z 1
a(uh , vh ) = a(x)u0h vh0 dx and f (vh ) = f vh dx.
0 0
h Nnodes ku − uh kL2 |u − uh |H 1
1/4 5 1.791646 × 10−2 2.480392 × 10−1
1/8 9 4.502711 × 10−3 1.247556 × 10−1
1/16 17 1.127148 × 10−3 6.246947 × 10−2
1/32 33 2.818787 × 10−4 3.124619 × 10−2
1/64 65 7.047542 × 10−5 1.562452 × 10−2
1/128 128 1.761921 × 10−5 7.812440 × 10−3
1/256 257 4.404826 × 10−6 3.906243 × 10−3
Using linear regression (Figure.5), we can also see that the errors in Table.4 obey
ku − uh kL2 ≈ 0.2870h1.9987 ,
ku − uh kH 1 ≈ 0.9935h0.9986 .
These linear regressions indicate that the finite element method for this problem can converge in
the optimal rates, which are second order in L2 norm and first order in H 1 norm.
l i n e ar r e gr e ssi on f or L 2 n or m e r r or
l i n e ar r e gr e ssi on f or H 1 n or m e r r or
L 2 n or m e r r or
H 1 n or m e r r or
log(error)
log(h)
3 3
(b) Problem: a(x)=1, ue = x 2 and f = − 4√ x
.
h Nnodes ku − uh kL2 |u − uh |H 1
1/4 5 7.625472 × 10−3 1.022294 × 10−1
1/8 9 2.029299 × 10−3 5.585353 × 10−2
1/16 17 5.324774 × 10−4 3.011300 × 10−2
1/32 33 1.378846 × 10−4 1.607571 × 10−2
1/64 65 3.523180 × 10−5 8.517032 × 10−3
1/128 128 8.876332 × 10−6 4.485323 × 10−3
1/256 257 2.203920 × 10−6 2.350599 × 10−3
Using linear regression (Figure.6), we can also see that the errors in Table.4 obey
ku − uh kL2 ≈ 0.1193h1.9593 ,
ku − uh kH 1 ≈ 0.3682h0.9081 .
These linear regressions indicate that the finite element method for this problem can converge,
but not in the optimal rates.
l i n e ar r e gr e ssi on f or L 2 n or m e r r or
l i n e ar r e gr e ssi on f or H 1 n or m e r r or
L 2 n or m e r r or
H 1 n or m e r r or
log(error)
log(h)
We can not use the uniform mesh to compute this problem. Since if we can use the uniform mesh,
then π1 should be the node point, that is to say
1 1
nh = n = ,
Nelem π
i.e.
nπ = Nelem , n, Nelem ∈ Z.
Using linear regression, we can also see that the errors (Figure.7) in Table.5 obey
−1.0798
ku − uh kH 1 ≈ 0.6454Nelem .
These linear regressions indicate that the adaptive finite element method for this problem can
converge in the optimal rates, which is first order in H 1 norm.
0.018
L 2 n or m e r r or
H 1 n or m e r r or
0.016
0.014
0.012
0.01
0.008
0.006
0.004
0.002
0
1 2 3 4 5 6 7 8 9 10
Using linear regression, we can also see that the errors (Figure.8) in Table.6 obey
−0.9706
ku − uh kH 1 ≈ 0.1825Nelem .
These linear regressions indicate that the adaptive finite element method for this problem can
converge in the optimal rates, which is first order in H 1 norm.
0.1
H 1 n or m e r r or
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
1 2 3 4 5 6 7 8 9 10
Heat Equation
Problem 3
1. Partition
I consider the following partition for finite element method:
0 1
x0 x1 xN −1 xN
= k
θµUi−1 − (2θµ − 1)Uik + θµUi+1
k
+ θτ fik + (1 − θ)τ fik+1 .
Since U (0) = U (1) = 0So, the θ-scheme can be written as the following matrix form
AU k+1 = BU k + F,
where
2(1 − θ)µ + 1 −(1 − θ)µ
−(1 − θ)µ 2(1 − θ)µ + 1 −(1 − θ)µ
.. .. ..
A=
. . . ,
−(1 − θ)µ 2(1 − θ)µ + 1 −(1 − θ)µ
−(1 − θ)µ 2(1 − θ)µ + 1
−(2θµ − 1) θµ
θµ −(2θµ − 1) θµ
.. .. ..
B=
. . . ,
θµ −(2θµ − 1) θµ
θµ −(2θµ − 1)
U k+1 (x1 ) U k (x1 )
U k+1
(x2 ) k
U (x2 )
.. ..
U k+1 =
.
, Uk =
.
,
U k+1 (xN −2 ) U k (xN −2 )
U k+1Page
(xN −1231) of 236 U k (xN −1 )
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 232
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 3 (continued)
f k (x1 ) f k+1 (x1 )
.. ..
. .
F = θτ f k (xi ) + (1 − θ)τ f k+1 (xi ) .
.. ..
. .
f k (xN −1 ) f k+1 (xN −1 )
From the Table(7)-(9), we can conclude that when µ < 0.5, Implicit Euler method, Explicit
Euler method and Crank-Nicolson method are convergent with optimal order in spacial, which
are second order in L∞ , L2 norm and first order in H 1 norm.
√
(b) Numerical results for θ-Method for τPage
= 232
h of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 233
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 3 (continued)
√
Table 10: Error norms for the Implicit Euler method with τ = h
√
Table 11: Error norms for the Explicit Euler method with τ = h
√
Table 12: Error norms for the Crank-Nicolson method with τ = h
From the Table(10)-(12), we can conclude that Implicit Euler method and Crank-Nicolson method
are unconditional stable, while when µ > 12 Explicit Euler method is not stable.
(c) Numerical results for θ-Method for τ = h
Table 14: Error norms for the Explicit Euler method with τ = h
From the Table(13)-(15), we can conclude that Implicit Euler method and Crank-Nicolson method
are unconditional stable, while when µ > 12 Explicit Euler method is not stable.
(d) Numerical results for θ-Method for τ = h2
Table 16: Error norms for the Implicit Euler method with τ = h2
From the Table(16)-(18), we can conclude that Implicit Euler method and Crank-Nicolson method
are unconditional stable, while when µ > 21 Explicit Euler method is not stable. Moreover, by
using linear regression (Figure.10) for Implicit Euler method errors, we can see that the errors in
Table.16 obey
ku − uh kL2 ≈ 0.9435h2.0580 ,
ku − uh kH 1 ≈ 7.2858h1.0137 .
These linear regressions indicate that the finite element method for this problem can converge in
the optimal rates, which are second order in L2 norm and first order in H 1 norm.
l i n e ar r e gr e ssi on f or L 2 n or m e r r or
l i n e ar r e gr e ssi on f or H 1 n or m e r r or
L 2 n or m e r r or
H 1 n or m e r r or
log(error)
log(h)
Figure 10: linear regression for L2 and H 1 norm errors of Implicit Euler method with τ = h2
Similarly, by using linear regression (Figure.11) for Crank-Nicolson Method, we can also see that
the errors in Table.18 obey
ku − uh kL2 ≈ 0.9382h2.0574 ,
ku − uh kH 1 ≈ 7.2445h1.0131 .
These linear regressions indicate that the finite element method for this problem can converge in
the optimal rates, which are second order in L2 norm and first order in H 1 norm.
Page 235 of 236
Wenqiang Feng Prelim Exam note for Numerical Analysis Page 236
Wenqiang Feng MATH 572 ( TTH 12:40pm): Computational Assignment #2 Problem 3 (continued)
l i n e ar r e gr e ssi on f or L 2 n or m e r r or
l i n e ar r e gr e ssi on f or H 1 n or m e r r or
L 2 n or m e r r or
H 1 n or m e r r or
log(error)
log(h)
Figure 11: linear regression for L2 and H 1 norm errors of Crank-Nicolson method with τ = h2