Modul 5: Metode Langsung untuk SPL

Modul 5: Metode Langsung untuk SPL

Kembali ke Metode Numerik

Outline

  1. Operasi matriks pada Python
  2. Review SPL: sistem persamaan linier (penjelasan tanpa kode)
  3. Memperoleh matriks diperbesar \(\tilde{A}\) dari \(A\textbf{x}=\textbf{b}\), dan sebaliknya
  4. Eliminasi Gauss dan substitusi balik
  5. Partial Pivoting
  6. Scaled Partial Pivoting
  7. (Materi pengayaan) Faktorisasi LU
import numpy as np

1. Operasi matriks pada Python

Sebelum masuk ke materi metode numerik untuk sistem persamaan linier (SPL), mari kita bahas lebih lanjut tentang operasi matriks menggunakan numpy di Python.

Mengingat kembali, tanpa numpy, matriks dalam Python bisa dituliskan sebagai list dua dimensi (list di dalam list).

matriks_manual = [ [1, 2, 3], [4, 5, 6] ]
print(matriks_manual)
[[1, 2, 3], [4, 5, 6]]

Ada beberapa keunggulan array numpy dibandingkan dengan list dua dimensi yang dibuat secara manual seperti itu. Cara membuatnya adalah memasukkan suatu list dua dimensi ke dalam np.array, seperti berikut:

# sebelumnya, sudah dibuat list dua dimensi bernama "matriks_manual"
matriks_numpy = np.array(matriks_manual)
print(matriks_numpy)
[[1 2 3]
 [4 5 6]]

Pada kode di atas, kita telah membuat list dua dimensi di variabel terpisah, sebelum memasukkannya di dalam np.array. Namun, tentu saja, kita bisa langsung membuat list dua dimensinya di dalam np.array:

matriks_baru = np.array([ [1,2,3], [4,5,6] ])
print(matriks_baru)
[[1 2 3]
 [4 5 6]]

Perhatikan bahwa tiap list di dalam list adalah baris pada matriks. Misalnya, ada list di dalam list, [1,2,3] yang menjadi baris pertama, diikuti dengan list di dalam list, [4,5,6] yang menjadi baris berikutnya. Kedua list tersebut merupakan bagian dari satu list besar (perhatikan, di dalam np.array itu diawali dan diakhiri kurung siku, karena sebenarnya np.array menerima input berupa list di dalam list).

Hal ini akan penting nantinya ketika ingin menerima input berupa matriks dari user (pengguna).

Sebenarnya, np.array bisa saja menerima input berupa list biasa (bisa dikatakan satu dimensi), di mana outputnya akan berupa array satu dimensi. Selain itu, numpy bisa membuat beberapa jenis array/matriks istimewa. Contohnya, array/matriks yang berisi angka nol semua, dengan np.zeros:

baris_nol = np.zeros(5) # lima elemen
print(baris_nol)
[0. 0. 0. 0. 0.]
matriks_nol = np.zeros( (3,2) ) # tiga baris, dua kolom
print(matriks_nol)
[[0. 0.]
 [0. 0.]
 [0. 0.]]

Perhatikan bahwa, untuk array berdimensi dua (matriks), ada kurung di dalam kurung (seolah-olah, input yang diterima adalah semacam “koordinat”), tidak seperti untuk array biasa (satu dimensi) yang langsung dimasukkan banyaknya elemen tanpa ada kurung lagi.

Selain nol semua, numpy juga bisa membuat array/matriks yang berisi angka 1 semua, dengan cara yang serupa, dengan np.ones.

print(np.ones(4))
[1. 1. 1. 1.]
print(np.ones( (2, 5) ))
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]

Untuk angka selain nol dan satu, kita tinggal membuat array/matriks yang berisi satu semua, kemudian dikalikan dengan apapun angka itu.

# matriks berisi 7 semua
print(7 * np.ones( (2, 5) ))
[[7. 7. 7. 7. 7.]
 [7. 7. 7. 7. 7.]]

Kemudian, kita juga bisa membuat matriks diagonal (yang tentunya merupakan matriks persegi), dengan elemen diagonal sesuai yang kita inginkan, menggunakan np.diag.

elemen_diagonal = np.array([5, 4, 3, 2])
print(np.diag(elemen_diagonal))
[[5 0 0 0]
 [0 4 0 0]
 [0 0 3 0]
 [0 0 0 2]]

Artinya, untuk membuat matriks identitas, kita bisa menerapkan np.diag pada np.ones.

print(np.diag(np.ones(4)))
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

Sebenarnya, dari numpy sudah ada fungsi khusus untuk membuat matriks identitas, yaitu np.identity.

print(np.identity(3))
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Sama seperti list, pada matriks juga dapat dilakukan indexing dan slicing. Indexing pada matriks juga dimulai dari 0. Matriks adalah array 2-D, sehingga indeks akan terdiri dari [i, j] di mana i menyatakan indeks baris dan j menyatakan indeks kolom.

A = np.array([[1,2,3], [4,5,6]]) #mendefinisikan matriks A 2x3
B = np.array([[-1,0,1], #mendefinisikan matriks B 2x3
              [0,0,1]])
C = np.array([[1,0,1], #mendefinisikan matriks B 3x3
              [0,1,1],
              [1,1,1]])
print(A)
[[1 2 3]
 [4 5 6]]
print(B)
[[-1  0  1]
 [ 0  0  1]]
print(C)
[[1 0 1]
 [0 1 1]
 [1 1 1]]
print(A[0]) #menampilkan baris pertama (indeks 0) dari matriks A
[1 2 3]
#menampilkan baris pertama (indeks 0), kolom kedua (indeks 1) dari matriks A
print(A[0, 1])
2
# tampilkan baris pertama (indeks 0),
# mulai dari kolom kedua (indeks 1) dan seterusnya
print(A[0, 1:])
[2 3]
# tampilkan baris pertama (indeks 0),
# tampilkan semua kolom sampai sebelum kolom ketiga (sebelum indeks 2)
print(A[0, :2])
[1 2]
# tampilkan nilai pada semua baris,
# tapi melihat kolom kedua (indeks 1) saja
print(A[:, 1])
[2 5]
# tampilkan semua baris,
# mulai dari kolom kedua (indeks 1) dan seterusnya
print(A[:, 1:])
[[2 3]
 [5 6]]
