• 沒有找到結果。

# Direct Methods for Solving Linear Systems

N/A
N/A
Protected

Academic year: 2022

Share "Direct Methods for Solving Linear Systems"

Copied!
89
0
0

(1)

logo

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

E-mail: min@math.ntnu.edu.tw

September 12, 2015

1 / 89

(2)

logo

### Outline

1 Linear systems of equations

2 Pivoting Strategies

3 Matrix factorization

4 Special types of matrices

(3)

logo

### Linear systems of equations

Three operations to simplify the linear system:

1 (λEi) → (Ei): Equation Ei can be multiplied by λ 6= 0 with the resulting equation used in place of Ei.

2 (Ei+ λEj) → (Ei): Equation Ej can be multiplied by λ 6= 0 and added to equation Eiwith the resulting equation used in place of Ei.

3 (Ei) ↔ (Ej): Equation Ei and Ej can be transposed in order.

Example 1

E1 : x1 + x2 + 3x4 = 4,

E2 : 2x1 + x2 − x3 + x4 = 1, E3 : 3x1 − x2 − x3 + 2x4 = −3, E4 : −x1 + 2x2 + 3x3 − x4 = 4.

3 / 89

(4)

logo

### Solution:

(E2− 2E1) → (E2), (E3− 3E1) → (E3)and (E4+ E1) → (E4):

E1: x1 + x2 + 3x4 = 4,

E2: − x2 − x3 − 5x4 = −7, E3: − 4x2 − x3 − 7x4 = −15, E4: 3x2 + 3x3 + 2x4 = 8.

(E3− 4E2) → (E3)and (E4+ 3E2) → (E4):

E1: x1 + x2 + 3x4 = 4,

E2: − x2 − x3 − 5x4 = −7,

E3: 3x3 + 13x4 = 13,

E4: − 13x4 = −13.

(5)

logo

Backward-substitution process:

1 E4 ⇒ x4= 1

2 Solve E3for x3: x3=1

3(13 − 13x4) =1

3(13 − 13) = 0.

3 E2gives

x2= −(−7 + 5x4+ x3) = −(−7 + 5 + 0) = 2.

4 E1gives

x1= 4 − 3x4− x2= 4 − 3 − 2 = −1.

5 / 89

(6)

logo

Solve linear systems of equations









a11x1+ a12x2+ · · · + a1nxn = b1

a21x1+ a22x2+ · · · + a2nxn = b2

...

an1x1+ an2x2+ · · · + annxn = bn

Rewrite in the matrix form

Ax = b, (1)

where

A =

a11 a12 · · · a1n a21 a22 · · · a2n

... ... . .. ... an1 an2 · · · ann

 , b =

 b1 b2

... bn

, x =

 x1 x2

... xn

 and [A, b] is called the augmented matrix.

(7)

logo

### Gaussian elimination with backward substitution

The augmented matrix in previous example is

1 1 0 3 4

2 1 −1 1 1

3 −1 −1 2 −3

−1 2 3 −1 4

.

(E2− 2E1) → (E2), (E3− 3E1) → (E3)and (E4+ E1) → (E4):

1 1 0 3 4

0 −1 −1 −5 −7

0 −4 −1 −7 −15

0 3 3 2 8

.

(E3− 4E2) → (E3)and (E4+ 3E2) → (E4):

1 1 0 3 4

0 −1 −1 −5 −7

0 0 3 13 13

0 0 0 −13 −13

.

7 / 89

(8)

logo

### The general Gaussian elimination procedure

Provided a116= 0, for each i = 2, 3, . . . , n,



Ei ai1 a11

E1



→ (Ei).

Transform all the entries in the first col. below the diagonal are zero. Denote the new entry in the ith row and jth col. by aij. For i = 2, 3 . . . , n − 1, providedaii 6= 0,



Ejaji

aiiEi



→ (Ej), ∀ j = i + 1, i + 2, . . . , n.

Transform all the entries in the ith column below the diagonal are zero.

Result anupper triangularmatrix:

a11 a12 · · · a1n b1

0 a22 · · · a2n b2

... . .. ... ... ... 0 · · · 0 ann bn

.

(9)

logo

