Linear Equation. A b x. Computational Physics by Agus Naba, Ph.D.- UB

Linear Equation

Ax = b ; x = ? A : m x n matrix b : m - vector x : n - vector Computational Physics by Agus Naba, Ph.D.- UB

1

Example R1

R2 i1 V1 i2

R3

i3

V2

R5

R4 Computational Physics by Agus Naba, Ph.D.- UB

2

i R + (i − i )R + (i − i )R + V 1

1

1

2

2

1

3

3

(i − i )R + (i − i )R − V (i − i )R + i R + (i − i )R 2

3

1

1

2

3

2

3

3

4

5

3

Computational Physics by Agus Naba, Ph.D.- UB

2

2

5

1

=0 =0 =0

3

R + R + R  −R   − R 1

2

3

2

3

−R

−R

2

R +R −R 2

5

5

  i  − V       −R i = V     R + R + R  i   0  3

3

1

5

2

4


5

1

2

3

4

Existence & Uniqueness Existence and Uniqueness of a solution Ax=b depend on whether the matrix A is singular or nonsingular. Nonsingular Matrix satisfies the following: •A has an inverse, i.e., A-1 such that AA-1=I •Det(A) ≠0 •Rank(A)=n •For any vector z ≠0, Az ≠0 Computational Physics by Agus Naba, Ph.D.- UB

5

Jika matriks A adalah nonsingular, maka A punya inverse A-1, dan Ax=b selalu mempunyai solusi unik, untuk setiap b Jika matriks A adalah singular, maka solusi Ax=b bergantung pada b, bisa ada atau bisa tidak ada solusi.


6

Example 1

2 x + 3x = b 1

2

1

5x + 4 x = b 1

2


2 7

 2 3  x   b  Ax =  = = b      5 4  x  b  1

1

2

2

Solusi x adalah unik karena A adalah nonsingular, berapapun b.


8

Example 2

2 x + 3x = b 1

2

1

4x + 6x = b 1

2


2 9

 2 3  x   b  Ax =  = = b      4 6  x  b  1

1

2

2

x

Karena A adalah singular, solusi mungkin ada mungkin tidak !, bergantung pada b. Kalaupun ada solusi, pasti solusinya tidak unik


10

Jika b=[4 7]T, tidak ada solusi untuk x. Jika b=[4 8]T, maka solusinya:

γ   x=  (4 − 2 γ ) / 3 Computational Physics by Agus Naba, Ph.D.- UB

11

Problem Transformations Ax = b ; x = A b −1

MAz = Mb - > Transformasi z = (MA ) Mb = A M Mb = A b z=x −1

−1

−1

−1

Transformasi tidak mengubah solusi, malah bahkan bisa mempermudah menemukan solusi Computational Physics by Agus Naba, Ph.D.- UB

12

Example: Permutation 0 0 1   v   v  1 0 0 v  =  v       0 1 0 v  v  1

3

2

1

3

2

P Matrik permutasi P selalu nonsingular, dan berlaku P-1 = PT.

APz = b z = (AP ) b = P A b = P x −1

−1

−1

x = Pz Computational Physics by Agus Naba, Ph.D.- UB

−1

13

Example: Diagonal Scaling Matriks diagonal D= {dij} : semua elemen dij= 0 untuk i≠j.

ADz = b z = (AD) b = D A b = D x −1

−1

−1

−1

x = Dz Computational Physics by Agus Naba, Ph.D.- UB

14

Triangular Linear Systems 2  x   3  1 2 0 − 4 − 6   x  =  − 6       0 0 − 1  x   1  1

2

3

Triangular Matrix x3=-1; x2=3; x1=-1 Jika A adalah matrix triangular, solusi lebih mudah ditemukan ! →Lakukan transformasi matrik A menjadi matriks triangular ! Computational Physics by Agus Naba, Ph.D.- UB

15

Triangular Matrix Types • Lower Triangular Matrix L={lij}: semua elemen diatas elemen diagonal bernilai nol, yaitu lij = 0 untuk i < j • Upper Triangular Matrix U={lij}: semua elemen dibawah elemen diagonal bernilai nol, yaitu lij = 0 untuk i > j NB: Matrix L dan U dapat dipermutasikan menjadi U dan L dengan matrix permutasi yang sesuai Computational Physics by Agus Naba, Ph.D.- UB