# tampilkan baris pertama (indeks nol),
# tapi kolom pertama dari belakang (hitung mundur)
print(A[0, -1])
3
# tampilkan baris pertama dari belakang,
# kolom pertama dari belakang
print(A[-1, -1])
6

Operasi dasar seperti penjumlahan dan pengurangan dapat dilakukan secara langsung seperti halnya penjumlahan/pengurangan bilangan.

print(A+B)
print(A-B)
[[0 2 4]
 [4 5 7]]
[[2 2 2]
 [4 5 5]]

Operasi perkalian skalar dapat menggunakan tanda bintang atau asterisk (*), dan urutannya boleh ditukar.

print(3*A) # 3 dikali A
print(B*4) # B dikali 4
[[ 3  6  9]
 [12 15 18]]
[[-4  0  4]
 [ 0  0  4]]

Apabila dua matriks dikalikan begitu saja dengan tanda bintang, maka perkalian akan dilakukan secara broadcasting, yaitu per elemen.

print(A)
print(B)
print(A*B)
[[1 2 3]
 [4 5 6]]
[[-1  0  1]
 [ 0  0  1]]
[[-1  0  3]
 [ 0  0  6]]

Perkalian matriks yang biasa kita kenal di aljabar linier tidak seperti itu. Numpy menyediakan fungsi khusus untuk perkalian matriks yang seperti di aljabar linier, yaitu np.matmul (matrix multiplication). Tentu saja, ada syarat ukuran matriks, yaitu \(m \times n\) dan \(n \times p\).

Kode berikut ini akan gagal karena tidak memenuhi syarat.

np.matmul(A,B)
ValueError: ignored

Perkalian A dengan B tidak dapat dilakukan dan muncul error message. Cek ukuran dari matriks dengan menggunakan np.shape.

print(np.shape(A)) #Ukuran matriks A
print(np.shape(B)) #Ukuran matriks B
print(np.shape(C)) #Ukuran matriks C
(2, 3)
(2, 3)
(3, 3)

Baik A dan B memiliki ukuran 2x3, sehingga AB tidak terdefinisi. Namun, apabila kita men-transpose B, kita dapat melakukan perkalian \(AB^T\). Untuk mentranspose matriks, gunakan np.transpose

np.matmul(A, np.transpose(B)) #A B^T
array([[2, 3],
       [2, 6]])

Sebagai tambahan, numpy juga bisa menghitung dot product (perkalian dot, yaitu hasil kali titik) antara dua array satu dimensi, menggunakan np.dot

vektor1 = np.array([1, -5, 0])
vektor2 = np.array([-3, 7, 10])
print(np.dot(vektor1, vektor2))
-38

Seandainya kita menggunakan np.dot dengan dua matriks, maka numpy akan mengartikannya sebagai np.matmul

np.dot(A, np.transpose(B)) #A B^T
array([[2, 3],
       [2, 6]])

Selebihnya bisa dibaca di dokumentasi numpy:

https://numpy.org/doc/stable/reference/generated/numpy.dot.html

Terakhir, numpy memiliki beberapa fungsi khusus lainnya untuk aljabar linier, yang menariknya mengharuskan penulisan “linalg” (linear algebra; aljabar linier), karena memang merupakan bagian khusus di dalam numpy. Contohnya adalah determinan dan invers.

D = np.array([[2, -3], [-2, 5]])

print(D)
print(np.linalg.det(D)) # det(D), yaitu determinan dari matriks D
print(np.linalg.inv(D)) # D^-1, yaitu invers dari matriks D

print(np.linalg.det(np.linalg.inv(D))) # det(D^-1)
[[ 2 -3]
 [-2  5]]
4.0
[[1.25 0.75]
 [0.5  0.5 ]]
0.24999999999999994

Jangan lupa, apabila ada hasil yang sedikit aneh, seperti 1/4 = 0.2499999…, itu disebabkan oleh kelemahan floating-point precision yang dibahas di pertemuan pertama kuliah Metode Numerik. Python tidak kebal terhadap masalah tersebut.

Selain itu, apabila keseluruhan matriks berisi bilangan bulat, bisa saja dilakukan integer division, di mana semua hasil pembagian itu dibulatkan ke bawah. Hal ini tentu sangat berbahaya jika ada operasi pembagian dalam metode numerik. Untuk menghindari masalah tersebut, array bisa dikonversi menjadi float semua, menggunakan .astype(float)

# berisi bilangan bulat semua
arraybulat = np.array([5, 4])
print(arraybulat)

# arraybulat[0] = 5//4 = floor(5/4) = floor(1.25) = 1
arraybulat[0] = arraybulat[0]/arraybulat[1]
print(arraybulat)
[5 4]
[1 4]
arraybulat = np.array([5, 4])
arrayfloat = arraybulat.astype(float)
print(arrayfloat)

# mencoba hal yang sama,
# kali ini tidak ada integer division sehingga 5/4 = 1.25
arrayfloat[0] = arrayfloat[0]/arrayfloat[1]
print(arrayfloat)
[5. 4.]
[1.25 4.  ]

2. Review SPL: sistem persamaan linier (penjelasan tanpa kode)

Suatu sistem persamaan linier (SPL) adalah kumpulan beberapa persamaan linier dalam beberapa variabel \(x_1, x_2, \dots, x_n\), misal sebanyak \(m\) persamaan. Idealnya, banyaknya variabel sama dengan banyaknya persamaan, yaitu \(n=m\). (Praktikum Metode Numerik akan membahas SPL dengan \(n=m\).)

Bentuk umum SPL bisa dituliskan sebagai berikut:

\[ \begin{align} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2n} x_n = b_2 \\ a_{31} x_1 + a_{32} x_2 + &\dots + a_{3n} x_n = b_3 \\ &\vdots \\ a_{m1} x_1 + a_{m2} x_2 + &\dots + a_{mn} x_n = b_m \end{align} \]