The process of Gaussian elimination result in a sequence of matrices as follows:

A = A(1) → A(2)→ · · · → A(n)=upper triangular matrix The matrix A(k)has the following form:

A(k)=

a(1)11 · · · a(1)1,k−1 a(1)1k · · · a(1)1j · · · a(1)1n

... . .. ... ... ... ...

0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k · · · a(k−1)k−1,j · · · a(k−1)k−1,n 0 · · · 0 a(k)kk · · · a(k)kj · · · a(k)kn

... ... ... ... ...

0 · · · 0 a(k)ik · · · a(k)ij · · · a(k)in

... ... ... ... ...

0 · · · 0 a(k)nk · · · a(k)nj · · · a(k)nn

9 / 89

(10)

logo

The entries of A(k)are produced by the formula

a(k)ij =





a(k−1)ij , for i = 1, . . . , k − 1, j = 1, . . . , n;

0, for i = k, . . . , n, j = 1, . . . , k − 1;

a(k−1)ija

(k−1) i,k−1

a(k−1)k−1,k−1 × a(k−1)k−1,j, for i = k, . . . , n, j = k, . . . , n.

The procedure will fail if one of the elements a(1)11, a(2)22, . . . , a(n)nn is zero.

a(i)ii is called the pivot element.

(11)

logo

### Backward substitution

The new linear system is triangular:

a11x1 + a12x2 + · · · + a1nxn = b1, a22x2 + · · · + a2nxn = b2,

...

annxn = bn Solving the nth equation for xngives

xn= bn ann

.

Solving the (n − 1)th equation for xn−1and using the value for xnyields

xn−1= bn−1− an−1,nxn

an−1,n−1

. In general,

xi= biPn

j=i+1aijxj aii

, ∀ i = n − 1, n − 2, . . . , 1.

11 / 89

(12)

logo

Algorithm 1 (Backward Substitution)

Suppose thatU ∈ Rn×n isnonsingular upper triangularand b ∈ Rn. This algorithm computes the solution ofU x = b.

For i = n, . . . , 1 tmp = 0

For j = i + 1, . . . , n

tmp = tmp + U (i, j) ∗ x(j) End for

x(i) = (b(i) − tmp)/U (i, i) End for

(13)

logo

Example 2

Solve system of linear equations.

6 −2 2 4

12 −8 6 10

3 −13 9 3

−6 4 1 −18

 x1 x2

x3

x4

=

 12 34 27

−38

Solution:

1st step Use 6 as pivot element, the first row as pivot row, and multipliers 2,12, −1are produced to reduce the system to

6 −2 2 4

0 −4 2 2

0 −12 8 1

0 2 3 −14

 x1 x2

x3

x4

=

 12 10 21

−26

13 / 89

(14)

logo

2nd step Use −4 as pivot element, the second row as pivot row, and multipliers 3, −12 are computed to reduce the system to

6 −2 2 4

0 −4 2 2

0 0 2 −5

0 0 4 −13

 x1

x2 x3

x4

=

 12 10

−9

−21

3rdstep Use 2 as pivot element, the third row as pivot row, and multipliers 2 is found to reduce the system to

6 −2 2 4

0 −4 2 2

0 0 2 −5

0 0 0 −3

 x1 x2 x3

x4

=

 12 10

−9

−3

(15)

logo

4thstep The backward substitution is applied:

x4 = −3

−3 = 1, x3 = −9 + 5x4

2 = −9 + 5 2 = −2, x2 = 10 − 2x4− 2x3

−4 = 10 − 2 + 4

−4 = −3, x1 = 12 − 4x4− 2x3+ 2x2

6 = 12 − 4 + 4 − 6

6 = 1.

This example is done since a(k)kk 6= 0 for all k = 1, 2, 3, 4.

How to do if a(k)kk = 0for some k?

15 / 89

(16)

logo

Example 3

Solve system of linear equations.

1 −1 2 −1 2 −2 3 −3

1 1 1 0

1 −1 4 3

 x1 x2

x3

x4

=

−8

−20

−2 4

Solution:

1st step Use 1 as pivot element, the first row as pivot row, and multipliers 2, 1, 1 are produced to reduce the system to