16

Forward Substitution • Dilakukan dalam memecahkan problem Lx = b dengan persamaan berikut:

x = b / l 1

1

11

(

i −1

x = b − ∑ l x i

i

j =1

ij

j

)/ l

ii

i = 2 ,L , n .

-------------------------Pseudocode------------------------------for j = 1 to n {loop over columns} if ljj =0 then stop {stop if matrix is singular} xj=bj / ljj {compute solution component} for i=j+1 to n bi=bi – lij xj {update right-hand side} end end Computational Physics by Agus Naba, Ph.D.- UB

17

Backward Substitution • Dilakukan dalam memecahkan problem Ux = b dengan persamaan berikut:

x =b /l n

n

nn

(

n

)

x = b − ∑u x / u i

i

j =i +1

ij

j

ii

i = n − 1,L ,1.

-------------------------Pseudocode------------------------------for j = n to 1 {loop over columns} if uij =0 then stop {stop if matrix is singular} xj=bj / ujj {compute solution component} for i=1 to j-1 bi=bi – uij xj {update right-hand side} end end Computational Physics by Agus Naba, Ph.D.- UB

18

Elementary Elimination Matrix (Gaussian Transformation) Dipakai untuk mentransformasi sembarang matriks menjadi matriks triangular

1 M  0 M a= 0 M  0 m =a /a k

i

i

0 L 0  a   a  M O M  M   M      L 1 0 L 0  a  a    =   L −m 1 L 0 a   0  O M M O M  M   M      L − m 0 L 1  a   0  , i = k + 1,L ,n disebut multiplier L O

k

0 M

1

1

k

k

k +1

k +1

n

n


19

Properties of M • Mk : lower triangular matrix, nonsingular • Mk=I-mekT, dimana m=[0,…,0,mk+1,..,mn]T (multiplier vector) dan ek adalah kolom ke k matriks identitas • Mk-1 = I+mekT adalah sama dg Mk kecuali tanda elemenelemen di bawah diagonal adalah dibalik • Jika Mj, j>k, adalah matrik elementer yang lain sbgmn Mk,dengan multiplier vector t, maka MkMj=I-mekT -tekT +mekT tekT =I-mekT -mekT Note : ekTt= 0; Computational Physics by Agus Naba, Ph.D.- UB

20

Example a = [2 4 − 2]

T

Cari M1 dan M2 !


21

1 0  Ma= −2 1   1 0 1 0  Ma= 0 1  0 0.5 1

2

0  2   2      0 4 = 0     1 − 2 0 0  2  2      0 4 = 4     1 − 2 0


22

L =M 1

−1 1

; L =M 2

−1 2

0 0 0 0 1 1     M M = − 2 1 0 ;L L = 2 1 0      1 0.5 1 − 1 − 0.5 1 1

2

1

2


23

Gaussian Elimination Jika matrix Gaussian Elimination sudah ditemukan, maka Ax=b bisa dengan mudah ditransformasi menjadi bentuk upper triangular M1Ax=M1b → Kolom PERTAMA matrix A bernilai nol semua kecuali baris pertama M2M1Ax=M2M1b → Kolom KEDUA matrix M1A bernilai nol semua kecuali baris kedua M3M2M1Ax=M3M2M1b → Kolom KETIGA matrix M2M1A bernilai nol semua kecuali baris ketiga Computational Physics by Agus Naba, Ph.D.- UB

24

MAx=Mn-1…M3M2M1Ax=Mn-1...M3M2M1b M=Mn-1…M3M2M1 M-1=L U=MA -> A= M-1U A=LU Computational Physics by Agus Naba, Ph.D.- UB

25

LU Factorization A=LU Ax=b → LUx=b → x=? Forward substitution: Ly=b Barkward substitution: Ux=y Computational Physics by Agus Naba, Ph.D.- UB

26

Algorithm of LU Factorization for k=1 to n-1 if akk=0 then stop for i=k+1 to n mik=aik/akk end for j=k to n for i=k+1 to n aij=aij-mikakj end end end

{Loop over columns} {stop if pivot is zero} {compute multipliers for current column}

