Download Review Practice Problems on Linear Algebra | MATH 371 and more Study notes Mathematics in PDF only on Docsity!
1
section 3.0 : review of linear algebra
a 11 x 1 + a 12 x 2 + · · · + a (^1) n x (^) n = b (^1) a 21 x 1 + a 22 x 2 + · · · + a (^2) n x (^) n = b (^2) ... a (^) n 1 x 1 + a (^) n 2 x 2 + · · · + a (^) nn x (^) n = b (^) n
: system of linear equations for x 1 ,... , x (^) n
We can write the system in 3 other forms.
∑n j=
a (^) ij x (^) j = b (^) i , i = 1 : n , i : row index , j : column index
a 11 a 12 · · · a (^1) n a 21 a 22 · · · a (^2) n ... ... ... a (^) m 1 a (^) m 2 · · · a (^) mn
x (^1) x (^2) ... x (^) n
^ =
b (^1) b (^2) ... b (^) m
- Ax = b , A = (aij )
basic problem : Given A, b, find x.
solution : x = b/A : no, but x = b\A does work in Matlab
thm : The following conditions are equivalent.
- The equation Ax = b has a unique solution for any vector b.
- A is invertible, i.e. there exists A−^1 such that AA −^1 = I
- det A "= 0
- The equation Ax = 0 has the unique solution x = 0.
- The columns of A are linearly independent.
- The eigenvalues of A are nonzero.
note
- If A is invertible, then x = A−^1 b (because then Ax = A(A−^1 b) = (AA −^1 )b = Ib = b), but this is not the best way to compute x in practice.
- There are two types of methods for solving linear systems, direct methods and iterative methods. We will begin with direct methods.
- Thurs 9/18 2
section 3.1 : Gaussian elimination for solving Ax = b
First consider the special case in which A is upper triangular.
a 11 x 1 + a 12 x 2 + · · · + a (^1) n x (^) n = b (^1)
a 22 x 2 + · · · + a (^2) n x (^) n = b (^2)
... a (^) n− 1 ,n− 1 x (^) n− 1 + a (^) n− 1 ,n x (^) n = b (^) n− 1 a (^) nn x (^) n = b (^) n
Then x (^) n = b (^) n /ann
x (^) n− 1 = (b (^) n− 1 − a (^) n− 1 ,n x (^) n )/an− 1 ,n− 1 ... x 1 = (b 1 − (a 12 x 2 + · · · + a (^1) n x (^) n ))/a 11
back substitution
- x (^) n = b (^) n /ann
- for i = n − 1 : −1 : 1
- sum = b (^) i
- for j = i + 1 : n
- sum = sum − a (^) ij · x (^) j
- x (^) i = sum/aii
operation count
divisions = n
mults = # adds = 12 n(n − 1) = 12 n 2 − 12 n ∼ 12 n 2 for large n
pf
mults = 1 + 2 + · · · + (n − 1) =
n∑− 1 i=
i = S
2 S =
n∑− 1 i=
i +
n∑− 1 i=
(n − i) =
n∑− 1 i=
(i + (n − i)) =
n∑− 1 i=
n = n(n − 1)
⇒ S = 12 n(n − 1) ok
Hence the leading order term in the operation count for back substitution is n 2.
note : Similar considerations hold if A is lower triangular.
4
ex
2 x 1 − x 2 = 1
−x 1 + 2x 2 − x 3 = 0
−x 2 + 2x 3 = 1
m^21 =^ −^1 /^2 m 31 = 0
m 32 = − 1 /(3/2) = − 2 / 3
x 3 = 1 , x 2 = ( 12 − (−1) · 1)/ 32 = 1 , x 1 = (1 − (−1) · 1)/2 = 1 check : ok
general n × n case
reduction to upper triangular form
- for k = 1 : n − 1
- for i = k + 1 : n
- m (^) ik = a (^) ik /akk
- for j = k + 1 : n
- a (^) ij = a (^) ij − m (^) ik · a (^) kj
- b (^) i = b (^) i − m (^) ik · b (^) k
note
The element a (^) kk in step k is called a pivot. In the previous example the pivots are 2, 32 , 43 (these are the diagonal elements in the last step).
operation count
The leading order term comes from line 5.
k = 1 ⇒ 2(n − 1)^2 ops k = 2 ⇒ 2(n − 2)^2 ops ... k = n − 2 ⇒ 2 · 22 ops k = n − 1 ⇒ 2 · 12 ops
n∑− 1 k=
k 2 = 2 · (n − 1)n(2n − 1) 6
Hence the operation count for Gaussian elimination is 23 n 3.
- Tues 9/23 5
note ∑^ n
k=
k =
n(n + 1) 2
∑n k=
k 2 =
n(n + 1)(2n + 1) 6
pf
- already done
- n 3 = n 3 − (n − 1)^3 + (n − 1)^3 + · · · − 23 + 2^3 − 13 + 1^3 =
∑n k=
(k 3 − (k − 1)^3 )
k 3 − (k − 1)^3 = k 3 − (k 3 − 3 k 2 + 3k − 1) = 3k 3 − 3 k + 1
n 3 =
∑n k=
(3k 2 − 3 k + 1) = 3
∑n k=
k 2 − 3
∑n k=
k +
∑n k=
1 = 3S − 3
n(n + 1) 2
S =
∑n k=
k 2 =
( n 3 +
n(n + 1) − n
)
n
( n 2 +
(n + 1) − 1
)
n
( n 2 +
n +
)
n ·
(2n 2 + 3n + 1) =
n ·
(2n + 1)(n + 1) ok
ex : electric circuit for charging a car battery (page 129, problem 13)
12 V
100 V
15 Ω^ DC generator
Apply Kirchoff’s voltage law and current law to determine the currents.
- The sum of the voltage changes around any closed loop is zero.
⇒ 4 I 2 + 12 − 15 I 3 = 0 , 15 I 3 − 100 + 10I 1 = 0 (using Ohm’s law V = IR)
- The sum of the currents flowing into a junction equals the sum flowing out.
⇒ I 1 = I 2 + I 3
I 1
I 2
I 3
=
Hence we can’t apply Gaussian elimination because the 1st pivot is zero.
- Thurs 9/25 7
section 3.3 : vector and matrix norms
def : A vector norm ||x|| is a function satisfying the following properties.
- ||x|| ≥ 0 and ||x|| = 0 ⇔ x = 0
- ||αx|| = |α| · ||x||
- ||x + y|| ≤ ||x|| + ||y|| : triangle inequality
ex
||x|| 2 = (x T^ · x)^1 /^2 =
∑^ n i=
x (^2) i
1 / 2 : Euclidean length
||x||∞ = max{|x (^) i | : i = 1,... , n} , pf...
ex : x =
( 1 2
) ⇒ ||x|| 2 =
5 , ||x||∞ = 2
note : Given a vector norm ||x||, we can define a matrix norm by ||A|| = max x"=
||Ax|| ||x||
The following properties are satisfied.
- ||A|| ≥ 0 and ||A|| = 0 ⇔ A = 0
- ||αA|| = |α| · ||A||
- ||A + B|| ≤ ||A|| + ||B||
- ||Ax|| ≤ ||A|| · ||x||
- ||AB|| ≤ ||A|| · ||B|| , pf :...
thm : ||A||∞ = max x"=
||Ax||∞ ||x||∞ = max i
∑ j
|a (^) ij | : max row sum
ex : A =
( 1 2 0 2
) ⇒ ||A||∞ = 3
x =
( 1 0
) ⇒ Ax =
( 1 0
) ⇒
||Ax||∞ ||x||∞
x =
( 0 1
) ⇒ Ax =
( 2 2
) ⇒
||Ax||∞ ||x||∞
x =
( 1 1
) ⇒ Ax =
( 3 2
) ⇒ ||Ax||∞ ||x||∞
= 3 : largest possible ratio
thm : ||A|| 2 = max x"=
||Ax|| 2 ||x|| 2
= max{
λ : λ is an eigenvalue of AT A}
ex : A T A =
( 1 0 2 2
) ( 1 2 0 2
)
( 1 2 2 8
) ⇒ ||A|| 2 = 2. 9208 (Matlab)
8
section 3.4 : error analysis
Ax = b , A : invertible
x : exact solution , x = A−^1 b
x ˜ : approximate solution , e = x − x˜ : error , r = b − Ax˜ : residual
note : Ae = r , pf : Ae = A(x − x˜) = Ax − Ax˜ = b − Ax˜ = r ok
Then e = 0 if and only if r = 0. However, if ||r|| is small, there’s no guarantee that ||e|| is also small.
ex (
- 01 0. 99 2
- 99 1. 01 2
) ⇒ x =
( 1 1
)
x ˜ 1 =
(
- 01
- 01
) ⇒ e 1 =
( − 0. 01 − 0. 01
) , ||e 1 ||∞ = 0. 01 , r 1 =
( − 0. 02 − 0. 02
) , ||r 1 ||∞ = 0. 02
x˜ 2 =
( 2 0
) ⇒ e 2 =
( − 1 1
) , ||e 2 ||∞ = 1 , r 2 =
( − 0. 02
- 02
) , ||r 2 ||∞ = 0. 02
Hence we see that ||r|| is small in both cases, while ||e|| is small in case 1 and 100 times larger in case 2. How large can ||e|| be?
thm :
||e|| ||x||
≤ κ(A)
||r|| ||b||
, where κ(A) = ||A|| · ||A −^1 || : condition number
ex
A =
(
- 01 0. 99
- 99 1. 01
) ⇒ ||A||∞ = 2
A−^1 =
( a b c d
)− 1
ad − bc
( d −b −c a
)
(
- 01 − 0. 99 − 0. 99 1. 01
)
(
- 25 − 24. 75 − 24. 75 25. 25
) ⇒ ||A −^1 ||∞ = 50 ⇒ κ∞ (A) = 100 ok
pf
- ||b|| = ||Ax|| ≤ ||A|| · ||x|| ⇒
||x||
||A||
||b||
Ae = r ⇒ e = A−^1 r ⇒ ||e|| = ||A −^1 r|| ≤ ||A−^1 || · ||r||
||e|| ||x||
= ||e|| ·
||x||
≤ ||A−^1 || · ||r|| ·
||A||
||b||
= κ(A) ·
||r|| ||b||
ok
10
section 3.5 : LU factorization
matrix form of Gaussian elimination
We consider the 3 × 3 case (but the general n × n case is similar).
a 11 a 12 a (^13) a 21 a 22 a (^23) a 31 a 32 a (^33)
step 1 : eliminate variable x 1 from eqs. 2 and 3
m 21 =
a (^21) a (^11) , m 31 =
a (^31) a (^11)
−m 21 1 0 −m 31 0 1
a 11 a 12 a (^13) a 21 a 22 a (^23) a 31 a 32 a (^33)
=
a 11 a 12 a (^13) 0 a 22 a (^23) 0 a 32 a (^33)
step 2 : eliminate variable x 2 from eq. 3
m 32 =
a (^32) a (^22)
0 −m 32 1
a 11 a 12 a (^13) 0 a 22 a (^23) 0 a 32 a (^33)
=
a 11 a 12 a (^13) 0 a 22 a (^23) 0 0 a (^33)
=^ U^ :^ upper triangular
⇒ E 2 E 1 A = U ⇒ A = E 1 − 1 E 2 − 1 U
E 1 =
−m 21 1 0 −m 31 0 1
, E^1 − 1 =
m 21 1 0 m 31 0 1
, check :^ E 1 E^ − 1 1 =
E 2 =
0 −m 32 1
, E^ − 1 2 =
0 m 32 1
,^ check^...
E 1 − 1 E 2 − 1 =
m 21 1 0 m 31 0 1
0 m 32 1
=
m 21 1 0 m 31 m 32 1
=^ L^ :^ lower triangular
final result : A = LU
- Thurs 10/2 11
ex
→
→
m 21 = −^12 m 32 = 3 −/^12 = − (^23) m 31 = 02 = 0
check :
LU =
=
=^ A^ ok
note
To solve Ax = b.
- factor A = LU
- solve Ly = b by forward substitution
- solve U x = y by back substitution
check : Ax = LU x = Ly = b ok
ex
A =
, b^ =
Using Gaussian elimination before we found x 1 = x 2 = x 3 = 1, but now we’ll use LU factorization.
Ly = b ⇒
y (^1) y (^2) y (^3)
=
⇒
y (^1) y (^2) y (^3)
=
1 (^24) 3
U x = y ⇒
x (^1) x (^2) x (^3)
=
1 (^24) 3
⇒
x (^1) x (^2) x (^3)
=
ok
question : So what’s the point of LU factorization?
answer : Some applications require solving Ax = b for a given matrix A and a sequence of different vectors b (e.g. in a time-dependent problem). Once the LU factorization of A is known, we can apply forward and back substitution to the sequence of vectors b - it’s not necessary to repeat the LU factorization.
- Tues 10/7 13
section 8.1 : 2-point boundary value problem
Find y(x) on 0 ≤ x ≤ 1 satisfying the differential equation −y ′′^ + q(x)y = r(x) subject to boundary conditions y(0) = α, y(1) = β, where q(x), r(x) are given. The equation models (for example) a steady state convection-reaction-diffusion system, where y(x) represents a velocity or temperature profile.
finite-difference scheme choose an integer n ≥ 1 and define h = (^) n+1^1 : mesh size
set x (^) i = ih for i = 0, 1 ,... , n + 1 : mesh points (note : x 0 = 0 , x (^) n+1 = 1)
y(x (^) i ) = y (^) i : exact solution , q (^) i = q(x (^) i ) , r (^) i = r(x (^) i )
recall : D+ y (^) i = y (^) i+1 − y (^) i h
, D (^) − y (^) i = y (^) i − y (^) i− 1 h
D+ D− y (^) i = D+ (D− y (^) i ) = D+
( y (^) i − y (^) i− 1 h
)
h
(D (^) + y (^) i − D+ y (^) i− 1 )
h
( y (^) i+1 − y (^) i h
( y (^) i − y (^) i− 1 h
))
y (^) i+1 − 2 y (^) i + y (^) i− 1 h^2 ≈ y ′′^ (x (^) i )
question : How accurate is the approximation?
y (^) i+1 = y(x (^) i+1 ) = y(x (^) i + h) : expand in a Taylor series about x = x (^) i
y (^) i+1 = y (^) i + hy (^) i′ + h^ 2 2 y^
′′ i +^ h 3 3! y^
′′′ i +^ h 4 4! y^
(4) i +^ h^
5 5! y^
(5) i +^ O(h^6 )
y (^) i− 1 = y (^) i − hy (^) i′ + h^ 2 2 y^
′′ i −^
h 3 3! y^
′′′ i +^
h 4 4! y^
(4) i −^ h 5 5! y^
(5) i +^ O(h^6 )
D+ D− y (^) i = y (^) i+1 − 2 y (^) i + y (^) i− 1 h^2
= y (^) i′′ + h^2 12
y (^) i(4) + O(h^4 ) : 2nd order accurate ︸ ︷︷ ︸ ↑ ︸ ︷︷ ︸ discrete exact discretization approximation value error
wi : numerical solution , wi ≈ y (^) i , w 0 = α , wn+1 = β
−
( (^) w i+1 −^2 wi +^ wi− 1 h^2
)
- q (^) i wi = r (^) i , i = 1,... , n 1 h 2
( −wi+1 +
( 2 + q (^) i h^2
) wi − wi− 1
) = r (^) i
i = 1 ⇒ (^) h^1
( −w 2 +
( 2 + q 1 h^2
) w 1 − α
) = r (^1)
i = n ⇒ (^) h^1
( −β +
( 2 + q (^) n h^2
) wn − wn− 1
) = r (^) n
14
h^2
2 + q 1 h^2 − 1 − 1 2 + q 2 h^2 − 1
......... ......... ......... − 1 2 + q (^) n− 1 h^2 − 1 − 1 2 + q (^) n h^2
w 1 w 2 ... ... ... wn− 1 wn
r 1 + α/h^2 r (^2) ... ... ... r (^) n− 1 r (^) n + β/h^2
A (^) h wh = r (^) h , Ah : tridiagonal , symmetric questions
- Is Ah invertible for all h, q(x), r(x)?
- Can wh be computed efficiently?
- Does wh → y (^) h as h → 0, i.e. does the scheme converge? If so, what is the order of accuracy?
LU factorization for a tridiagonal system (Thomas algorithm)
b 1 c (^1) a 2 b 2 c (^2)
......... ...... (^) c n− 1 a (^) n b (^) n
l 2 1
...... ...... ln 1
u 1 c (^1) u 2 c (^2)
...... ... (^) c n− 1 u (^) n
find L, U b 1 = u 1 ⇒ u 1 = b (^1) ak = lk u (^) k− 1 ⇒ lk = a (^) k /uk− 1 b (^) k = lk c (^) k− 1 + u (^) k ⇒ u (^) k = b (^) k − lk c (^) k− 1
} for k = 2 : n
solve Lz = r z 1 = r (^1) lk zk− 1 + zk = r (^) k ⇒ zk = r (^) k − lk zk− 1 for k = 2 : n
solve U w = z u (^) n wn = zn ⇒ wn = zn /un u (^) k wk + c (^) k wk+1 = zk ⇒ wk = (zk − c (^) k wk+1 )/uk for k = n − 1 : − 1 : 1
note operation count = O(n) memory = O(n) if vectors are used instead of full matrices
- Thurs 10/9 16
review of two important properties of norms
first recall : ||x||∞ = max{|x 1 |,... , |x (^) n |}
||A||∞ = max x"=
||Ax||∞ ||x||∞
= max i
∑ j
|a (^) ij |
property 1 : ||Ax||∞ ≤ ||A||∞ · ||x||∞
property 2 : ||AB||∞ ≤ ||A||∞ · ||B||∞
pf
- the proof is a bit abstract, so instead we’ll give an example (which shows the idea of the proof)
x =
( 1 − 2
) , A =
( 3 − 4 1 0
)
Ax =
( 3 − 4 1 0
) ( 1 − 2
)
( (3)(1) + (−4)(−2) (1)(1) + (0)(−2)
)
( 11 1
)
||x||∞ = 2 , ||Ax||∞ = 11 , ||A||∞ = max{| 3 | + |− 4 |, | 1 | + | 0 |} = 7
||Ax||∞ = 11 ≤ 14 = 7 · 2 = ||A||∞ · ||x||∞ ok
physical meaning of property 1
If x is the input to a system and Ax is the output, then the norm of the output is bounded by ||A||∞ times the norm of the input.
- here we give the actual proof
||AB||∞ = max x"=
||ABx||∞ ||x||∞
: by definition
≤ max x"=
||A||∞ · ||Bx||∞ ||x||∞
: by property 1
≤ max x"=
||A||∞ · ||B||∞ · ||x||∞ ||x||∞ : by property 1 again
= ||A||∞ · ||B||∞ ok
note
it follows that ||A 2 ||∞ ≤ ||A||^2 ∞ , ||A k^ ||∞ ≤ ||A||k ∞ for any k ≥ 1
17
section 3.8 : iterative methods
Gaussian elimination is an example of a direct method for solving Ax = b, in the sense that the exact solution is obtained after a finite number of steps. In practice, the O(n 3 ) operation count is a serious obstacle when n is large (and storage can be an issue too). Here we consider an alternative class of methods called iterative methods which generate a sequence of approximate solutions x (^) k such that limk→∞ x (^) k = x. As we shall see, iterative methods have some advantages over direct methods.
Ax = b ⇔ x = Bx + c : equivalent linear system
x (^) k+1 = Bx (^) k + c : fixed-point iteration
B : iteration matrix
Jacobi method
A = L + D + U : splitting
D = diag(a 11 ,... , a (^) nn ) , assume a (^) ii "= 0, i = 1 : n
L =
a 21 0 ......... ......... a (^) n 1 · · · · · · a (^) n,n− 1 0
, U =
0 a 12 · · · · · · a (^1) n 0...^ ...
...... ... ... (^) a n− 1 ,n 0
Ax = b ⇔ (L + D + U )x = b
⇔ Dx = −(L + U )x + b
⇔ x = −D −^1 (L + U )x + D −^1 b , B (^) J = −D −^1 (L + U )
Dx (^) k+1 = −(L + U )x (^) k + b
component form
a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 ⇒ a 11 x ( 1 k +1)= −
( a 12 x ( 2 k )+ a 13 x ( 3 k)
)
a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 ⇒ a 22 x ( 2 k +1)= −
( a 21 x ( 1 k )+ a 23 x ( 3 k)
)
a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 ⇒ a 33 x ( 3 k +1)= −
( a 31 x ( 1 k )+ a 32 x ( 2 k)
)
19
Gauss-Seidel method
A = L + D + U : as before
Ax = b ⇔ (L + D + U )x = b
⇔ (L + D)x = −U x + b
⇔ x = −(L + D)−^1 U x + (L + D)−^1 b , B (^) GS = −(L + D)−^1 U
(L + D)x (^) k+1 = −U x (^) k + b : solve by forward elimination
component form
a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 ⇒ a 11 x ( 1 k +1)= −
( a 12 x ( 2 k )+ a 13 x ( 3 k)
)
a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 ⇒ a 22 x ( 2 k +1)= −
( a 21 x ( 1 k +1)+ a 23 x ( 3 k)
)
a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 ⇒ a 33 x ( 3 k +1)= −
( a 31 x ( 1 k +1)+ a 32 x ( 2 k+1)
)
Hence x ( i k+1)is used as soon as it’s computed, in contrast with Jacobi.
ex
2 x 1 − x 2 = 1 ⇒ 2 x ( 1 k +1)= x ( 2 k )+ 1
−x 1 + 2x 2 = 1 ⇒ 2 x ( 2 k +1)= x ( 1 k +1)+ 1
k x ( 1 k ) x ( 2 k) 0 0 0 1 1/2 3/ 2 7/8 15/ 3 31/32 63/
Hence Gauss-Seidel converges faster than Jacobi.
||e 0 ||∞ = 1 , ||e 1 ||∞ = 12 , ||e 2 ||∞ = 18 , ||e 3 ||∞ = 321 ⇒ ||ek+1 ||∞ = 14 ||ek ||∞
A =
( 2 − 1 − 1 2
) ⇒ BGS = −(L + D)−^1 U = − (^14)
( 2 0 1 2
) ( 0 − 1 0 0
)
( (^0 ) (^0 )
)
||BGS ||∞ = 12
- Since ||BGS ||∞ < 1, the theorem implies that Gauss-Seidel converges.
- The error norm decreases by a factor of 14 in each step, which is less than 1 2 =^ ||BGS^ ||∞^.
- Thurs 10/23 20
summary
A =
( 2 − 1 − 1 2
)
BJ =
( (^0 ) 1 2 0
) ⇒ ||BJ ||∞ = 12 , ||ek+1 ||∞ = 12 ||ek ||∞
BGS =
( (^0 ) (^0 )
) ⇒ ||BGS ||∞ = 12 , ||ek+1 ||∞ = 14 ||ek ||∞
We see that ||ek+1 ||∞ ≤ ||B||∞ · ||ek ||∞ in both cases (as required by the previous theorem), but the bound is not sharp in the case of GS. To explain the actual behavior of the error, we need to examine the eigenvalues of the iteration matrix.
def
If Ax = λx, where x "= 0 is a vector (real or complex) and λ is a scalar (real or complex), then λ is an eigenvalue of A and x is a corresponding eigenvector.
ex
A =
( 0 1 1 0
) : permutation matrix , A
( 1 1
)
( 1 1
) , A
( 1 − 1
)
( − 1 1
)
⇒ λ = 1 is an e-value with corresponding e-vector x =
( 1 1
)
λ = − 1.................. ”................... x =
( 1 − 1
)
note
Ax = λx , x "= 0 ⇔ (A − λI)x = 0 , x "= 0 ⇔ det(A − λI) = 0
fA (λ) = det(A − λI) : characteristic polynomial of A
The e-values of A are the roots of the characteristic polynomial.
ex : fA (λ) = det(A − λI) = det
( −λ 1 1 −λ
) = λ 2 − 1 = 0 ⇒ λ = ± 1 ok
thm
If A is upper triangular, then the e-values of A are the diagonal elements.
pf
A =
a 11 · · · a (^1) n
... ... 0 ann
⇒^ A^ −^ λI^ =
a 11 − λ · · · a (^1) n
... ... 0 ann − λ
fA (λ) = det(A − λI) = (a 11 − λ) · · · (a (^) nn − λ) = 0 ⇒ λ = a (^) ii for some i ok