di mana koefisien \(a_{ij}\) adalah koefisien pada persamaan ke-i untuk variabel \(x_j\), dan ada konstanta \(b_i\) untuk tiap persamaan \(i = 1, 2, \dots, m\).

Umumnya, semua koefisien \(a_{ij}\) serta konstanta \(b_i\) sudah diketahui nilainya, dan ingin dicari nilai-nilai \(x_j\) yang bersama memenuhi semua persamaan sekaligus, disebut solusi dari SPL tersebut.

Apabila \(n=m\), bentuk umum SPL menjadi

\[ \begin{align} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2n} x_n = b_2 \\ a_{31} x_1 + a_{32} x_2 + &\dots + a_{3n} x_n = b_3 \\ &\vdots \\ a_{n1} x_1 + a_{n2} x_2 + &\dots + a_{nn} x_n = b_n \end{align} \]

yang dapat dituliskan dalam bentuk perkalian matriks-vektor:

\[ \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ a_{31} & a_{32} & \dots & a_{3n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ b_3 \\ \vdots \\ b_n \end{pmatrix} \]

Matriks koefisien \(a_{ij}\) bisa ditulis \(A\), vektor kolom \(x_j\) bisa ditulis \(\textbf{x}\), dan vektor kolom \(b_i\) bisa ditulis \(\textbf{b}\), agar bentuk perkalian matriks-vektor di atas bisa diringkas: \(A\textbf{x}=\textbf{b}\). Dalam hal ini, \(\textbf{x}\) adalah vektor solusi.

Dengan demikian, notasi \(a_{ij}\) bisa juga diartikan sebagai elemen matriks \(A\) pada baris ke-i, kolom ke-j.

Kita bisa “menggabungkan” vektor \(\textbf{b}\) menjadi kolom baru (kolom paling kanan) di matriks \(A\), sehingga dari yang tadinya berukuran \(n \times n\) menjadi berukuran \(n \times \left(n+1\right)\):

\[ \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} & b_1 \\ a_{21} & a_{22} & \dots & a_{2n} & b_2 \\ a_{31} & a_{32} & \dots & a_{3n} & b_3 \\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1} & a_{n2} & \dots & a_{nn} & b_n \end{pmatrix} \]

Matriks baru ini biasa disebut augmented matrix atau matriks diperbesar, dan biasa ditulis \(\tilde{A}\).

Kita juga bisa menuliskan \(a_{i,\left(n+1\right)} = b_i\) untuk \(i = 1, 2, \dots, n\), agar konsisten dengan notasi \(a_{ij}\) yaitu elemen matriks pada baris ke-i, kolom ke-j.

\[ \tilde{A} = \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} & a_{1,\left(n+1\right)} \\ a_{21} & a_{22} & \dots & a_{2n} & a_{2,\left(n+1\right)} \\ a_{31} & a_{32} & \dots & a_{3n} & a_{3,\left(n+1\right)} \\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1} & a_{n2} & \dots & a_{nn} & a_{n,\left(n+1\right)} \end{pmatrix} \]

3. Memperoleh matriks diperbesar dari \(A\textbf{x}=\textbf{b}\), dan sebaliknya

Misalkan kita punya SPL seperti berikut:

\[ \begin{array}{rcrcrcrc} x_1 & - & x_2 & + & 2x_3 & - & x_4 & = & -8 \\ 2x_1 & - & 2x_2 & + & 3x_3 & - & 3x_4 & = & -20 \\ x_1 & + & x_2 & + & x_3 & & & = & -2 \\ x_1 & - & x_2 & + & 4x_3 & + & 3x_4 & = & 4 \end{array} \]

Matriks koefisien \(A\) dan vektor \(\textbf{b}\) dari SPL di atas adalah

\[ A = \begin{pmatrix} 1 & -1 & 2 & -1 \\ 2 & -2 & 3 & -3 \\ 1 & 1 & 1 & 0 \\ 1 & -1 & 4 & 3 \end{pmatrix}, \hspace{0.5cm} \textbf{b} = \begin{pmatrix} -8 \\ -20 \\ -2 \\ 4 \end{pmatrix}\]

Sehingga, matriks diperbesar \(\tilde{A}\) dari SPL di atas adalah

\[ \tilde{A} = \begin{pmatrix} 1 & -1 & 2 & -1 & -8 \\ 2 & -2 & 3 & -3 & -20 \\ 1 & 1 & 1 & 0 & -2 \\ 1 & -1 & 4 & 3 & 4 \end{pmatrix} \]

Apabila kita hanya memiliki matriks koefisien \(A\) dan vektor \(\textbf{b}\), kita bisa saja memperoleh matriks diperbesar \(\tilde{A}\) menggunakan numpy:

# [[1,-1,2,-1, -8],[2,-2,3,-3,-20],[1,1,1,0,-2], [1,-1,4,3,4]]

A_koef = np.array([
    [1,-1,2,-1],
    [2,-2,3,-3],
    [1,1,1,0],
    [1,-1,4,3]
])

# buat vektor baris dua dimensi, lalu transpos
# karena diperlukan vektor kolom dua dimensi
b = np.transpose(np.array([[-8,-20,-2,4]]))

# matriks diperbesar diperoleh dengan meletakan vektor b "di sebelah kanan" A,
# atau sama saja "ditumpuk" secara horizontal, bisa dengan numpy.hstack
A_diperbesar = np.hstack((A_koef, b))

print(A_diperbesar)
[[  1  -1   2  -1  -8]
 [  2  -2   3  -3 -20]
 [  1   1   1   0  -2]
 [  1  -1   4   3   4]]

Tentunya, kita juga bisa memperoleh matriks koefisien \(A\) dan vektor \(\textbf{b}\) secara pemrograman apabila hanya diketahui matriks diperbesar \(\tilde{A}\).

A_diperbesar = np.array([
    [1,-1,2,-1,-8],
    [2,-2,3,-3,-20],
    [1,1,1,0,-2],
    [1,-1,4,3,4]
])

# n adalah banyaknya baris
n = np.shape(A_diperbesar)[0]