{apply transformation to remaining submatrix}


27

Example

x

1

+ 2x

+ 2x

=

+ 4x + 6x

+ 2x + 4x

= 6 = 10

2

4x 4x

1

1

2

2

3

3

3

3

Cari M1 dan M2 lalu temukan L dan U, kemudian pecahkan x1,x2,dan x3 ! Computational Physics by Agus Naba, Ph.D.- UB

28

Problem of LU Factorization (LUF) Metode LUF tidak bisa dipakai jika elemen diagonal matrix A bernilai nol/sangat kecil, meskipun A adalah nonsingular. Masalah ini diatasi dengan melakukan pivoting, yaitu menukar baris matrix yang elemen diagonalnya nol/sangat kecil dgn baris yang lain yang elemen diagonalnya tidak nol./besar Computational Physics by Agus Naba, Ph.D.- UB

29

Example 1 0 1  A=  1 0  

non-singular BUT no LU factorization

1 1 1 0 1 1 A= = = LU     1 1 1 1 0 0 non-singular and has LU factorization


30

Example 2 ε 1 A=  1 1  

0<ε< εmach

0  1  1 0 M= ; L =    − 1 / ε 1 1 / ε 1 1  ε 1  ε U= =   0 1 − 1 / ε  0 − 1 / ε  Computationalarithmetic Physics In floating-point ! by Agus Naba, Ph.D.- UB

31

1  ε 1   1 0  ε LU =  = ≠A     1 / ε 1 0 − 1 / ε  1 0


32

Example 2 (contd.) ε 1 A=  1 1  

1 1 A =  ε 1   c

 1 0 1 0 M= ; L =    − ε 1 ε 1  1  1 1  1 1 1 U= = =    0 1 − ε  0 1 − ε  0 1 Computationalarithmetic Physics In floating-point ! by Agus Naba, Ph.D.- UB

33

1 0 1 1 1 1 LU =  = =A     ε 1 0 1 ε 1


c

34

Ax=b

1 2 2   x   3   4 4 2  x  =  6       4 6 2  x  10 1

2

3


35

0 1 0  P = 1 0 0   0 0 1 1

 4 4 2  x   6  P Ax = 1 2 2  x  =  3  = P b      4 6 4  x  10 1

1

2

1

3


36

0 1 0  M = 1 0 0   0 0 1 1

4 4 2   x   6  M P Ax = 0 1 1.5  x  = 1.5 = M P b      0 2 2   x   4  1

1

1

2

1

1

3


37

1 0 0  P = 0 0 1    0 1 0 2

4 4 2   x   6  P M P Ax = 0 2 2   x  =  4  = P M P b      0 1 1.5  x  1.5 1

2

1

1

2

2

1

1

3


38

0 1 0  M = 1 0 0   0 0 1 2

4 4 2   x   6  M P M P Ax = 0 2 2   x  =  4  = M P M P b      0 0 0.5  x  − 0.5 1

2

2

1

1

2

2

2

1

1

3


39

L = M = (M P M P ) = P L P L −1

−1

2

2

1

1

T

1

T

1

2

2

0.25 0.5 1   = 1 0 0   1 0  1


40

0.25 0.5 1 4 4 2      LU = 1 0 0 0 2 2 =A    1 0 0 0 0.5  1


41

LU Factorization by Gaussian Elimination with Partial Pivoting for k=1 to n-1 find index p such that |apk| > |aik| for k ≤ i ≤ n if p ≠ k then interchange rows k and p if akk=0 then continue with next k for i=k+1 to n mik=aik/akk end for j=k to n for i=k+1 to n aij=aij-mikakj end end end

{Loop over columns} {search for pivot current column} {interchange rows, if necessary} {skip current column if it’s zero already} {compute multipliers for current column}

{apply transformation to remaining submatrix}


42

Cholesky Factorization for k=1 to n akk = sqrt(akk) for i=k+1 to n aik=aik/akk end for j=k+1 to n for i=k+1 to n aij=aij-aikajk end end end

{Loop over columns}


43

Banded System d1=b1 for i=2 to n mi=ai/di-1 di=bi-mici-1 end


44


45

Linear Equation. A b x. Computational Physics by Agus Naba, Ph.D.- UB

Recommend Documents