1 −1 2 −1

0 0 −1 −1

0 2 −1 1

0 0 2 4

 x1 x2

x3

x4

=

−8

−4 6 12

(17)

logo

2nd step Since a(2)22 = 0and a(2)32 6= 0, the operation

(E2) ↔ (E3)is performed to obtain a new system

1 −1 2 −1

0 2 −1 1

0 0 −1 −1

0 0 2 4

 x1

x2 x3

x4

=

−8 6

−4 12

3rdstep Use −1 as pivot element, the third row as pivot row, and multipliers −2 is found to reduce the system to

1 −1 2 −1

0 2 −1 1

0 0 −1 −1

0 0 0 2

 x1 x2

x3 x4

=

−8 6

−4 4

17 / 89

(18)

logo

4thstep The backward substitution is applied:

x4 = 4 2 = 2, x3 = −4 + x4

−1 = 2, x2 = 6 − x4+ x3

2 = 3,

x1 = −8 + x4− 2x3+ x2

1 = −7.

This example illustrates what is done if a(k)kk = 0for some k.

If a(k)pk 6= 0 for some p with k + 1 ≤ p ≤ n, then the operation (Ek) ↔ (Ep)is performed to obtain new matrix.

If a(k)pk = 0for each p, then the linear system does not have a unique solution and the procedure stops.

(19)

logo

Algorithm 2 (Gaussian elimination)

GivenA ∈ Rn×n andb ∈ Rn, this algorithm implements the Gaussian elimination procedure to reduceAtoupper triangular and modify the entries ofbaccordingly.

For k = 1, . . . , n − 1

Let p be the smallest integer with k ≤ p ≤ n and apk 6= 0.

If @ p, then stop.

If p 6= k, then perform (Ep) ↔ (Ek).

For i = k + 1, . . . , n t = A(i, k)/A(k, k) A(i, k) = 0

b(i) = b(i) − t × b(k) For j = k + 1, . . . , n

A(i, j) = A(i, j) − t × A(k, j) End for

End for

End for 19 / 89

(20)

logo

### Number of floating-point arithmetic operations

Eliminate kth column For i = k + 1, . . . , n

t = A(i, k)/A(k, k); b(i) = b(i) − t × b(k).

For j = k + 1, . . . , n

A(i, j) = A(i, j) − t × A(k, j) End for

End for

Multiplications/divisions

(n − k) + (n − k) + (n − k)(n − k) = (n − k)(n − k + 2) Additions/subtractions

(n − k) + (n − k)(n − k) = (n − k)(n − k + 1)

(21)

logo

Total number of operations for multiplications/divisions

n−1

X

k=1

(n − k)(n − k + 2) =

n−1

X

k=1

(n2− 2nk + k2+ 2n − 2k)

= (n2+ 2n)

n−1

X

k=1

1 − 2(n + 1)

n−1

X

k=1

k +

n−1

X

k=1

k2

= (n2 + 2n)(n − 1) − 2(n + 1)(n − 1)n

2 +(n − 1)n(2n − 1) 6

= 2n3+ 3n2− 5n

6 .

Total number of operations for additions/subtractions

n−1

X

k=1

(n − k)(n − k + 1) =

n−1

X

k=1

(n2− 2nk + k2+ n − k)

= (n2+ n)

n−1

X

k=1

1 − (2n + 1)

n−1

X

k=1

k +

n−1

X

k=1

k2= n3− n 3 .

21 / 89

(22)

logo

Backward substitution

x(n) = b(n)/U (n, n).

For i = n − 1, . . . , 1

tmp = U (i, i + 1) × x(i + 1) For j = i + 2, . . . , n

tmp = tmp + U (i, j) × x(j) End for

x(i) = (b(i) − tmp)/U (i, i) End for

Multiplications/divisions 1 +

n−1

X

i=1

[(n − i) + 1] =n2+ n 2 Additions/subtractions

n−1

X

i=1

[(n − i − 1) + 1] =n2− n 2

(23)

logo

The total number of arithmetic operations in Gaussian elimination with backward substitution is:

Multiplications/divisions 2n3+ 3n2− 5n

6 + n2+ n 2 = n3

3 + n2−n 3 ≈ n3