# gunakan nilai-nilai di semua baris sampai kolom ke-n untuk matriks koefisien
A_koef = A_diperbesar[:, :n]

# untuk vektor b, peroleh nilai pada kolom terakhir
b = A_diperbesar[:, n]

print("Matriks koefisien A:")
print(A_koef)
print("Vektor b:")
print(b)
Matriks koefisien A:
[[ 1 -1  2 -1]
 [ 2 -2  3 -3]
 [ 1  1  1  0]
 [ 1 -1  4  3]]
Vektor b:
[ -8 -20  -2   4]

4. Eliminasi Gauss dan substitusi balik

Seperti yang sudah dipelajari di aljabar linier, eliminasi Gauss adalah teknik menyelesaikan (SPL) dengan menerapkan operasi baris elementer (OBE) pada matriks diperbesarnya. OBE adalah beberapa operasi khusus bisa yang dilakukan pada satu/dua baris dalam suatu matriks diperbesar (misal baris \(E_i\) dan \(E_j\)), dengan sifat istimewa yaitu tidak akan mengubah nilai solusi SPL \(x_1, x_2, \dots, x_n\) sama sekali. OBE bisa berupa:

  1. Pertukaran baris: \((E_i) \leftrightarrow (E_j)\)
  2. Perkalian baris oleh skalar: \((\lambda E_i) \rightarrow (E_i)\)
  3. Penjumlahan suatu baris dengan kelipatan skalar dari baris lain \(( E_i + \lambda E_j ) \rightarrow (E_i)\)

(“E” adalah singkatan dari equation.)

Misalkan terdapat SPL yang dinyatakan dalam bentuk \(Ax = b\), di mana A adalah matriks berukuran \(n \times n\) dan \(\textbf{b}\) adalah vektor berukuran \(n \times 1\). Eliminasi Gauss bertujuan untuk mengubah SPL awal menjadi bentuk triangular:

\[ \begin{align} a_{11} x_1 + a_{12} x_2 + \dots + a_{1n}x_n &= a_{1,n+1} \\ a_{22} x_2 + \dots + a_{2n}x_n &= a_{2,n+1} \\ \vdots \\ a_{nn}x_n &= a_{n,n+1} \end{align} \]

Sebenarnya, semua koefisien yang terlihat “hilang” itu masih ada, hanya saja sudah berhasil diubah menjadi nol.

Kemudian, untuk mencari nilai \(x_1, x_2, \dots, x_n\) dari bentuk triangular tersebut, lakukan substitusi balik (back substitution). Dari persamaan terakhir diperoleh

\[x_n = \frac{a_{n,n+1}}{a_{nn}}\]

Substitusi \(x_n\) ke persamaan ke-(n-1) diperoleh \(x_{n-1}\). Substitusi \(x_n\) dari \(x_{n-1}\) ke persamaan ke-(n-2) diperoleh \(x_{n-2}\). Lakukan terus sampai mendapatkan \(x_1\).

Untuk menyelesaikan SPL menggunakan eliminasi Gauss dan substitusi balik (Gaussian elimination with backward substitution), SPL dapat ditulis sebagai matriks diperbesar. Kemudian, lakukan langkah-langkah berikut.

  1. Buat semua entri di bawah \(a_{ii}\) (untuk setiap kolom \(i = 1, 2, \dots, n\)) menjadi nol dengan melakukan operasi baris elementer: \[\left( E_j - \left( \frac{a_{ji}}{a_{ii}} \right)E_i \right) \rightarrow E_j\] untuk baris ke-j, dengan \(j = i+1, i+2, \dots, n\). Namun, kita bisa menuliskan \(m = \frac{a_{ji}}{a_{ii}}\) (m: multiplier; pengkali) agar bentuknya menjadi \[\left( E_j - mE_i \right) \rightarrow E_j\] atau sama saja \[E_j \leftarrow \left( E_j - mE_i \right)\]
  2. Lakukan substitusi balik, diawali rumus: \[x_n = \frac{a_{n,n+1}}{a_{nn}}\] Kemudian, menghitung mundur, untuk \(i = n-1, n-2, \dots, 2, 1\), hitung: \[x_i = \frac{a_{i,n+1} - \sum_{j=i+1}^{n}a_{ij}x_j}{a_{ii}}\] Fun fact: rumus itu diperoleh dengan melakukan pindah ruas pada persamaan di bentuk triangularnya, seperti berikut: \[a_{ii}x_{ii} + a_{i,i+1}x_{i,i+1} + a_{i,i+2}x_{i,i+2} + \dots + a_{in}x_{in} = b_i = a_{i,n+1}\] \[a_{ii}x_{ii} + \sum_{j=i+1}^{n}a_{ij}x_j = a_{i,n+1}\] \[a_{ii}x_{ii} = a_{i,n+1} - \sum_{j=i+1}^{n}a_{ij}x_j\] \[x_i = \frac{a_{i,n+1} - \sum_{j=i+1}^{n}a_{ij}x_j}{a_{ii}}\]

Namun, ketika melakukan eliminasi Gauss, apabila ada elemen diagonal \(a_{ii}\) yang bernilai nol, maka baris yang mengandung \(a_{ii}\) perlu ditukar dengan baris di bawahnya yang elemennya taknol pada kolom yamg sama (kolom ke-i), agar elemen diagonal yang baru menjadi taknol.

Implementasi Eliminasi Gauss dan Substitusi Balik

