import numpy as npModul 6 Praktikum Metode Numerik: Metode Iteratif untuk SPL
Modul 6 Praktikum Metode Numerik: Metode Iteratif untuk SPL
Kembali ke Metode Numerik
Ini adalah pertemuan terakhir praktikum Metode Numerik tahun ini.
Outline:
- Beberapa jenis norm vektor
- Masalah copy untuk array numpy
- Metode Jacobi, algoritma praktis
- Metode Gauss-Seidel, algoritma praktis
- Metode Relaksasi (SOR), algoritma praktis
- (Pengayaan) Metode Jacobi, dalam bentuk matriks (teoritis)
- (Pengayaan) Metode Gauss-Seidel, dalam bentuk matriks (teoritis)
- (Pengayaan) Metode SOR, dalam bentuk matriks (teoritis)
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 ]6. (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 xkmatriks_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]7. (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 xkmatriks_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]8. (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 xkmatriks_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 ]