3 Additions/subtractions

n3− n

3 +n2− n 2 = n3

3 + n2 2 −5n

6 ≈ n3 3

23 / 89

(24)

logo

Exercise

Page 368: 5, 10, 12, 15

(25)

logo

### Pivoting Strategies

If a(k)kk is small in magnitude compared to a(k)jk , then

|mjk| =

a(k)jk a(k)kk

> 1.

Round-off error introduced in the computation of a(k+1)j` = a(k)j` − mjka(k)k`, for ` = k + 1, . . . , n.

Error can be increased when performing the backward substitution for

xk = bk−Pn

j=k+1a(k)kjxj a(k)kk

with a small value of a(k)kk.

25 / 89

(26)

logo

Example 4 The linear system

E1 : 0.003000x1 + 59.14x2 = 59.17, E2 : 5.291x1 − 6.130x2 = 46.78,

has the exact solution x1 = 10.00and x2 = 1.000. Suppose Gaussian elimination is performed on this system using four-digit arithmetic with rounding.

a11= 0.0030is small and m21= 5.291

0.0030 = 1763.6¯6 ≈ 1764.

Perform (E2− m21E1) → (E2):

0.0030x1 + 59.14x2 = 59.17

− 104309.37¯6x2 = −104309.37¯6.

(27)

logo

Rounding with four-digit arithmetic:

Coefficient of x2:

−6.130 − 1764 × 59.14 = −6.130 − 104322.96

≈ −6.130 − 104300 = −104306.13

≈ −104300.

Right hand side:

46.78 − 1764 × 59.17 = 46.78 − 104375.88

≈ 46.78 − 104400 = −104353.22

≈ −104400.

New linear system:

0.0030x1 + 59.14x2 = 59.17

− 104300x2 ≈ −104400.

27 / 89

(28)

logo

Approximated solution:

x2 = 104400

104300 ≈ 1.001, x1 = 59.17 − 59.14 × 1.001

0.0030 = 59.17 − 59.19914 0.0030

≈ 59.17 − 59.20

0.0030 = −10.00.

This ruins the approximation to the actual value x1 = 10.00.

(29)

logo

### Partial pivoting

To avoid the pivot element small relative to other entries, pivoting is performed by selecting an element a(k)pq with a larger magnitude as the pivot.

Specifically, select pivoting a(k)pk with

|a(k)pk| = max

k≤i≤n|a(k)ik | and perform (Ek) ↔ (Ep).

This row interchange strategy is called partial pivoting.

29 / 89

(30)

logo

Example 5

Reconsider the linear system

E1 : 0.003000x1 + 59.14x2 = 59.17, E2 : 5.291x1 − 6.130x2 = 46.78.

Find pivoting with

max{|a11|, |a21|} = 5.291 = |a21|.

Perform (E2) ↔ (E1):

E1 : 5.291x1 − 6.130x2 = 46.78, E2 : 0.003000x1 + 59.14x2 = 59.17.

The multiplier for new system is m21= a21

a11

= 0.0005670.

(31)

logo

The operation (E2− m21E1) → (E2)reduces the system to 5.291x1 − 6.130x2 = 46.78,

59.14x2 ≈ 59.14.

The four-digit answers resulting from the backward substitution are the correct values x1 = 10.00and x2= 1.000.

31 / 89

(32)

logo

Example 6 The linear system

E1: 30.00x1 + 591400x2 = 591700, E2: 5.291x1 − 6.130x2 = 46.78,

is the same as that in previous example except that all the entries in the first equation have been multiplied by 104. The pivoting is a11= 30.00and the multiplier

m21= 5.291

30.00 = 0.1764 leads to the system

30.00x1 + 591400x2 = 591700

− 104300x2 ≈ −104400,

which has inaccurate solution x2 ≈ 1.001 and x1 ≈ −10.00.

(33)

logo

### Scaled partial pivoting

Define a scale factor si as si = max

1≤j≤n|aij|, for i = 1, . . . , n.

If si= 0for some i, then the system has no unique solution.

In the ith column, choose the least integer p ≥ i with

|api|

sp = max

i≤k≤n

|aki| sk and perform (Ei) ↔ (Ep)if p 6= i.

