(Pertemuan 9) Sistem Persamaan Linear Matriks: LU Factorization, Jacobi, Gauss-Seidel

LU Factorization method, Jacobi method, Gauss-Seidel method, SOR method

Offline di Departemen Matematika

Author

Aslab Mata Kuliah Praktikum Metode Numerik

Published

May 19, 2025

Kembali ke Metode Numerik

Outline:

(Pengayaan) Faktorisasi LU
(Pengayaan) Metode Jacobi, dalam bentuk matriks (teoritis)
(Pengayaan) Metode Gauss-Seidel, dalam bentuk matriks (teoritis)
(Pengayaan) Metode SOR, dalam bentuk matriks (teoritis)
Diskusi

Praktikum Metode Numerik PTA 2024-2025
Departemen Matematika FMIPA Universitas Indonesia

Modul praktikum ini diawali dengan pembahasan tentang beberapa jenis norm vektor. Kemudian, metode yang dibahas di modul kali ini, utamanya hanyalah “versi praktis” untuk algoritma metode Jacobi, metode Gauss-Seidel, dan metode relaksasi (SOR). Metode Gauss-Seidel adalah perbaikan dari metode Jacobi, dan metode SOR adalah perbaikan dari metode Gauss-Seidel. Istilah “versi praktis” di sini maksudnya agar dibedakan dari bentuk matriks \(T\textbf{x}+\textbf{c}\) (sebagai materi pengayaan) untuk metode-metode tersebut.

Inti sari dari ketiga metode tersebut adalah perumuman dari metode fixed-point (dari bab 2, metode numerik untuk root-finding), yang tadinya dilakukan untuk satu variabel/persamaan saja, menjadi dilakukan untuk beberapa variabel/persamaan sekaligus (yang kebetulan membentuk SPL). Langkah paling pertama dalam mempersiapkan penyelesaian SPL dengan metode-metode tersebut adalah seperti melakukan pindah ruas ke sebelah kanan untuk semua suku kecuali variabel pada diagonal, mirip dengan ide substitusi balik. Langkah ini tersirat ketika menuliskan bentuk umum metode-metode tersebut dalam bentuk sumasi. Selain itu, seperti metode fixed-point, diperlukan tebakan awal untuk nilai tiap variabel.

Untuk perumuman metode fixed-point yang lebih umum lagi, yaitu untuk sistem persamaan yang tidak harus linier (tidak harus berbentuk SPL), dibahas di bab 10.1 di buku Burden. Bab 8, 9, dan 10 di buku Burden dibahas di mata kuliah pilihan program studi S1 Matematika yang bernama “Matematika Numerik”, dengan prasyarat Metode Numerik.

Pembahasan teoritis di kelas (perkuliahan) Metode Numerik juga mencakup pembahasan metode Jacobi, metode Gauss-Seidel, dan metode SOR dalam bentuk matriks, dengan bentuk umum \(T\textbf{x}+\textbf{c}\). Bentuk matriks untuk metode-metode tersebut tidak menjadi fokus di praktikum (bahkan di buku Burden, akhir halaman 452, juga disebut bahwa bentuk matriks tersebut biasanya hanya digunakan untuk pembahasan teoritis), tetapi tetap disediakan di modul praktikum ini sebagai materi pengayaan.

1. (Materi pengayaan) Faktorisasi LU

Untuk mengurangi banyaknya operasi pada penyelesaian SPL dengan matriks (serta untuk beberapa alasan lainnya), faktorisasi matriks seringkali dilakukan. Ada bermacam-macam faktorisasi matriks, namun yang paling umum digunakan adalah faktorisasi LU (juga disebut dekomposisi LU). Pada faktorisasi LU, matriks A ditulis ulang (difaktorisasi) sebagai perkalian (bukan penjumlahan) antara matriks segitiga bawah L (lower triangular) dan matriks segitiga atas U (upper triangular):

\[A = LU\]

Ada tiga metode yang paing sering digunakan untuk faktorisasi LU, yaitu 1. Metode Doolittle 2. Metode Crout 3. Metode Cholesky

Perbedaan di antara ketiga metode tersebut adalah pada bentuk matriks \(L\) dan \(U\) yang akan diperoleh, lebih tepatnya pada ketentuan untuk elemen diagonalnya akan seperti apa.