def EliminasiGauss(matriks_input):
    # konversi matriks_input jadi matriks baru yang isinya float semua,
    # karena apabila ada bilangan bulat, bisa jadi dilakukan integer division
    # yang bisa sangat memperparah error
    matriks = matriks_input.astype(float)

    # memperoleh ukuran baris dari matriks diperbesar
    n = np.shape(matriks)[0]
    # Ingat bahwa ukuran matriks diperbesar adalah n x (n+1)

    for i in range (n): # untuk kolom ke-i (dari kolom awal sampai ke-n)
        # Saat ini, kita sedang melakukan eliminasi Gauss untuk kolom ke-i.
        # Semua nilai koefisien di bawah elemen diagonal akan dibuat nol

        # Sebelum mengeliminasi, kita perlu memastikan elemen diagonal taknol.
        # Kalau misalnya nol, kita perlu melihat baris-baris berikutnya
        # untuk bertukar baris agar elemen diagonal menjadi taknol

        # Variabel p ("pivot") akan digunakan untuk melihat baris.
        # Kita lihat dulu baris ke-i
        p = i
        # sehingga, saat ini, matriks[p,i] adalah elemen diagonal.
        # Ingat, elemen diagonal harusnya taknol.

        # Kalau ternyata nilai elemen tersebut adalah nol,
        # lanjut melihat di bawahnya (mencari calon baris yang bisa ditukar),
        # dan kalau masih nol, lihat ke bawahnya lagi, dan seterusnya
        while p<n and matriks[p,i]==0:
            p += 1
        # tapi jangan sampai keluar dari matriks (melewati baris terakhir),
        # makanya dibuat syarat p<n
        
        # Kalau sudah keluar dari matriks, artinya semua elemen di bawah
        # diagonal, bahkan termasuk elemen diagonal, itu nol semua.
        # Sayangnya, SPL tidak bisa diselesaikan
        if p == n:
            return "SPL tidak memiliki solusi unik."
        # Namun, kalau bisa diselesaikan, lanjut...
        else:
            # Tadinya, p melihat baris ke-i.
            # Kalau p sudah pindah ke bawahnya (sudah tidak sama dengan i),
            # artinya elemen diagonal saat ini bernilai nol, dan perlu ditukar 
            # dengan baris di bawahnya yang nilainya taknol (yaitu yang sedang
            # ditunjuk oleh indeks p). Maka tukarlah
            if p != i:
                matriks[[p,i], :] = matriks[[i,p], :]
                # syntax khusus numpy untuk menukar baris ke-i dan
                # baris ke-p, di mana semua nilai per kolom masih sama,
                # maksudnya tidak ada kolom yang ditukar, sehingga ditulis :
            
            # Ada pertukaran maupun tidak, yang pasti, sekarang elemen diagonal
            # sudah aman, sudah pasti taknol. Mari lanjut ke proses eliminasi.
            # Lakukan untuk tiap baris ke-j, yaitu untuk semua baris di bawah
            # elemen diagonal.
            for j in range (i+1, n):
                # Melakukan proses eliminasi dengan OBE (sesuai rumus di atas)
                m = matriks[j,i]/matriks[i,i] # m: "multiplier" atau pengkali
                matriks[j] = matriks[j] - m * matriks[i]
                #   (E_j) <- (   E_j    - m *    E_i   )
    
    # Setelah semua itu dilakukan untuk tiap kolom, eliminasi Gauss selesai
    return matriks
def SubstitusiBalik(matriks_input):
    # jaga-jaga
    matriks = matriks_input.astype(float)

    # memperoleh ukuran baris dari matriks diperbesar
    n = np.shape(matriks)[0]

    # vektor solusi, sementara isi dengan nol dulu
    solution = np.zeros(n)

    # lakukan dulu yang paling mudah, yaitu untuk baris paling bawah
    solution[n-1] = matriks[n-1, n]/matriks[n-1, n-1]

    # untuk baris-baris di atasnya, kita lakukan for loop, menghitung mundur,
    # terapkan rumus substitusi balik
    for i in range (n-2, -1, -1):
        # hitung sumasi, simpan langsung ke matriks[i,n]
        # agar langsung dijumlahkan ke b_i yaitu a_{i,n+1}
        for j in range(i+1, n):
            matriks[i,n] = matriks[i,n] - matriks[i,j] * solution[j]
        # peroleh solusi menggunakan rumus (dan memanfaatkan hasil sumasi)
        solution[i] = matriks[i,n]/matriks[i,i]
    return solution
aug_matriks = np.array(eval(input('Masukkan matriks diperbesar dari SPL yang akan diselesaikan: ')))
# mengubah input Anda ke dalam array numpy (matriks)
print("Berikut matriks yang dimasukkan:")
print(aug_matriks)

triangular_form = EliminasiGauss(aug_matriks)
if type(triangular_form) == type(""):
    # kalau output berupa string, artinya SPL tidak bisa diselesaikan
    print(triangular_form)
else:
    # Namun, kalau eliminasi Gauss berhasil, lanjut ke substitusi balik
    solution = SubstitusiBalik(triangular_form)
    print('Solusi dari SPL tersebut adalah: ')
    for i in range(len(solution)):
        print('x{0} = {1}'.format(i+1, solution[i]))
Masukkan matriks diperbesar dari SPL yang akan diselesaikan: [[1,-1,2,-1, -8],[2,-2,3,-3,-20],[1,1,1,0,-2], [1,-1,4,3,4]]
Berikut matriks yang dimasukkan:
[[  1  -1   2  -1  -8]
 [  2  -2   3  -3 -20]
 [  1   1   1   0  -2]
 [  1  -1   4   3   4]]
Solusi dari SPL tersebut adalah: 
x1 = -7.0
x2 = 3.0
x3 = 2.0
x4 = 2.0

Contoh penggunaan langsung (tanpa perlu menerima input):

matriks_diperbesar = np.array([
    [0.003, 59.14, 59.17],
    [5.291, -6.13, 46.78]
])

# langsung print vektor solusi
print(SubstitusiBalik(EliminasiGauss(matriks_diperbesar)))
[10.  1.]
matriks_diperbesar = np.array([
    [4.0, -1, 0, -1, 0, 0, 0, 0, 0, 25],
    [-1, 4, -1, 0, -1, 0, 0, 0, 0, 50],
    [0, -1, 4, 0, 0, -1, 0, 0, 0, 150],
    [-1, 0, 0, 4, -1, 0, -1, 0, 0, 0],
    [0, -1, 0, -1, 4, -1, 0, -1, 0, 0],
    [0, 0, -1, 0, -1, 4, 0, 0, -1, 50],
    [0, 0, 0, -1, 0, 0, 4, -1, 0, 0],
    [0, 0, 0, 0, -1, 0, -1, 4, -1, 0],
    [0, 0, 0, 0, 0, -1, 0, -1, 4, 25]
])
print(matriks_diperbesar)

