ENGR20005
Numerical Methods in Engineering
Workshop 4 Solutions
3.5 (a) The system of equations may be written as
4 −1 −2 x 1
1 −2 3y=4 (1)
111z1 (b) Writing Eq. (1) in an augmented matrix, we have
4 −1 −2 1
1 −2 3 4 (2)
1111
Using Gaussian elimination, Eq. (2) may be put into row-echelon form by
14 r 1
4 −1 −2 1r2−14r11 −1 −1 1
004 1−1−1 1
1r3 424
4 15 −→0 1 −2 −7
6 7
Then using back substitution, the solution to Eq. (1) is given by
x 4 7
7 y = −3
6 7
1
z
r3−1r1 424 4 7715
1−234−→0−4 2 4
1111 0 −47r21−1−1 1
533 424
r3+5r2 4 2 4
7 15
−→0 1 −2−7
001
24 7
3.6 (a)
(c) Using Crout’s method, the [L] and [U] matrices are given by
4 0 0 [L]=1−74 0
1 45 4
1 −14 −21 [U]=0 1 −2
001 (d) To solve Eq. (1), we first solve the system
[L]{r} = {c} using forward substitution, which gives
r1 1 4
7 r2 = −15
6 7
[U]{x} = {r} x 4
r3
Then, using back substitution, we may find the vector {x} using
which gives
and hence
Decomposing the matrix
r3
x 65
into the [M] and [N] matrices for
2
7 7
y = −3
6 7
r1 1 2
7 r2 = −22
z
(e) If we now replace the right hand side of Eq. (1) with {c} = {2, 6, 7}, we have
28 14
y = −29
73 28
20 −3 [A]= −10 2
z
73 28
The Point Jacobi method:
The Gauss–Seidel method:
20 0 [M]= 0 2
0 3 [N]= 10 0
20 0 [M]= −10 2
0 3 [N]= 0 0
3.7 (a)
The Point Jacobi method:
k123 1 0.1500 0.8250 0.1875
1 5.5000 1.2500 4.6250 The Gauss–Seidel method:
k123 1 0.1500 0.1875 0.2156
1 1.2500 1.4375 1.5781
If the matrix [A] has a Cholesky decomposition, the we only have to store a single
(b) Performing three iterates for each method, we have
matrix [L]. Whereas for an LU decomposition, we need to store both [L] and [U]. (b) Expanding [L][L]T and equating it with [A] gives
l2 ll ll
11
[L][L]T = l l 21 11
a11 a21 a31 = a21 a22 a32
a31 a32 a33 (c) Equating each term of the matrix we have
l2 11
=a = a21
(3) (4) (5) (6) (7) (8)
l l
l2 +l2 21 22
l l +l32l
l l +l32l
31 21 22
l2 +l2 +l2
31 11
31 21
22
31
32
33
21 11
31 11
l21 l11 l2 +l2
=a = a31
= a32
21 22 l31 l11 l31l21 + l32l22
22
l2 +l2 +l2 31 32 33
=a
33
3
11
Solving Eq. (3) to (8) sequentially, we have
l11 = √a11 (9)
l21 = a21 (10) l11
l22 = a22 − l21 (11) l31 = a31/l11 (12)
l32 = (a32 − l31l21)/l22 (13) l33 = a33 − l321 − l322 (14)
which gives the matrix
√a11 0 0
[L] = a21/l11 a22 − l21 0 (15)
a31/l11 (a32 − l31l21)/l22 a33 − l321 − l322
Eq. (15) shows us that the diagonal terms are formed by taking the square root
of a2ii minus the sum of the squares in the ith row, thus we have
j−1
l j j = ± a j j − l 2
As for the off diagonal terms, we see that the terms are given by aij minus the sum of the product of the ith and jth terms
1 j−1 lij = l aij −likljk
jj k=1
Note that from Eq. (15), we see that we can either find [L] by starting from the top left corner and moving rowwise, which is called the Cholesky–Banachiewicz algorithm, or columnwise, the Cholesky–Crout algorithm.
(d),(e) An implementation of the Cholesky–Crout algorithm is given in Q4 7 sol.m.
4
jk k=1