import numpy as np(Pertemuan 8) Sistem Persamaan Linear: Jacobi, Gauss-Seidel
Jacobi method, Gauss-Seidel method, SOR method
Offline di Departemen Matematika
Kembali ke Metode Numerik
Materi Pembahasan:
Beberapa jenis norm vektor
Masalah copy untuk array numpy
Metode Jacobi, algoritma praktis
Metode Gauss-Seidel, algoritma praktis
Metode Relaksasi (SOR), algoritma praktis
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. Beberapa jenis norm vektor
Ada beberapa jenis norm vektor, yaitu semacam ukuran “besar” (magnitude) atau “panjang” untuk vektor. Salah satu kegunaannya adalah membandingkan “ukuran” antara dua vektor, yang mana yang lebih besar/kecil. Tiga jenis norm yang terkenal adalah:
- Euclidean norm, sering disebut 2-norm atau \(L_2\)-norm, dan perhitungannya seperti menggunakan teorema Pythagoras. Penulisan: \(||\textbf{v}||_2\)
- Infinity norm (norm tak hingga), terkadang disebut \(L_{\infty}\)-norm, yaitu sama saja maksimum dari semua mutlak elemen vektor. Penulisan: \(||\textbf{v}||_{\infty}\)
- Manhattan distance atau Taxicab norm, sering disebut 1-norm atau \(L_1\)-norm, yaitu menjumlahkan mutlak tiap elemen vektor. Penulisan: \(||\textbf{v}||_1\)
Numpy bisa menghitung beberapa jenis norm, termasuk ketiga jenis norm di atas, menggunakan numpy.linalg.norm(vektor, jenis_norm), di mana vektor dibuat dengan numpy.array.
vektor_kecil = np.array([3,-4])
print(vektor_kecil)[ 3 -4]Contoh norm-2, dengan option ord=2:
np.linalg.norm(vektor_kecil, ord=2)5.0Output adalah 5, karena \(||\textbf{v}||_2=\sqrt{3^2+\left(-4\right)^2}=\sqrt{25}=5\).
Sebenarnya, “ord” tidak harus ditulis:
np.linalg.norm(vektor_kecil, 2)5.0Contoh norm-infinity, dengan option ord=numpy.inf:
np.linalg.norm(vektor_kecil, np.inf)4.0Output adalah 4, karena \(||\textbf{v}||_{\infty} = \max \left( |3|, |-4| \right) = \max \left( 3, 4 \right) = 4\).
Contoh norm-1:
np.linalg.norm(vektor_kecil, 1)7.0Output adalah 7, karena \(||\textbf{v}||_1 = |3| + |-4| = 3+4=7\).
Masing-masing jenis norm memiliki kelebihan dan kekurangan. (Untuk ke depannya, kita akan menggunakan norm-infinity, sesuai buku Burden). Apapun jenis norm yang Anda gunakan, untuk perhitungan apapun, pastikan Anda konsisten selalu menggunakan jenis norm yang sama dari awal sampai akhir perhitungan.
Untuk jenis norm lainnya, bisa baca lebih lanjut di dokumentasi numpy (pada keterangan option “ord”): https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html
2. Masalah copy untuk array numpy
Ada satu hal yang perlu dibahas sebelum melanjutkan ke pembahasan metode iteratif untuk SPL.
Salah satu kekurangan numpy (dan Python secara umum) adalah bahwa kita tidak bisa meng-copy suatu array (ataupun list) dengan assignment. Melakukan assignment seolah-olah hanya membuat “sinonim”, sehingga perubahan pada salah satu array/list juga akan mengubah array/list yang satunya (istilahnya shallow copy). Perhatikan,
import numpy as nparray1 = np.array([9, 8, 7, 6])
print(array1)[9 8 7 6]# Apakah cara copy seperti ini?
array2 = array1print(array2)[9 8 7 6]Seandainya array2 diubah…
array2[2] = 15
print(array2)[ 9 8 15 6]… maka array1 juga mengalami perubahan yang sama.
print(array1)[ 9 8 15 6]Cara copy yang tepat untuk array maupun list adalah menggunakan akhiran .copy() seperti berikut ini.
array3 = array1.copy()
print(array3)[ 9 8 15 6]Sehingga, perubahan pada salah satu tidak akan mempengaruhi yang satunya lagi. Artinya, copy telah dilakukan secara sempurna (disebut deep copy).
array3[3] = 20
print(array3)
print(array1)[ 9 8 15 20]
[ 9 8 15 6]Untuk ke depannya, kita akan sering menggunakan .copy().
Fun fact: sepertinya, permasalahan shallow copy ini berhubungan erat dengan cara dibuatnya library numpy yang sebenarnya tersambung dengan bahasa pemrograman C, yang juga memiliki keanehan yang serupa untuk array.
3. Metode Jacobi, algoritma praktis
Misalkan vektor \(\textbf{x}^{(k)} = \left( x_1^{(k)}, x_2^{(k)}, \dots, x_n^{(k)} \right)^t\) (ditranspos agar berupa vektor kolom) adalah hasil aproksimasi pada iterasi ke-k untuk solusi SPL \(n\)-variabel \(A\textbf{x}=\textbf{b}\).
Metode Jacobi memiliki formula sebagai berikut:
\[x_i^{(k)} = \frac{1}{a_{ii}} \left[ \sum_{j=1,j\ne i}^{n} \left(-a_{ij}x_j^{(k-1)} \right) + b_i \right],\hspace{0.5cm} i = 1, 2, \dots, n \]
Kriteria pemberhentian iterasi bisa berupa * error mutlak: \(||\textbf{x}^{(k)}-\textbf{x}^{(k-1)}|| < \epsilon\) * error relatif: \(\frac{||\textbf{x}^{(k)}-\textbf{x}^{(k-1)}||}{||\textbf{x}^{(k)}||} < \epsilon\)
Pada kode berikut ini, kita akan menggunakan error mutlak dan norm tak hingga.
import numpy as np
def Jacobi(matriks, tebakan_awal, tol):
# banyaknya baris pada matriks, atau sama saja banyaknya variabel
n = np.shape(matriks)[0]
# definisikan vektor x0 sebagai tebakan awal
x0 = tebakan_awal.copy()
# error sementara (karena error belum diketahui), agar masuk while loop
err = tol + 1
while err > tol: # selama kriteria pemberhentian belum terpenuhi,
# anggap x0 sebagai x^(k-1) atau hasil iterasi sebelumnya,
# kemudian nilai yang baru akan disimpan pada vektor x^(k) berikut:
x = x0.copy()
# metode Jacobi untuk tiap i
for i in range(n):
b = matriks[i, n]
for j in range(n):
if j != i:
b = b - matriks[i,j] * x0[j]
# perhatikan bahwa selalu digunakan x0,
# yaitu nilai-nilai x^(k-1), yaitu dari iterasi sebelumnya
x[i] = b/matriks[i,i]
# update error mutlak
err = np.linalg.norm(x-x0, np.inf)
# memasuki iterasi selanjutnya,
# vektor x0 sekarang adalah vektor x yang barusan dihitung
x0 = x
# jika keluar while loop maka metode selesai, x(k) adalah vektor hasil akhir
return xmatriks_diperbesar = np.array(eval(input("Masukkan matriks diperbesar: "))).astype(float)
tebakan_awal = np.array(eval(input("Masukkan tebakan awal: ")), dtype=float)
# catatan: .astype(float) dan dtype=float melakukan hal yang sama,
# dengan cara penggunaan yang sedikit berbeda:
# - numpy.array(...).astype(float)
# - numpy.array(..., dtype=float)
# tidak ada salahnya apabila menggunakan salah satu saja (lebih baik konsisten)
toleransi = eval(input("Masukkan toleransi: "))
hasil = Jacobi(matriks_diperbesar, tebakan_awal, toleransi)
print("Hasil metode Jacobi adalah:")
print(hasil)Masukkan matriks diperbesar: [ [10, -1, 2, 0, 6], [-1, 11, -1, 3, 25], [2, -1, 10, -1, -11], [0, 3, -1, 8, 15] ]
Masukkan tebakan awal: [0,0,0,0]
Masukkan toleransi: 10**-4
Hasil metode Jacobi adalah:
[ 0.99998973 2.00001582 -1.00001257 1.00001924]4. Metode Gauss-Seidel, algoritma praktis
Metode Gauss-Seidel adalah modifikasi/perkembangan dari metode Jacobi, di mana semua nilai \(x_i\) yang digunakan untuk perhitungan adalah selalu yang terbaru. Artinya, ketika menghitung \(x_i^{(k)}\), meskipun nilai-nilai \(x_{i+1}, \dots, x_n\) yang digunakan adalah dari iterasi sebelumnya, nilai-nilai \(x_1, x_2, \dots, x_{i-1}\) yang digunakan adalah yang baru saja dihitung, yaitu dari iterasi yang sama. Hal ini tidak seperti metode Jacobi yang selalu menggunakan nilai-nilai dari iterasi sebelumnya.
Oleh karena itu, untuk metode Gauss-Seidel, penulisan sumasi dipisah menjadi dua bagian, yaitu satu bagian untuk penggunaan nilai-nilai dari iterasi yang sama \((k)\), dan satu bagian untuk penggunaan nilai-nilai dari iterasi sebelumnya \((k-1)\).
\[x_i^{(k)} = \frac{1}{a_{ii}} \left[ -\sum_{j=1}^{i-1} \left( a_{ij}x_j^{(k)} \right) - \sum_{j=i+1}^{n} \left( a_{ij}x_j^{(k-1)} \right) + b_i \right],\hspace{0.5cm} i = 1, 2, \dots, n \]
Akibat selalu menggunakan nilai-nilai terbaru, metode Gauss-Seidel cenderung lebih cepat konvergen memenuhi toleransi daripada metode Jacobi.
Secara algoritma, perubahan ini pun sebenarnya sangat kecil. Antara metode Jacobi dan metode Gauss-Seidel, perbedaannya hanya di satu baris saja…
import numpy as np
def GaussSeidel(matriks, tebakan_awal, tol):
# banyaknya baris pada matriks, atau sama saja banyaknya variabel
n = np.shape(matriks)[0]
# definisikan vektor x0 sebagai tebakan awal
x0 = tebakan_awal.copy()
# error sementara (karena error belum diketahui), agar masuk while loop
err = tol + 1
while err > tol: # selama kriteria pemberhentian belum terpenuhi,
# anggap x0 sebagai x^(k-1) atau hasil iterasi sebelumnya,
# kemudian nilai yang baru akan disimpan pada vektor x^(k) berikut:
x = x0.copy()
# metode Gauss-Seidel untuk tiap i
for i in range(n):
b = matriks[i, n]
for j in range(n):
if j != i:
# perubahan dari metode Jacobi hanya di baris berikut
b = b - matriks[i,j] * x[j]
# perhatikan bahwa selalu digunakan x,
# yaitu nilai-nilai x^(k), yaitu nilai-nilai terbaru;
x[i] = b/matriks[i,i]
# update error mutlak
err = np.linalg.norm(x-x0, np.inf)
# memasuki iterasi selanjutnya,
# vektor x0 sekarang adalah vektor x yang barusan dihitung
x0 = x
# jika keluar while loop maka metode selesai, x(k) adalah vektor hasil akhir
return xmatriks_diperbesar = np.array(eval(input("Masukkan matriks diperbesar: "))).astype(float)
tebakan_awal = np.array(eval(input("Masukkan tebakan awal: ")), dtype=float)
toleransi = eval(input("Masukkan toleransi: "))
hasil = GaussSeidel(matriks_diperbesar, tebakan_awal, toleransi)
print("Hasil Gauss-Seidel adalah:")
print(hasil)Masukkan matriks diperbesar: [ [10, -1, 2, 0, 6], [-1, 11, -1, 3, 25], [2, -1, 10, -1, -11], [0, 3, -1, 8, 15] ]
Masukkan tebakan awal: [0,0,0,0]
Masukkan toleransi: 10**-4
Hasil Gauss-Seidel adalah:
[ 1.00000836 2.00000117 -1.00000275 0.99999922]5. Metode Relaksasi (SOR), algoritma praktis
Metode relaksasi (relaxation methods) adalah metode untuk mempercepat kekonvergenan dari solusi yang dihasilkan oleh metode iteratif untuk SPL (dalam hal ini, metode Gauss-Seidel). Berdasarkan besar faktor relaksasi \(\omega\), metode relaksasi terbagi menjadi dua jenis, yaitu * metode under relaxation jika \(0 < \omega < 1\) * metode over relaxation jika \(\omega > 1\)
Sesuai buku Burden, pembahasan kita akan fokus ke metode over relaxation. Metode over relaxation yang diterapkan terus-menerus untuk tiap iterasi Gauss-Seidel disebut metode Successive Over-Relaxation (SOR).
Untuk sembarang nilai omega, rumus metode relaksasi sebagai modifikasi Gauss-Seidel bisa dituliskan sebagai berikut:
\[x_i^{(k)} = \left(1-\omega\right)x_i^{(k-1)} + \frac{\omega}{a_{ii}} \left[ -\sum_{j=1}^{i-1} \left( a_{ij}x_j^{(k)} \right) - \sum_{j=i+1}^{n} \left( a_{ij}x_j^{(k-1)} \right) + b_i \right],\hspace{0.5cm} i = 1, 2, \dots, n \]
Catatan: jika \(\omega = 1\), diperoleh metode Gauss-Seidel yang telah dibahas sebelumnya (tanpa relaksasi).
Lagi-lagi, perbedaan kode antara metode Gauss-Seidel dan metode SOR hanya di satu baris saja…
import numpy as np
def SOR(matriks, tebakan_awal, tol, omega=1):
# banyaknya baris pada matriks, atau sama saja banyaknya variabel
n = np.shape(matriks)[0]
# definisikan vektor x0 sebagai tebakan awal
x0 = tebakan_awal.copy()
# error sementara (karena error belum diketahui), agar masuk while loop
err = tol + 1
while err > tol: # selama kriteria pemberhentian belum terpenuhi,
# anggap x0 sebagai x^(k-1) atau hasil iterasi sebelumnya,
# kemudian nilai yang baru akan disimpan pada vektor x^(k) berikut:
x = x0.copy()
# metode Gauss-Seidel untuk tiap i
for i in range(n):
b = matriks[i, n]
for j in range(n):
if j != i:
b = b - matriks[i,j] * x[j]
# perhatikan bahwa selalu digunakan x,
# yaitu nilai-nilai x^(k), yaitu nilai-nilai terbaru;
# bedanya dengan metode Gauss-Seidel hanya di baris berikut:
x[i] = (1-omega) * x0[i] + omega*b/matriks[i,i] # hasil SOR
# update error mutlak
err = np.linalg.norm(x-x0, np.inf)
# memasuki iterasi selanjutnya,
# vektor x0 sekarang adalah vektor x yang barusan dihitung
x0 = x
# jika keluar while loop maka metode selesai, x(k) adalah vektor hasil akhir
return xmatriks_diperbesar = np.array(eval(input("Masukkan matriks diperbesar: "))).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 = SOR(matriks_diperbesar, tebakan_awal, omega, toleransi)
print("Hasil SOR adalah:")
print(hasil)Masukkan matriks diperbesar: [ [4, 3, 0, 24], [3, 4, -1, 30], [0, -1, 4, -24] ]
Masukkan tebakan awal: [0,0,0]
Masukkan faktor relaksasi (omega): 1.25
Masukkan toleransi: 10**-4
Hasil SOR adalah:
[ 2.99998919 4.00000321 -4.9999937 ]