The scale factors s1, . . . , snare computed only once and must also be interchanged when row interchanges are performed.

33 / 89

(34)

logo

Example 7

Apply scaled partial pivoting to the linear system E1: 30.00x1 + 591400x2 = 591700, E2: 5.291x1 6.130x2 = 46.78.

The scale factors s1and s2are

s1= max{|30.00|, |591400|} = 591400 and

s2= max{|5.291|, | − 6.130|} = 6.130.

Consequently,

|a11| s1

= 30.00

591400 = 0.5073 × 10−4,

|a21|

s2 = 5.291

6.130 = 0.8631, and the interchange (E1) ↔ (E2)is made.

(35)

logo

Applying Gaussian elimination to the new system 5.291x1 − 6.130x2 = 46.78, 30.00x1 + 591400x2 = 591700

produces the correct results: x1 = 10.00and x2 = 1.000.

35 / 89

(36)

logo

Exercise

Page 379: 2, 4, 6, 31

(37)

logo

### Matrix factorization

This equation has a unique solutionx = A−1bwhen the coefficient matrix A isnonsingular.

Use Gaussian elimination to factor the coefficient matrix into a product of matrices. The factorization is called LU-factorizationand has the formA = LU, whereLisunit lower triangularandU isupper triangular.

The solution to the original problemAx = LU x = bis then found by a two-step triangular solve process:

Ly = b, U x = y.

LU factorization requiresO(n3)arithmetic operations.

Forward substitution for solving a lower-triangular system Ly = brequiresO(n2). Backward substitution for solving an upper-triangular system U x = y requiresO(n2) arithmetic operations.

37 / 89

(38)

logo

A =

1 1 0 3

2 1 −1 1

3 −1 −1 2

−1 2 3 −1

⇒ A1 := L1A ≡

1 0 0 0

−2 1 0 0

−3 0 1 0 1 0 0 1

 A =

1 1 0 3

0 −1 −1 −5 0 −4 −1 −7

0 3 3 2

⇒ A2 := L2A1

1 0 0 0

0 1 0 0

0 −4 1 0

0 3 0 1

 A1 =

1 1 0 3

0 −1 −1 −5

0 0 3 13

0 0 0 −13

= L2L1A

(39)

logo

We have

A = L−11 L−12 A2= LR.

where L and R are lower and upper triangular, respectively.

Question

How to compute L−11 and L−12 ?

L1 =

1 0 0 0

−2 1 0 0

−3 0 1 0 1 0 0 1

= I −

 0 2 3

−1

 1 0 0 0 

L2 =

1 0 0 0

0 1 0 0

0 −4 1 0

0 3 0 1

= I −

 0 0 4

−3

 0 1 0 0 

39 / 89

(40)

logo

Since

 I −

 0 2 3

−1

 1 0 0 0 

 I +

 0 2 3

−1

 1 0 0 0 

= I,

we have

L−11 =

1 0 0 0

−2 1 0 0

−3 0 1 0 1 0 0 1

−1

=

1 0 0 0 2 1 0 0 3 0 1 0

−1 0 0 1

(41)

logo

Since

 I −

 0 0 4

−3

 0 1 0 0 

 I +

 0 0 4

−3

 0 1 0 0 

= I,

we have

L−12 =

1 0 0 0

0 1 0 0

0 −4 1 0

0 3 0 1

−1

=

1 0 0 0

0 1 0 0

0 4 1 0

0 −3 0 1

41 / 89

(42)

logo

By the fact

L−11 L−12 =

1 0 0 0 2 1 0 0 3 0 1 0

−1 0 0 1

1 0 0 0

0 1 0 0

0 4 1 0

0 −3 0 1

=

1 0 0 0

2 1 0 0

3 4 1 0

−1 −3 0 1

 ,

it holds that

1 1 0 3

2 1 −1 1

3 −1 −1 2

−1 2 3 −1

=

1 0 0 0

2 1 0 0

3 4 1 0

−1 −3 0 1

1 1 0 3

0 −1 −1 −5

0 0 3 13

0 0 0 −13

 .

(43)

logo

For a given vectorv ∈ Rnwithvk6= 0for some 1 ≤ k ≤ n, let

