Chapter 1 Numerical Linear Algebra
1.1 A Simple Lower Triangular System:
We will start by writing a function to solve a
lower triangular matrix system Lx = b.
1.1.1 A Lower Triangular Solver:
Here is a simple function to solve such a system.
LTriSol.m
1 function x = LTriSol(L,b)
2 %
3 % L is n x n Lower Triangular Matrix
4 % b is nx1 data vector
5 % Obtain x by forward substitution
6 %
7 n = length(b);
8 x = zeros(n,1);
9 for j=1:n-1
10 x(j) = b(j)/L(j,j);
11 b(j+1:n) = b(j+1:n) - x(j)*L(j+1:n,j);
12 end
13 x(n) = b(n)/L(n,n);
syntax highlighted by Code2HTML, v. 0.8.8b
Note that the syntax is very C like.
To use this function, we would enter the following
commands at the Matlab prompt. For now, we are assuming that you
are running Matlab in a local directory /LOCAL
and that your Matlab code LTriSol.m is also in this directory.
So we fire up Matlab and enter these commands:
LTriSol.prompts
1 albert:Source) matlab
2 ----------------------------------------------------------
3 Your MATLAB license will expire in 57 days.
4 Please contact your system administrator or
5 The MathWorks to renew this license.
6 ----------------------------------------------------------
7
8 < M A T L A B >
9 Copyright 1984-1999 The MathWorks, Inc.
10 Version 5.3.1.29215a (R11.1)
11 Oct 6 1999
12
13
14 To get started, type one of these: helpwin, helpdesk, or demo.
15 For product information, type tour or visit www.mathworks.com.
16
17 >> A = [2 0 0; 1 5 0; 7 9 8]
18
19 A =
20
21 2 0 0
22 1 5 0
23 7 9 8
24
25 >> b = [6; 2; 5]
26
27 b =
28
29 6
30 2
31 5
32
33 >> x = LTriSol(A,b)
34
35 x =
36
37 3.0000
38 -0.2000
39 -1.7750
syntax highlighted by Code2HTML, v. 0.8.8b
which solves the system as we wanted.
1.1.2 An Upper Triangular Solver:
Here is a simple function to solve a similar system
where A is upper triangular.
UTriSol.m
1 function x = UTriSol(U,b)
2 %
3 % U is nxn nonsingular Upper Triangular matrix
4 % b is nx1 data vector
5 % x is solved by back substitution
6 %
7 n = length(b);
8 x = zeros(n,1);
9 for j = n:-1:2
10 x(j) = b(j)/U(j,j);
11 b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
12 end
13 x(1) = b(1)/U(1,1);
syntax highlighted by Code2HTML, v. 0.8.8b
As usual, to use this function, we would enter the following
commands at the Matlab prompt. We are still assuming that you
are running Matlab in a local directory /LOCAL
and that your Matlab code UTriSol.m is also in this directory.
So we fire up Matlab and enter these commands:
UTriSol.prompts
1 >> C = [7 9 8; 0 1 5; 0 0 2]
2
3 C =
4
5 7 9 8
6 0 1 5
7 0 0 2
8
9 >> b = [6; 2; 5]
10
11 b =
12
13 6
14 2
15 5
16
17 >> x = UTriSol(C,b)
18
19 x =
20
21 11.5000
22 -10.5000
23 2.5000
syntax highlighted by Code2HTML, v. 0.8.8b
which solves the system as we wanted.
1.1.3 The LU Decomposition of A Without Pivoting:
Here is a simple function to solve a system
using the LU decomposition of A.
GE.m
1 function [L,U] = GE(A)
2 %
3 % A is nxn matrix
4 % L is nxn lower triangular
5 % U is nxn upper triangular
6 %
7 % We compute the LU decomposition of A using
8 % Gaussian Elimination
9 %
10
11 [n,n] = size(A);
12 for k=1:n-1
13 % find multiplier
14 A(k+1:n,k) = A(k+1:n,k)/A(k,k);
15 % zero out column
16 A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
17 end
18 L = eye(n,n) + tril(A,-1);
19 U = triu(A);
syntax highlighted by Code2HTML, v. 0.8.8b
We fire up Matlab and enter these commands:
LU.prompts
1 >> A = [17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; ...
2 11 18 25 2 9]
3
4 A =
5
6 17 24 1 8 15
7 23 5 7 14 16
8 4 6 13 20 22
9 10 12 19 21 3
10 11 18 25 2 9
11
12 >> [L,U] = GE(A);
13 >> L
14
15 L =
16
17 1.0000 0 0 0 0
18 1.3529 1.0000 0 0 0
19 0.2353 -0.0128 1.0000 0 0
20 0.5882 0.0771 1.4003 1.0000 0
21 0.6471 -0.0899 1.9366 4.0578 1.0000
22
23 >> U
24
25 U =
26
27 17.0000 24.0000 1.0000 8.0000 15.0000
28 0 -27.4706 5.6471 3.1765 -4.2941
29 0 0 12.8373 18.1585 18.4154
30 0 0 0 -9.3786 -31.2802
31 0 0 0 0 90.1734
32
33 >> b = [1; 3; 5; 7; 9]
34
35 b =
36
37 1
38 3
39 5
40 7
41 9
42
43 >> y = LTriSol(L,b)
44
45 y =
46
47 1.0000
48 1.6471
49 4.7859
50 -0.4170
51 0.9249
52
53 >> x = UTriSol(U,y)
54
55 x =
56
57 0.0103
58 0.0103
59 0.3436
60 0.0103
61 0.0103
62
63 >> c = A*x
64
65 c =
66
67 1.0000
68 3.0000
69 5.0000
70 7.0000
71 9.0000
72
73
syntax highlighted by Code2HTML, v. 0.8.8b
which solves the system as we wanted.
1.1.4 The LU Decomposition of A With Pivoting:
Here is a simple function to solve a system
using the LU decomposition of A with what is
called pivoting. This means we find the largest
absolute value entry in the column we are trying to
zero out and perform row interchanges to bring that
one to the pivot position.
GePiv.m
1 function [L,U,piv] = GePiv(A);
2 %
3 % A is nxn matrix
4 % L is nxn lower triangular matrix
5 % U is nxn upper triangular matrix
6 % piv is a nx1 integer vector to hold variable order
7 % permutations
8 %
9 [n,n] = size(A);
10 piv = 1:n;
11 for k=1:n-1
12 [maxc,r] = max(abs(A(k:n,k)));
13 q = r+k-1;
14 piv([k q]) = piv([q k]);
15 A([k q],:) = A([q k],:);
16 if A(k,k) ~=0
17 A(k+1:n,k) = A(k+1:n,k)/A(k,k);
18 A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
19 end
20 end
21 L = eye(n,n) + tril(A,-1);
22 U = triu(A);
23
24
syntax highlighted by Code2HTML, v. 0.8.8b
We fire up Matlab and enter these commands:
GePiv.prompts
1 >> A = [17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; ...
2 11 18 25 2 9]
3
4 A =
5
6 17 24 1 8 15
7 23 5 7 14 16
8 4 6 13 20 22
9 10 12 19 21 3
10 11 18 25 2 9
11
12 >> b = [1; 3; 5; 7; 9]
13
14 b =
15
16 1
17 3
18 5
19 7
20 9
21
22 >> [L,U,piv] = GePiv(A);
23 >> L
24
25 L =
26
27 1.0000 0 0 0 0
28 0.7391 1.0000 0 0 0
29 0.4783 0.7687 1.0000 0 0
30 0.1739 0.2527 0.5164 1.0000 0
31 0.4348 0.4839 0.7231 0.9231 1.0000
32
33 >> U
34
35 U =
36
37 23.0000 5.0000 7.0000 14.0000 16.0000
38 0 20.3043 -4.1739 -2.3478 3.1739
39 0 0 24.8608 -2.8908 -1.0921
40 0 0 0 19.6512 18.9793
41 0 0 0 0 -22.2222
42
43 >> piv
44
45 piv =
46
47 2 1 5 3 4
48
49 >> y = LTriSol(L,b(piv));
50 >> y
51
52 y =
53
54 3.0000
55 -1.2174
56 8.5011
57 0.3962
58 -0.2279
59
60 >> x = UTriSol(U,y);
61 >> x
62
63 x =
64
65 0.0103
66 0.0103
67 0.3436
68 0.0103
69 0.0103
70
71 >> c = A*x
72
73 c =
74
75 1.0000
76 3.0000
77 5.0000
78 7.0000
79 9.0000
syntax highlighted by Code2HTML, v. 0.8.8b
which solves the system as we wanted.