Pada bab 6.5 di buku Burden, dibahas metode Doolittle, di mana faktorisasi LU dilakukan dengan menggunakan eliminasi Gauss (sedangkan metode Cholesky dan metode Crout dibahas di bab 6.6, algoritma 6.6 dan 6.7). Berikut ini, kita hanya membahas metode Doolittle.

Jika eliminasi Gauss dapat dilakukan pada sistem \(A\overrightarrow{x}=\overrightarrow{b}\) tanpa melakukan pertukaran baris, maka \(A=LU\), di mana \(m_{ji} = \frac{a_{ji}^{(i)}}{a_{ii}^{(i)}}\),

\(L = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ m_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ m_{n1} & m_{n2} & \cdots & 1 \end{bmatrix}\)

\(U = \begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn}^{(n)} \end{bmatrix}\)

(Fun fact: perhatikan bahwa, untuk metode Doolittle, semua elemen diagonal matriks \(L\) adalah 1, sedangkan elemen diagonal matriks \(U\) tidak harus satu. Untuk metode Crout, terbalik: semua elemen diagonal matriks \(U\) harus 1, sedangkan elemen diagonal matriks \(L\) boleh selain 1.)

Implementasi Faktorisasi LU dengan Metode Doolittle

import numpy as np
matrix = np.array(eval(input('Masukkan matriks yang akan difaktorisasi: ')))

def LUFactorization(input_matrix):
    matrix = input_matrix.astype(float)

    n = np.shape(matrix)[0] #mengambil ukuran baris dari matriks
    L = np.identity(n) #mendefinisikan L sebagai matriks identitas nxn
    #operasi baris elementer
    for i in range(n):
        for j in range(i+1, n):
            m = matrix[j,i]/matrix[i,i]
            L[j,i] = m #Pasang elemen L_ji menjadi multiplisitas m = a_ji/a_ii
            matrix[j]= matrix[j]-m*matrix[i]
    return (L, matrix)

L = LUFactorization(matrix)[0] #mengambil L pada LUFactorization
U = LUFactorization(matrix)[1] #mengambil matrix pada LUFactorization

print("faktorisasi LU matriksnya adalah :")
print("L = \n{0}".format(L)) #print L
#print U
print("U = \n{0}".format(U))

print("Apabila dikalikan, hasilnya menjadi:")
LU = np.matmul(L,U) #Hasil perkalian L dan U
print("LU = \n{0}".format(LU))

Masukkan matriks yang akan difaktorisasi: [[1,2,3],[4,5,6],[7,8,9]]
faktorisasi LU matriksnya adalah :
L = 
[[1. 0. 0.]
 [4. 1. 0.]
 [7. 2. 1.]]
U = 
[[ 1.  2.  3.]
 [ 0. -3. -6.]
 [ 0.  0.  0.]]
Apabila dikalikan, hasilnya menjadi:
LU = 
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]

2. (Pengayaan) Metode Jacobi, dalam bentuk matriks (teoritis)

Secara konsep, metode iteratif untuk SPL bisa dianggap sebagai semacam perumuman dari metode fixed-point, yang tadinya hanya satu variabel/persamaan menjadi banyak variabel/persamaan. Bentuk sumasi untuk masing-masing metode memang terlihat agak berbeda satu sama lain (seperti tidak bisa disamakan atau dibuat bentuk umumnya), terutama antara metode Jacobi dan metode Gauss-Seidel. Namun, mengingat asal-usulnya sebagai perumuman metode fixed-point, dan berhubung sistem persamaan yang ingin diselesaikan bersifat linier (membentuk SPL), metode iteratif untuk SPL bisa dituliskan dalam suatu bentuk umum menggunakan matriks (bentuk matriks), yakni

\[\textbf{x}^{(k)} = T\textbf{x}^{(k-1)} + \textbf{c}\]

di mana isi matriks \(T\) dan vektor konstanta \(\textbf{c}\) ditentukan tergantung metode iteratif yang digunakan: apakah metode Jacobi, metode Gauss-Seidel, atau metode SOR.

Sekilas, bentuk umum tersebut memang terlihat lebih sederhana, seperti betapa sederhananya metode fixed-point. Namun, secara perhitungan, perkalian matriks bisa memakan waktu yang relatif lama, sehingga versi praktis yang telah dibahas sebelumnya lah yang lebih cocok untuk dibuat program maupun untuk perhitungan manual.