# langsung print vektor solusi
print(SubstitusiBalik(EliminasiGauss(matriks_diperbesar)))
[[  4.  -1.   0.  -1.   0.   0.   0.   0.   0.  25.]
 [ -1.   4.  -1.   0.  -1.   0.   0.   0.   0.  50.]
 [  0.  -1.   4.   0.   0.  -1.   0.   0.   0. 150.]
 [ -1.   0.   0.   4.  -1.   0.  -1.   0.   0.   0.]
 [  0.  -1.   0.  -1.   4.  -1.   0.  -1.   0.   0.]
 [  0.   0.  -1.   0.  -1.   4.   0.   0.  -1.  50.]
 [  0.   0.   0.  -1.   0.   0.   4.  -1.   0.   0.]
 [  0.   0.   0.   0.  -1.   0.  -1.   4.  -1.   0.]
 [  0.   0.   0.   0.   0.  -1.   0.  -1.   4.  25.]]
[18.75 37.5  56.25 12.5  25.   37.5   6.25 12.5  18.75]

5. Partial Pivoting

Pada eliminasi Gauss, pertukaran baris perlu dilakukan ketika elemen diagonal bernilai nol. Namun, kalaupun tidak ada yang tepat bernilai nol, masalah round-off error bisa saja menyebabkan hasil komputasi meleset jauh ketika elemen diagonal bernilai sangat kecil mendekati nol (yang kemudian digunakan sebagai pembagi dalam eliminasi Gauss). Penyebabnya, pembagian oleh bilangan yang sangat kecil bisa sangat sensitif. Contohnya, \(\frac{10}{0.00025} = 40000\), sedangkan \(\frac{10}{0.00020} = 50000\). Dengan begitu, apabila nilai diagonal yang sangat kecil itu salah sedikit (karena masalah floating-point precision atau sejenisnya), hasil akhir nantinya bisa menjadi sangat meleset.

Solusi yang paling sederhana adalah memodifikasi elminasi Gauss, yaitu agar selalu menukarkan elemen diagonal \(a_{ii}\) dengan elemen terbesar di kolom ke-\(i\) yang ada di bawahnya (tentu saja, keseluruhan baris ikut ditukar, bukan hanya dua nilai). Solusi ini disebut partial pivoting, dan apapun modifikasi pada eliminasi Gauss untuk menghindari masalah di atas disebut strategi pivoting.

Dalam menerapkan strategi pivoting, algoritma substitusi balik tetap sama persis, karena hanya algoritma eliminasi Gauss yang dimodifikasi.

argmax dan argmin

Sebelum membahas strategi pivoting, mari kita bahas argmax dan argmin. Kedua fungsi ini tersedia dari numpy. Perbedaannya dengan max dan min cukup sederhana: max dan min mengembalikan nilainya, sedangkan argmax dan argmin mengembalikan indeksnya.

arraykecil = np.array([17, 8, 27, 54, 34])

nilai_max = np.max(arraykecil)
nilai_argmax = np.argmax(arraykecil)

print(f"Nilai maksimum ada pada indeks {nilai_argmax}, yaitu {nilai_max}")

nilai_min = np.min(arraykecil)
nilai_argmin = np.argmin(arraykecil)

print(f"Nilai minimum ada pada indeks {nilai_argmin}, yaitu {nilai_min}")
Nilai maksimum ada pada indeks 3, yaitu 54
Nilai minimum ada pada indeks 1, yaitu 8

Dengan demikian, apapun konteksnya, apabila kita memerlukan indeks letaknya saja, kita bisa langsung menggunakan argmax atau argmin.

Fun fact: baik argmax maupun argmin dikenal dalam matematika, dalam pembahasan teoritis juga. Contohnya, jika \(f(x) = 2-x^2\),

\[\underset{x\in\mathbb{N}}{\arg\max} f(x) = 1\]

karena bilangan asli \(x\) yang membuat nilai \(f(x)\) paling besar (di antara semua pilihan bilangan asli lainnya) adalah \(x=1\).

Partial Pivoting

Untuk eliminasi Gauss, daripada memeriksa apakah elemen diagonal bernilai nol atau tidak, selalu pilihlah indeks \(p \ge k\) terkecil* sedemikian sehingga,

\[|a_{pk}^{(k)}| = \max_{k \le i \le n}|a_{ik}^{(k)}|\]

di mana \(k\) adalah indeks untuk baris/kolom dari elemen diagonal yang sedang diurus (\(a_{kk}\)).

*apabila ada lebih dari satu baris yang sama-sama memuat nilai terbesar, pilih saja yang pertama kali ditemukan

Intinya, kita mencari indeks \(p\) untuk baris yang memuat elemen maksimum (dari semua elemen di bawah elemen diagonal), agar baris tersebut bisa ditukar dengan baris ke-\(k\) yang memuat elemen diagonal yang sedang diurus.

(“Perpangkatan” dengan \((k)\) itu sebenarnya hanya menandakan bahwa, pada saat itu, kita sedang mengurus elemen diagonal pada baris/kolom ke-\(k\). Penulisan seperti itu cukup untuk pembahasan teoritis saja. Tujuannya hanya untuk menekankan bahwa, setelah tiap OBE, nilai koefisien bisa jadi berbeda, sehingga harus diberi label tambahan seperti itu untuk memperjelas, nilai koefisien pada tahapan mana yang dimaksud.)

Setelah indeks \(p\) diperoleh, barulah lakukan pertukaran baris \((E_k) \leftrightarrow (E_p)\)

Implementasi Partial Pivoting