`ik = vi

vk, i = k + 1, . . . , n,

`k=

0 · · · 0 `k+1,k · · · `n,k T

, and

Mk= I − `keTk =

1 · · · 0 0 · · · 0 ... . .. ... ... ... 0 · · · 1 0 · · · 0 0 · · · −`k+1,k 1 · · · 0 ... ... ... . .. ...

0 · · · −`n,k 0 · · · 1

 .

43 / 89

(44)

logo

Then one can verify that Mkv =

v1 · · · vk 0 · · · 0 T

.

Mk is called aGaussian transformation, the vector`kaGauss vector. Furthermore, one can verify that

Mk−1 = (I − `keTk)−1 = I + `keTk =

1 · · · 0 0 · · · 0 ... . .. ... ... ... 0 · · · 1 0 · · · 0 0 · · · `k+1,k 1 · · · 0 ... ... ... . .. ...

0 · · · `n,k 0 · · · 1

 .

(45)

logo

Given a nonsingular matrix A ∈ Rn×n, denote A(1)≡ [a(1)ij ] = A.

Ifa(1)11 6= 0, then

M1 = I − `1eT1, where

`1=

0 `21 · · · `n1

T

, `i1= a(1)i1

a(1)11, i = 2, . . . , n, can be formed such that

A(2) = M1A(1) =

a(1)11 a(1)12 · · · a(1)1n 0 a(2)22 · · · a(2)2n ... ... . .. ... 0 a(2)n2 · · · a(2)nn

 ,

where

a(2)ij = a(1)ij − `i1× a(1)1j , for i = 2, . . . , n and j = 2, . . . , n.

45 / 89

(46)

logo

In general, at the k-th step, we are confronted with a matrix A(k) = Mk−1· · · M2M1A(1)

=

a(1)11 a(1)12 · · · a(1)1,k−1 a(1)1k · · · a(1)1n 0 a(2)22 · · · a(2)2,k−1 a(2)2k · · · a(2)2n ... ... . .. ... ... ... 0 0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k · · · a(k−1)k−1,n 0 0 · · · 0 a(k)kk · · · a(k)kn

... ... ... ... . .. ... 0 0 · · · 0 a(k)kn · · · a(k)nn

 .

If the pivota(k)kk 6= 0, then the multipliers

`ik = a(k)ik a(k)kk

, i = k + 1, . . . , n,

(47)

logo

can be computed and the Gaussian transformation Mk = I − `keTk, where `k=

0 · · · 0 `k+1,k · · · `nk

T

, can be applied to the left of A(k)to obtain

A(k+1) = MkA(k)

=

a(1)11 a(1)12 · · · a(1)1,k−1 a(1)1k a(1)1,k+1 · · · a(1)1n 0 a(2)22 · · · a(2)2,k−1 a(2)2k a(2)2,k+1 · · · a(2)2n ... ... . .. ... ... ... ... 0 0 · · · a(k−1)k−1,k−1 a(k−1)k−1,k a(k−1)k−1,k+1 · · · a(k−1)k−1,n 0 0 · · · 0 a(k)kk a(k)k,k+1 · · · a(k)kn

... ... ... 0 a(k+1)k+1,k+1 · · · a(k+1)k+1,n

... ... ... ... ... ...

0 0 · · · 0 0 a(k+1)n,k+1 · · · a(k+1)nn

 ,

47 / 89

(48)

logo

in which

a(k+1)ij = a(k)ij − `ika(k)kj , (2) for i = k + 1, . . . , n, j = k + 1, . . . , n. Upon the completion,

U ≡ A(n)= Mn−1· · · M2M1A isupper triangular. Hence

A = M1−1M2−1· · · Mn−1−1 U ≡ LU,

(49)

logo

where

L ≡ M1−1· · · Mn−1−1 = (I − `1eT1)−1(I − `2eT2)−1· · · (I − `n−1eTn−1)−1

= (I + `1eT1)(I + `2eT2) · · · (I + `n−1eTn−1)

= I + `1eT1 + `2eT2 + · · · + `n−1eTn−1

=

1 0 0 · · · 0

`21 1 0 · · · 0

`31 `32 1 · · · 0 ... ... ... . .. ...

`n1 `n2 `n3 · · · 1

isunit lower triangular. This matrix factorization is called the LU-factorizationofA.

49 / 89

(50)

logo

Algorithm 3 (LU Factorization)

Given anonsingular squarematrixA ∈ Rn×n, this algorithm computes aunit lower triangularmatrixLand anupper triangularmatrixU such thatA = LU. The matrix A is overwrittenby L and U .

For k = 1, . . . , n − 1 For i = k + 1, . . . , n

A(i, k) = A(i, k)/A(k, k) For j = k + 1, . . . , n

A(i, j) = A(i, j) − A(i, k) × A(k, j) End for

End for End for

(51)

logo

### Forward Substitution

When a linear systemLx = bislower triangularof the form

`11 0 · · · 0

`21 `22 · · · 0 ... ... . .. ...

`n1 `n2 · · · `nn

 x1 x2

... xn

=

 b1 b2

... bn

 ,

where all diagonals`ii6= 0, xican be obtained by the following procedure

x1 = b1/`11,

x2 = (b2− `21x1)/`22,

x3 = (b3− `31x1− `32x2)/`33, ...