Meskipun demikian, bentuk umum di atas masih ada kegunaannya, khususnya untuk mempermudah pembahasan teoritis seperti analisis error. Berikut ini, kita tetap akan membahas bentuk matriks untuk ketiga metode tersebut sebagai materi pengayaan.

Sebelumnya, dari SPL \(A\textbf{x}=\textbf{b}\), kita bisa “memecah” matriks koefisien \(A\) menjadi tiga bagian, yaitu \(A = (-L_{neg}) + D + (-U_{neg})\) atau sama saja \(A = D - L_{neg} - U_{neg}\): * Matriks \((-L_{neg})\) adalah matriks segitiga bawah menggunakan elemen matriks \(A\) yang berada di bawah diagonal, sedangkan sisanya nol. * Matriks \(D\) adalah matriks diagonal yang menggunakan elemen diagonal matriks \(A\), sedangkan sisanya nol. * Matriks \((-U_{neg})\) adalah matriks segitiga atas yang menggunakan elemen \(A\) yang berada di atas diagonal, sedangkan sisanya nol.

Perhatikan bahwa matriks \((-L_{neg})\) dan \((-U_{neg})\) dituliskan dengan tanda minus. Sebenarnya, nilai elemen segitiga bawah/atas yang disimpan ke matriks \(L_{neg}\) dan \(U_{neg}\) ini adalah negatif dari nilai aslinya di matriks \(A\). Sehinga, matriks segitiga bawah/atas yang memuat nilai aslinya bisa ditulis dengan minus: \((-L_{neg})\) dan \((-U_{neg})\). Keterangan “neg” maksudnya negatif, sehingga minus negatif menjadi kembali positif atau menjadi nilai aslinya. Hati-hati, pembahasan di buku Burden tidak melibatkan keterangan “neg”, sehingga langsung ditulis misalnya \(A=D-L-U\).

\[A = \begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{bmatrix} = \begin{bmatrix} a_{11} & 0 & \dots & 0 \\ 0 & a_{22} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & a_{nn} \end{bmatrix} - \begin{bmatrix} 0 & 0 & \dots & 0 \\ -a_{21} & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -a_{n1} & -a_{n2} & \dots & 0 \end{bmatrix} - \begin{bmatrix} 0 & -a_{12} & \dots & -a_{1n} \\ 0 & 0 & \dots & -a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 \end{bmatrix} \\ A = D - L_{neg} - U_{neg} \]

Dengan demikian, kita bisa menyusun fungsi untuk memisahkan matriks koefisien \(A\) menjadi \(D - L_{neg} - U_{neg}\).

def PisahDLnegUneg(matriks_A):
    # memperoleh ukuran n x n dari matriks A, ambil banyaknya baris aja
    n = np.shape(matriks_A)[0]

    # buat dulu matriks D, Lneg dan Uneg, ukuran n x n, sementara nol semua
    D = np.zeros((n,n))
    Lneg = np.zeros((n,n))
    Uneg = np.zeros((n,n))

    # double for loop melihat tiap elemen di matriks A...
    for i in range(n): # baris ke-i
        for j in range(n): # kolom ke-j
            if i == j: # jika elemen diagonal...
                # ... maka simpan ke matriks D
                D[i, j] = matriks_A[i, j]
            elif i > j: # jika lebih ke bawah daripada ke kanan...
                # ... maka simpan ke matriks Lneg (karena segitiga bawah)
                Lneg[i, j] = -matriks_A[i, j] # (jangan lupa dibuat negatif)
            else: # selain itu (berarti segitiga atas)
                # simpan ke matriks Uneg, jangan lupa dibuat negatif
                Uneg[i, j] = -matriks_A[i, j]
    
    # return tiga matriks sekaligus sebagai satu kesatuan
    return (D, Lneg, Uneg)

# Contoh
matriks_koef = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

D, Lneg, Uneg = PisahDLnegUneg(matriks_koef)
print(D)
print(Lneg)
print(Uneg)

[[1. 0. 0.]
 [0. 5. 0.]
 [0. 0. 9.]]
[[ 0.  0.  0.]
 [-4.  0.  0.]
 [-7. -8.  0.]]
[[ 0. -2. -3.]
 [ 0.  0. -6.]
 [ 0.  0.  0.]]