def PartialPivoting(matriks_input):
    # konversi matriks_input jadi matriks baru yang isinya float semua,
    # karena apabila ada bilangan bulat, bisa jadi dilakukan integer division
    # yang bisa sangat memperparah error
    matriks = matriks_input.astype(float)

    # n adalah banyaknya baris dari matriks diperbesar
    n = np.shape(matriks)[0]

    # untuk tiap kolom ke-i kecuali dua kolom terakhir
    for i in range(n-1):
        # Kumpulkan semua nilai yang ada di bawah elemen diagonal (pivot)
        below_pivot = abs(matriks[i:,i])
        # yaitu semua elemen pada baris di bawah elemen diagonal,
        # tetapi pada kolom yang sama.
        # Dibuat nilai mutlak karena yang diperhatikan hanya besarnya,
        # apakah dekat dengan nol atau tidak, bukan positif/negatifnya

        # Memilih indeks baris yang memuat elemen maksimum, sebagai "pivot" baru
        pivot_row = np.argmax(below_pivot)
        # Nilai yang disimpan itu sebenarnya adalah indeks "pergeseran" ke bawah
        # Misalnya, apabila variabel pivot_row bernilai 2, artinya baris yang
        # dipilih ada pada indeks (i+2), atau secara umum ditulis pivot_row+i

        # jika nilai pada baris yang akan di-pivot itu juga nol (padahal sudah
        # maksimum), maka sebenarnya semua nilai yang bisa ditukar itu nol semua
        # sehingga SPL tidak mungkin bisa diselesaikan
        if matriks[i,pivot_row+i] == 0:
            return "Tidak ada solusi unik"
        else:
            # Apabila taknol, lakukan pertukaran baris
            matriks[[pivot_row+i,i], :]= matriks[[i,pivot_row+i], :]

        # melanjutkan eliminasi Gauss seperti biasa
        for j in range(i+1,n):
            m = matriks[j,i]/matriks[i,i]
            matriks[j] = matriks[j]-m*matriks[i]
    
    return matriks
matriks = np.array(eval(input('Masukkan matriks yang akan dipivotkan: ')))

triangular_form = PartialPivoting(matriks)

print("Triangular matriksnya adalah :\n {0}".format(triangular_form))
Masukkan matriks yang akan dipivotkan: [[51, -18, 21, -96, -93], [84, -69, 69, 67, -6], [-42, 50, 14, -80, 51], [2, 8, 7, 3, 6]]
Triangular matriksnya adalah :
 [[ 8.40000000e+01 -6.90000000e+01  6.90000000e+01  6.70000000e+01
  -6.00000000e+00]
 [ 7.10542736e-15  2.38928571e+01 -2.08928571e+01 -1.36678571e+02
  -8.93571429e+01]
 [-4.60949996e-15  0.00000000e+00  6.20538117e+01  4.21674141e+01
   1.05968610e+02]
 [-1.84336495e-15  0.00000000e+00  0.00000000e+00  4.71963193e+01
   1.86585489e+01]]

Contoh penggunaan langsung (tanpa perlu menerima input):

matriks_diperbesar = np.array([
    [51, -18, 21, -96, -93],
    [84, -69, 69, 67, -6],
    [-42, 50, 14, -80, 51],
    [2, 8, 7, 3, 6]
])
print(matriks_diperbesar)

matriks_triangular = PartialPivoting(matriks_diperbesar)
print(matriks_triangular)

# vektor solusi
print(SubstitusiBalik(matriks_triangular))
[[ 51 -18  21 -96 -93]
 [ 84 -69  69  67  -6]
 [-42  50  14 -80  51]
 [  2   8   7   3   6]]
[[ 8.40000000e+01 -6.90000000e+01  6.90000000e+01  6.70000000e+01
  -6.00000000e+00]
 [ 7.10542736e-15  2.38928571e+01 -2.08928571e+01 -1.36678571e+02
  -8.93571429e+01]
 [-4.60949996e-15  0.00000000e+00  6.20538117e+01  4.21674141e+01
   1.05968610e+02]
 [-1.84336495e-15  0.00000000e+00  0.00000000e+00  4.71963193e+01
   1.86585489e+01]]
[-1.74956515 -0.22002462  1.4390443   0.39533907]
matriks_diperbesar = np.array([
    [4.0, -1, 0, -1, 0, 0, 0, 0, 0, 25],
    [-1, 4, -1, 0, -1, 0, 0, 0, 0, 50],
    [0, -1, 4, 0, 0, -1, 0, 0, 0, 150],
    [-1, 0, 0, 4, -1, 0, -1, 0, 0, 0],
    [0, -1, 0, -1, 4, -1, 0, -1, 0, 0],
    [0, 0, -1, 0, -1, 4, 0, 0, -1, 50],
    [0, 0, 0, -1, 0, 0, 4, -1, 0, 0],
    [0, 0, 0, 0, -1, 0, -1, 4, -1, 0],
    [0, 0, 0, 0, 0, -1, 0, -1, 4, 25]
])
print(matriks_diperbesar)

# langsung print vektor solusi
print(SubstitusiBalik(PartialPivoting(matriks_diperbesar)))
[[  4.  -1.   0.  -1.   0.   0.   0.   0.   0.  25.]
 [ -1.   4.  -1.   0.  -1.   0.   0.   0.   0.  50.]
 [  0.  -1.   4.   0.   0.  -1.   0.   0.   0. 150.]
 [ -1.   0.   0.   4.  -1.   0.  -1.   0.   0.   0.]
 [  0.  -1.   0.  -1.   4.  -1.   0.  -1.   0.   0.]
 [  0.   0.  -1.   0.  -1.   4.   0.   0.  -1.  50.]
 [  0.   0.   0.  -1.   0.   0.   4.  -1.   0.   0.]
 [  0.   0.   0.   0.  -1.   0.  -1.   4.  -1.   0.]
 [  0.   0.   0.   0.   0.  -1.   0.  -1.   4.  25.]]
[18.75 37.5  56.25 12.5  25.   37.5   6.25 12.5  18.75]

6. Scaled Partial Pivoting

Sebelumnya, kita hanya mencari baris dengan nilai terbesar yang berada di bawah elemen diagonal, untuk menghindari kemungkinan elemen diagonal terlalu kecil. Namun, apabila elemen diagonal setelah pertukaran itu menjadi sangat besar, sama saja semua elemen lainnya menjadi relatif sangat kecil, sehingga bisa timbul masalah yang sama.

Oleh karena itu, ada baiknya kita memodifikasi (lagi) syarat pemilihan baris untuk ditukar dengan elemen diagonal, yaitu memilih semacam “pertengahan”, daripada sekedar memilih yang paling besar. Pada scaled partial pivoting, kita