xn = (bn− `n1x1− `n2x2− · · · − `n,n−1xn−1)/`nn.

51 / 89

(52)

logo

The general formulation for computing xi is xi =

bi

i−1

X

j=1

`ijxj

`ii, i = 1, 2, . . . , n.

Algorithm 4 (Forward Substitution)

Suppose thatL ∈ Rn×n isnonsingular lower triangularand b ∈ Rn. This algorithm computes the solution of Lx = b.

For i = 1, . . . , n tmp = 0

For j = 1, . . . , i − 1

tmp = tmp + L(i, j) ∗ x(j) End for

x(i) = (b(i) − tmp)/L(i, i) End for

(53)

logo

Example 8

E1 : x1 + x2 + 3x4 = 4,

E2 : 2x1 + x2 − x3 + x4 = 1, E3 : 3x1 − x2 − x3 + 2x4 = −3, E4 : −x1 + 2x2 + 3x3 − x4 = 4.

Solution:

The sequence {(E2− 2E1) → (E2), (E3− 3E1) → (E3), (E4− (−1)E1) → (E4), (E3− 4E2) → (E3),

(E4− (−3)E2) → (E4)}converts the system to the triangular system

x1 + x2 + 3x4 = 4,

− x2 − x3 − 5x4 = −7, 3x3 + 13x4 = 13,

− 13x4 = −13.

53 / 89

(54)

logo

LU factorization of A:

A =

1 1 0 3

2 1 −1 1

3 −1 −1 2

−1 2 3 −1

=

1 0 0 0

2 1 0 0

3 4 1 0

−1 −3 0 1

1 1 0 3

0 −1 −1 −5

0 0 3 13

0 0 0 −13

= LU.

(55)

logo

Solve Ly = b:

1 0 0 0

2 1 0 0

3 4 1 0

−1 −3 0 1

 y1

y2 y3 y4

=

 8 7 14

−7

which implies that y1 = 8,

y2 = 7 − 2y1 = −9, y3 = 14 − 3y1− 4y2= 26, y4 = −7 + y1+ 3y2 = −26.

55 / 89

(56)

logo

Solve U x = y:

1 1 0 3

0 −1 −1 −5

0 0 3 13

0 0 0 −13

 x1 x2

x3

x4

=

 8

−9 26

−26

which implies that x4 = 2,

x3 = (26 − 13x4)/3 = 0,

x2 = (−9 + 5x4+ x3)/(−1) = −1, x1 = 8 − 3x4− x2 = 3.

(57)

logo

### Partial pivoting

At the k-th step, select pivoting a(k)pk with

|a(k)pk| = max

k≤i≤n|a(k)ik |

and perform (Ek) ↔ (Ep). That is, choose a permutation matrix

Pk=

Ik−1 0 0 0 0

0 0 0 1 0

0 0 Ip−k−1 0 0

0 1 0 0 0

0 0 0 0 In−p

so that

(PkA(k))kk

= max

k≤i≤n

(A(k))ik and

A(k+1)= M(k)PkA(k).

57 / 89

(58)

logo

LetP1, . . . , Pk−1 be thepermutationschosen andM1, . . . Mk−1 denote theGaussian transformationsperformed in the first k − 1steps. At the k-th step, a permutation matrixPk is chosen so that

|(PkMk−1· · · M1P1A)kk| = max

k≤i≤n|(Mk−1· · · M1P1A)ik| . As a consequence,|`ij| ≤ 1for i = 1, . . . , n, j = 1, . . . , i. Upon completion, we obtain an upper triangular matrix

U ≡ Mn−1Pn−1· · · M1P1A. (3) Since anyPkissymmetricandPkTPk= Pk2 = I, we have

Mn−1Pn−1· · · M2P2M1P2· · · Pn−1Pn−1· · · P2P1A = U, therefore,

Pn−1· · · P1A = (Mn−1Pn−1· · · M2P2M1P2· · · Pn−1)−1U.

(59)

logo

In summary,Gaussian elimination with partial pivotingleads to theLU factorization

P A = LU, (4)

where

P = Pn−1· · · P1 is a permutation matrix, and

L ≡ (Mn−1Pn−1· · · M2P2M1P2· · · Pn−1)−1

= Pn−1· · · P2M1−1P2M2−1· · · Pn−1Mn−1−1 . Since

Pj =

Ij−1 0 0 0 0

0 0 0 1 0

0 0 Ip−j−1 0 0

0 1 0 0 0

0 0 0 0 In−p

, `j =

 0

... 0

`j+1,j ...

`nj

 ,

59 / 89

(60)

logo

it implies that for i < j,

eTi Pj = eTi , eTi `j = 0, Pj`i =

0 · · · 0 `˜i+1,i · · · `˜n,i T

≡ ˜`i,

P2M1−1P2 = P2(I + `1eT1)P2 = I + ˜`1eT1

P2M1−1P2M2−1 = (I + ˜`1eT1)(I + `2eT2) = I + ˜`1eT1 + `2eT2,

P3 P2M1−1P2M2−1 P3= I + ˆ`1eT1 + ˜`2eT2

⇒ · · ·

Therefore,Lisunit lower triangular.

(61)

logo

Algorithm 5 (LU -factorization with Partial Pivoting)

Given a nonsingular A ∈ Rn×n, this algorithm finds apermutationP, and computes aunit lower triangularLand anupper triangularU such thatP A = LU. A isoverwrittenby L and U , and P is not formed. Aninteger arraypis instead used for storing the row/column indices.

p(1 : n) = 1 : n For k = 1, . . . , n − 1

m = k

For i = k + 1, . . . , n

If|A(p(m), k)| < |A(p(i), k)|, thenm = i End For

` = p(k); p(k) = p(m); p(m) = ` For i = k + 1, . . . , n

A(p(i), k) = A(p(i), k)/A(p(k), k) For j = k + 1, . . . , n

A(p(i), j) = A(p(i), j) − A(p(i), k)A(p(k), j) End For

End For

End For 61 / 89

An additional senior teacher post, to be offset by a post in the rank of Certificated Master/Mistress or Assistant Primary School Master/ Mistress as appropriate, is provided

An additional senior teacher post, to be offset by a post in the rank of CM or Assistant Primary School Master/Mistress (APSM) as appropriate, is provided to each primary

An additional senior teacher post, to be offset by a post in the rank of CM or Assistant Primary School Master/Mistress (APSM) as appropriate, is provided to each primary

An additional senior teacher post, to be offset by a post in the rank of Certificated Master/Mistress or Assistant Primary School Master/Mistress as appropriate, is provided to

An additional senior teacher post, to be offset by a post in the rank of APSM, is provided to each primary special school/special school with primary section that operates six or

An additional senior teacher post, to be offset by a post in the rank of APSM, is provided to each primary special school/special school with primary section that operates six or

An additional senior teacher post, to be offset by a post in the rank of CM or APSM as appropriate, is provided to each primary special school/special school with

An additional senior teacher post, to be offset by a post in the rank of CM or Assistant Primary School Master/Mistress (APSM) as appropriate, is provided to each primary