Selanjutnya, kita bisa menuliskan matriks \(T_j\) dan vektor konstanta \(\textbf{c}_j\) untuk metode Jacobi sebagai berikut:

\[T_j = D^{-1}\left(L_{neg}+U_{neg}\right), \hspace{0.5cm} \textbf{c}_j = D^{-1}\textbf{b}\]

sehingga rumus iterasi metode Jacobi menjadi

\[\textbf{x}^{(k)} = T_j\textbf{x}^{(k-1)} + \textbf{c}_j\]

def JacobiTeoritis(matriks_koefisien, vektor_b, tebakan_awal, tol):
    # pisahkan dulu
    D, Lneg, Uneg = PisahDLnegUneg(matriks_koefisien)

    # susun matriks T_j dan vektor konstanta c_j
    D_invers = np.linalg.inv(D)
    Tj = np.matmul( D_invers, Lneg+Uneg )
    cj = np.matmul( D_invers, vektor_b )

    # iterasi pertama

    # x^(k-1), salin dari tebakan awal
    xk_1 = tebakan_awal.copy()

    # x^(k), rumus metode Jacobi bentuk matriks
    xk = np.matmul( Tj, xk_1 ) + cj

    # iterasi kedua dan seterusnya dalam while loop

    while np.linalg.norm(xk_1 - xk, np.inf) > tol: # kriteria pemberhentian
        # yang sebelumnya menjadi x^(k) itu sekarang menjadi x^(k-1)
        xk_1 = xk

        # lakukan iterasi untuk memperoleh x^(k) yang baru
        xk = np.matmul( Tj, xk_1 ) + cj

    # jika sudah keluar while loop, toleransi sudah terpenuhi
    return xk

matriks_koef = np.array(eval(input("Masukkan matriks koefisien A: "))).astype(float)
vektor_b = np.array(eval(input("Masukkan vektor b: "))).astype(float)
tebakan_awal = np.array(eval(input("Masukkan tebakan awal: ")), dtype=float)

toleransi = eval(input("Masukkan toleransi: "))

hasil = JacobiTeoritis(matriks_koef, vektor_b, tebakan_awal, toleransi)
print("Hasil metode Jacobi (teoritis) adalah:")
print(hasil)

Masukkan matriks koefisien A: [ [10, -1, 2, 0], [-1, 11, -1, 3], [2, -1, 10, -1], [0, 3, -1, 8] ]
Masukkan vektor b: [6, 25, -11, 15]
Masukkan tebakan awal: [0, 0, 0, 0]
Masukkan toleransi: 10**-4
Hasil metode Jacobi (teoritis) adalah:
[ 0.99998973  2.00001582 -1.00001257  1.00001924]

3. (Pengayaan) Metode Gauss-Seidel, dalam bentuk matriks (teoritis)

Untuk metode Gauss-Seidel, kita definisikan matriks \(T_g\) dan vektor konstanta \(\textbf{c}_g\) sebagai berikut:

\[T_g = \left( D - L_{neg} \right)^{-1} U_{neg}, \hspace{0.5cm} \textbf{c}_g = \left( D - L_{neg} \right)^{-1} \textbf{b}\]

Sehingga, rumus iterasi untuk metode Gauss-Seidel bentuk matriks bisa ditulis:

\[\textbf{x}^{(k)} = T_g \textbf{x}^{(k-1)} + \textbf{c}_g\]

def GaussSeidelTeoritis(matriks_koefisien, vektor_b, tebakan_awal, tol):
    # pisahkan dulu
    D, Lneg, Uneg = PisahDLnegUneg(matriks_koefisien)

    # susun matriks T_g dan vektor konstanta c_g
    DminusLneg_invers = np.linalg.inv(D - Lneg)
    Tg = np.matmul( DminusLneg_invers, Uneg )
    cg = np.matmul( DminusLneg_invers, vektor_b )

    # iterasi pertama

    # x^(k-1), salin dari tebakan awal
    xk_1 = tebakan_awal.copy()

    # x^(k), rumus metode Gauss-Seidel bentuk matriks
    xk = np.matmul( Tg, xk_1 ) + cg

    # iterasi kedua dan seterusnya dalam while loop

    while np.linalg.norm(xk_1 - xk, np.inf) > tol: # kriteria pemberhentian
        # yang sebelumnya menjadi x^(k) itu sekarang menjadi x^(k-1)
        xk_1 = xk

        # lakukan iterasi untuk memperoleh x^(k) yang baru
        xk = np.matmul( Tg, xk_1 ) + cg

    # jika sudah keluar while loop, toleransi sudah terpenuhi
    return xk