Scaled Partial Pivoting

Definisikan

\[s_i = \max_{k \le i \le n} |a_{ij}|\]

Pilih \(p \le k\) terkecil sedemikian sehingga

\[\frac{|a_{pk}^{(k)}|}{s_k} = \max_{k \le i \le n} \frac{a_{ik}^{(k)}}{s_i}\]

Kemudian, lakukan operasi \(\left( E_k \right) \leftrightarrow \left( E_p \right)\)

Implementasi Scaled Partial Pivoting

def ScaledPartialPivoting(matriks_input):
    # konversi matriks_input jadi matriks baru yang isinya float semua,
    # karena apabila ada bilangan bulat, bisa jadi dilakukan integer division
    # yang bisa sangat memperparah error
    matriks = matriks_input.astype(float)

    # memperoleh banyaknya baris pada matriks
    n = np.shape(matriks)[0]

    # menentukan scalar tiap kolom dibandingan masing-masing baris yang paling besar
    s = np.array([max(abs(matriks[i,:n])) for i in range(n)])

    # Apabila ada scalar yang nol, semua nilai pada baris tersebut nol,
    # sehingga SPL tidak bisa diselesaikan
    if 0 in s:
        return "tidak ada solusi unik"
    # Kalau bisa diselesaikan, lanjut...
    for i in range(n-1):
        below_pivot = abs(matriks[i:,i])/s[i:]
        pivot_row = np.argmax(below_pivot)
        if matriks[i,pivot_row+i] == 0:
            return "Tidak ada solusi unik"
        else:
            matriks[[pivot_row+i,i], :] = matriks[[i,pivot_row+i],  :]
            s[pivot_row+i],s[i]=s[i],s[pivot_row+i]

        # lanjut eleminasi Gauss
        for j in range(i+1,n):
            m = matriks[j,i]/matriks[i,i]
            matriks[j] = matriks[j]-m*matriks[i]
    return matriks
matriks = np.array(eval(input('Masukkan matriks yang akan dipivotkan: ')))

triangular_form = ScaledPartialPivoting(matriks)

print("Triangular matriksnya adalah :\n {0}".format(triangular_form))
Masukkan matriks yang akan dipivotkan: [[51, -18, 21, -96, -93], [84, -69, 69, 67, -6], [-42, 50, 14, -80, 51], [2, 8, 7, 3, 6]]
Triangular matriksnya adalah :
 [[ 8.40000000e+01 -6.90000000e+01  6.90000000e+01  6.70000000e+01
  -6.00000000e+00]
 [ 0.00000000e+00  9.64285714e+00  5.35714286e+00  1.40476190e+00
   6.14285714e+00]
 [ 0.00000000e+00  1.77635684e-15  3.98888889e+01 -4.87580247e+01
   3.81259259e+01]
 [ 7.10542736e-15  1.52153128e-15  0.00000000e+00 -1.81922748e+02
  -7.19211699e+01]]

Contoh penggunaan langsung (tanpa perlu menerima input):

matriks_diperbesar = np.array([
    [51, -18, 21, -96, -93],
    [84, -69, 69, 67, -6],
    [-42, 50, 14, -80, 51],
    [2, 8, 7, 3, 6]
])
print(matriks_diperbesar)

matriks_triangular = ScaledPartialPivoting(matriks_diperbesar)
print(matriks_triangular)

# vektor solusi
print(SubstitusiBalik(matriks_triangular))
[[ 51 -18  21 -96 -93]
 [ 84 -69  69  67  -6]
 [-42  50  14 -80  51]
 [  2   8   7   3   6]]
[[ 8.40000000e+01 -6.90000000e+01  6.90000000e+01  6.70000000e+01
  -6.00000000e+00]
 [ 0.00000000e+00  9.64285714e+00  5.35714286e+00  1.40476190e+00
   6.14285714e+00]
 [ 0.00000000e+00  1.77635684e-15  3.98888889e+01 -4.87580247e+01
   3.81259259e+01]
 [ 7.10542736e-15  1.52153128e-15  0.00000000e+00 -1.81922748e+02
  -7.19211699e+01]]
[-1.74956515 -0.22002462  1.4390443   0.39533907]
matriks_diperbesar = np.array([
    [4.0, -1, 0, -1, 0, 0, 0, 0, 0, 25],
    [-1, 4, -1, 0, -1, 0, 0, 0, 0, 50],
    [0, -1, 4, 0, 0, -1, 0, 0, 0, 150],
    [-1, 0, 0, 4, -1, 0, -1, 0, 0, 0],
    [0, -1, 0, -1, 4, -1, 0, -1, 0, 0],
    [0, 0, -1, 0, -1, 4, 0, 0, -1, 50],
    [0, 0, 0, -1, 0, 0, 4, -1, 0, 0],
    [0, 0, 0, 0, -1, 0, -1, 4, -1, 0],
    [0, 0, 0, 0, 0, -1, 0, -1, 4, 25]
])
print(matriks_diperbesar)

# langsung print vektor solusi
print(SubstitusiBalik(ScaledPartialPivoting(matriks_diperbesar)))
[[  4.  -1.   0.  -1.   0.   0.   0.   0.   0.  25.]
 [ -1.   4.  -1.   0.  -1.   0.   0.   0.   0.  50.]
 [  0.  -1.   4.   0.   0.  -1.   0.   0.   0. 150.]
 [ -1.   0.   0.   4.  -1.   0.  -1.   0.   0.   0.]
 [  0.  -1.   0.  -1.   4.  -1.   0.  -1.   0.   0.]
 [  0.   0.  -1.   0.  -1.   4.   0.   0.  -1.  50.]
 [  0.   0.   0.  -1.   0.   0.   4.  -1.   0.   0.]
 [  0.   0.   0.   0.  -1.   0.  -1.   4.  -1.   0.]
 [  0.   0.   0.   0.   0.  -1.   0.  -1.   4.  25.]]
[18.75 37.5  56.25 12.5  25.   37.5   6.25 12.5  18.75]

7. (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.]]