matriks_koef = np.array(eval(input("Masukkan matriks koefisien A: "))).astype(float)
vektor_b = np.array(eval(input("Masukkan vektor b: "))).astype(float)
tebakan_awal = np.array(eval(input("Masukkan tebakan awal: ")), dtype=float)

toleransi = eval(input("Masukkan toleransi: "))

hasil = GaussSeidelTeoritis(matriks_koef, vektor_b, tebakan_awal, toleransi)
print("Hasil metode Gauss-Seidel (teoritis) adalah:")
print(hasil)

Masukkan matriks koefisien A: [ [10, -1, 2, 0], [-1, 11, -1, 3], [2, -1, 10, -1], [0, 3, -1, 8] ]
Masukkan vektor b: [6, 25, -11, 15]
Masukkan tebakan awal: [0, 0, 0, 0]
Masukkan toleransi: 10**-4
Hasil metode Gauss-Seidel (teoritis) adalah:
[ 1.00000836  2.00000117 -1.00000275  0.99999922]

4. (Pengayaan) Metode SOR, dalam bentuk matriks (teoritis)

Untuk metode SOR, diberikan suatu nilai omega, kita definisikan matriks \(T_{\omega}\) dan vektor konstanta \(\textbf{c}_{\omega}\) sebagai berikut:

\[T_{\omega} = \left( D-\omega L \right)^{-1}\left[ (1-\omega)D + \omega U \right], \hspace{0.5cm} \textbf{c}_{\omega} = \omega \left( D-\omega L\right)^{-1} \textbf{b}\]

Sehingga, rumus iterasi untuk metode SOR bentuk matriks bisa ditulis:

\[\textbf{x}^{(k)} = T_{\omega} \textbf{x}^{(k-1)} + \textbf{c}_{\omega}\]

def SORTeoritis(matriks_koefisien, vektor_b, tebakan_awal, omega, tol):
    # pisahkan dulu
    D, Lneg, Uneg = PisahDLnegUneg(matriks_koefisien)

    # susun matriks T_omega dan vektor konstanta c_omega
    DminusomegaLneg_invers = np.linalg.inv( D - omega * Lneg)
    T_omega = np.matmul ( DminusomegaLneg_invers, (1-omega)*D + omega*Uneg )
    c_omega = omega * np.matmul( DminusomegaLneg_invers, vektor_b )

    # iterasi pertama

    # x^(k-1), salin dari tebakan awal
    xk_1 = tebakan_awal.copy()

    # x^(k), rumus metode SOR bentuk matriks
    xk = np.matmul( T_omega, xk_1 ) + c_omega

    # iterasi kedua dan seterusnya dalam while loop

    while np.linalg.norm(xk_1 - xk, np.inf) > tol: # kriteria pemberhentian
        # yang sebelumnya menjadi x^(k) itu sekarang menjadi x^(k-1)
        xk_1 = xk

        # lakukan iterasi untuk memperoleh x^(k) yang baru
        xk = np.matmul( T_omega, xk_1 ) + c_omega

    # jika sudah keluar while loop, toleransi sudah terpenuhi
    return xk

matriks_koef = np.array(eval(input("Masukkan matriks koefisien A: "))).astype(float)
vektor_b = np.array(eval(input("Masukkan vektor b: "))).astype(float)
tebakan_awal = np.array(eval(input("Masukkan tebakan awal: ")), dtype=float)
omega = eval(input("Masukkan faktor relaksasi (omega): "))
toleransi = eval(input("Masukkan toleransi: "))

hasil = SORTeoritis(matriks_koef, vektor_b, tebakan_awal, omega, toleransi)
print("Hasil metode SOR (teoritis) adalah:")
print(hasil)

Masukkan matriks koefisien A: [ [4, 3, 0], [3, 4, -1], [0, -1, 4] ]
Masukkan vektor b: [24, 30, -24]
Masukkan tebakan awal: [0, 0, 0]
Masukkan faktor relaksasi (omega): 1.25
Masukkan toleransi: 10**-4
Hasil metode SOR (teoritis) adalah:
[ 2.99998919  4.00000321 -4.9999937 ]