Nama : Renata Amalia Putri
NIM : 220605110074
Kelas : C (Reguler)
Mata Kuliah : Kalkulus
Dosen Pengampu : Prof. Dr. Suhartono, M.kom
Jurusan : Teknik Informatika
Universitas : UIN Maulana Malik Ibrahim Malang
1. Operasi Vektor
Apabila menambahkan atau mengurangkan nilai elemen vektor dengan suatu skalar (konstanta yang hanya memiliki besaran), maka operasi penjumlahan atau pengurangan akan dilakukan pada setiap elemen vektor. Misal, jika melakukan penjumlahan pada vektor u dan v, maka operasi akan terjadi pada masing-masing elemen dnegan indeks yang sama.
Berikut contoh penerapannya pada bahasa pemrograman R :
u <- seq(4,6)
v <- seq(8,10)
# penjumlahan
u+v
## [1] 12 14 16
# pengurangan
u-v
## [1] -4 -4 -4
Contoh sintaks vektor pada R apabila salah satu vektor memiliki panjang yang berbeda :
x <- seq(2,5)
u+x
## Warning in u + x: longer object length is not a multiple of shorter object
## length
## [1] 6 8 10 9
Berdasarkan pada contoh tersebut, R akan mengeluarkan output peringatan/warning yang menunjukkan operasi dilakukan pada vektor dengan panjang berbeda. Tetapi R akan tetap melakukan perhitungan dengan menjumlahkan kembali vektor u yang belum dijumlahkan dengan vektor x hingga seluruh elemen vektor u dilakukan operasi penjumlahan.
Operasi lain yang dapat dilakukan pada vektor adalah menghitung inner product dan panjang vektor. Berikut adalah contoh bagaimana cara menghitung inner product dan panjang vektor menggunakan R :
# inner product
u%*%v
## [,1]
## [1,] 137
# panjang vektor u
sqrt(sum(u*u))
## [1] 8.774964
1.1 Operasi Matriks
Misalkan ada dua buah matriks A dan B, jika kedua matriks tersebut saling dijumlahkan maka penjumlahan dua buah matriks hanya dapat dilakukan pada matriks dengan ukuran yang seragam atau sama.
Berikut contoh operasi penjumlahan pada matriks :
A <- matrix(1:9,3)
B <- matrix(10:18,3)
C <- matrix(1:6,3)
# penjumlahan dengan skalar
A+1
## [,1] [,2] [,3]
## [1,] 2 5 8
## [2,] 3 6 9
## [3,] 4 7 10
# penjumlahan A+B
A+B
## [,1] [,2] [,3]
## [1,] 11 17 23
## [2,] 13 19 25
## [3,] 15 21 27
Operasi perhitungan lain pada matriks adalah operasi perkalian matriks. Untuk perkalian matriks, jumlah kolom matriks sebelah kiri harus sama dengan jumlah baris pada matriks sebelah kanan.
Pada bahasa R perkalian matriks dilakukan menggunakan operator **(%*%). Berikut adalah contoh penerapan perkalian matriks pada R**:
# Perkalian matriks
A%*%B
## [,1] [,2] [,3]
## [1,] 138 174 210
## [2,] 171 216 261
## [3,] 204 258 312
2. Operasi Baris Elementer
Ada tiga buah operasi dasar pada baris matriks operasi baris elementer. Ketiga operasi tersebut adalah :
Ketiga proses tersebut akan terjadi secara berulang, khususnya jika ingin mengerjakan sistem persamaan linier menggunakan algoritma eliminasi Gauss.
Fungsi row scalling pada R dapat dituliskan pada sintaks berikut :
scale_row <- function(m, row, k){
m[row, ] <- m[row, ]*k
return(m)
}
Berikut adalah contoh penerapannya :
# membuat matriks A
(A <- matrix(1:50, nrow=10))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 11 21 31 41
## [2,] 2 12 22 32 42
## [3,] 3 13 23 33 43
## [4,] 4 14 24 34 44
## [5,] 5 15 25 35 45
## [6,] 6 16 26 36 46
## [7,] 7 17 27 37 47
## [8,] 8 18 28 38 48
## [9,] 9 19 29 39 49
## [10,] 10 20 30 40 50
# lakukan scaling pada row 4 dengan nilai 12
scale_row(m=A, row=4, 12)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 11 21 31 41
## [2,] 2 12 22 32 42
## [3,] 3 13 23 33 43
## [4,] 48 168 288 408 528
## [5,] 5 15 25 35 45
## [6,] 6 16 26 36 46
## [7,] 7 17 27 37 47
## [8,] 8 18 28 38 48
## [9,] 9 19 29 39 49
## [10,] 10 20 30 40 50
Row swapping merupakan proses yang berulang, perlu menyimpan terlebih dahulu baris matriks pertama ke dalam sebuah objek. Baris matriks pertama selanjutnya diganti dengan baris matriks kedua, sedangkan baris matriks kedua selanjutnya akan diganti dengan baris matriks pertama yang telah disimpan terlebih dahulu ke dalam sebuah objek. Fungsi row swapping pada R dapat dituliskan pada sintaks berikut :
swap_row <- function(m, row1, row2){
row_tmp <- m[row1, ]
m[row1, ] <- m[row2, ]
m[row2, ] <- row_tmp
return(m)
}
Berikut adalah contoh penerapan fungsi swap_row() :
# pertukarkan baris 3 dengan baris 9
swap_row(m=A, row1=3, row2=9)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 11 21 31 41
## [2,] 2 12 22 32 42
## [3,] 9 19 29 39 49
## [4,] 4 14 24 34 44
## [5,] 5 15 25 35 45
## [6,] 6 16 26 36 46
## [7,] 7 17 27 37 47
## [8,] 8 18 28 38 48
## [9,] 3 13 23 33 43
## [10,] 10 20 30 40 50
Pada proses row replacement, proses perhitungan dilakukan dengan melakukan penjumlahan suatu baris matriks dengan baris matriks lainnya dengan terlebih dahulu melakukan row scalling terhadap matriks lainnya. Berikut adalah fungsi replace_row() yang ditulis pada R :
replace_row <- function(m, row1, row2, k){
m[row2, ] <- m[row2, ] + m[row1, ]*k
return(m)
}
Berikut adalah contoh penerapan fungsi replace_row() :
replace_row(m=A, row1=2, row2=4, k=-6)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 11 21 31 41
## [2,] 2 12 22 32 42
## [3,] 3 13 23 33 43
## [4,] -8 -58 -108 -158 -208
## [5,] 5 15 25 35 45
## [6,] 6 16 26 36 46
## [7,] 7 17 27 37 47
## [8,] 8 18 28 38 48
## [9,] 9 19 29 39 49
## [10,] 10 20 30 40 50
3. Eliminasi Gauss
Eliminasi Gauss merupakan sebuah cara untuk mencari penyelesaian sistem persamaan linier. Ide dasar dari eliminasi Gauss adalah melakukan operasi matematika pada baris matriks dan melanjutkannya sampai hanya tersisa satu variabel saja.
3.1 Row Echelon Form
Sebuah matriks merupakan row echelon form jika matriks tersebut memenuhi beberapa kondisi :
Contoh soal 3.1 Selesaikan sistem persamaan berikut :
x1 + x2 + x3 = 6
x1 + 2x2 - x3 = 2
2x1 + x2 + 2x3 = 10
Jawaban : Augmented matrix sistem persamaan linier tersebut adalah sebagai berikut :
[1 1 1 6]
Operasi baris elementer selanjutnya dilakukan pada matriks tersebut. Pada langkah pertama, baris ke-2 dikurangkan dengan baris ke-1 (B1 - B2) dan baris ke-3 dikurangkan oleh dua kali baris ke-1 (B1 - 2B1).
[1 1 1 6
1 2 -1 2
2 1 2 10]
B2 - B1 => B3 - 2B1
[1 1 1 6 0 1 -2 -4 0 -1 0 -2]
Selanjutnya, hasil dari langkah pertama menjadi input dari langkah selanjutnya. Pada langkah selanjutnya operasi baris elementer kembali dilanjutkan. Baris ke-3 dikurangkan dengan baris ke-2 (B3 - B2).
[1 1 1 6
0 1 -2 -4
0 -1 0 -2]
B3 - B2 =>
[1 1 1 6 0 1 -2 -4 0 0 -2 -6]
Setelah diperoleh matriks row echelon form selanjutnya penyelesaian persamaan :
x3 = -6/-2 = 3
x2 = 1/1 (-4-(2)3) = 2
x1 = 1/1 (6-2-3) = 1
Menyusun fungsi pada R untuk menyelesaikan sistem persamaan linier menggunakan matriks row echelon form. Fungsi yang akan dibentuk hanya sampai pad algoritma ke-4. Proses substitusi akan dilakukan secara manual. Berikut adalah sintaks yang digunakan :
ref_matrix <- function(a){
m <- nrow(a)
n <- ncol(a)
piv <- 1
# cek elemen diagonal apakah bernilai nol
for(row_curr in 1:m){
if(piv <= n){
i <- row_curr
while(a[i, piv] == 0 && i < m){
i <- i+1
if(i > m){
i <- row_curr
piv <- piv+1
if(piv > n)
return(a)
}
}
# jika diagonal bernilai nol, lakukan row swapping
if(i != row_curr)
a <- swap_row(a, i, row_curr)
# proses triangulasi untuk membentuk matriks segitiga atas
for(j in row_curr:m)
if(j != row_curr){
c <- a[j, piv]/a[row_curr, piv]
a <- replace_row(a, row_curr, j, -c)
}
piv <- piv+1
}
}
return(a)
}
Dengan menggunakan fungsi ref_matrix(), dapat membentuk matriks Row echelon form pada Contoh soal 3.1 sebagai berikut :
am <- c(1,1,2,
1,2,1,
1,-1,2,
6,2,10)
(m <- matrix(am, nrow=3))
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 6
## [2,] 1 2 -1 2
## [3,] 2 1 2 10
ref_matrix(m)
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 6
## [2,] 0 1 -2 -4
## [3,] 0 0 -2 -6
Contoh soal 3.2 Dengan menggunakan fungsi ref_matrix(), buatlah matriks row echelon form dari sistem persamaan linier berikut :
2x1 + x2 - x3 = 1
3x1 + 2x2 - 2x3 = 1
x1 - 5x2 + 4x3 = 3
Jawaban : Augmented matrix dari sistem persamaan tersebut adalah sebagai berikut :
[2 1 -1 1]
Penyelesaian matriks tersebut adalah sebagai berikut :
(m <- matrix(c(2,3,1,
1,2,-5,
-1,-2,4,
1,1,3), nrow=3))
## [,1] [,2] [,3] [,4]
## [1,] 2 1 -1 1
## [2,] 3 2 -2 1
## [3,] 1 -5 4 3
ref_matrix(m)
## [,1] [,2] [,3] [,4]
## [1,] 2 1.0 -1.0 1.0
## [2,] 0 0.5 -0.5 -0.5
## [3,] 0 0.0 -1.0 -3.0
Proses lebih lanjut akan menghasilkan penyelesaian sebagai berikut :
x1 = 1
x2 = 2
x3 = 3
3.2 Eliminasi Gauss-Jordan
Metode eliminasi Gauss-Jordan membentuk matriks menjadi bentuk reduced row echelon form. Metode ini merupakan pengembangan metode eliminasi Gauss, dimana sebelah kiri augmented matrix diubah menjadi matriks diagonal.
Contoh soal 3.3 Selesaikan sistem persamaan berikut :
x1 + x2 = 3
2x1 + 4x2 = 8
Jawaban :
Augmented Matrix dari persamaan linier tersebut adalah sebagai berikut :
[1 1 3]
Operasi baris elementer selanjutnya dilakukan pada matriks tersebut.
[1 1 3
2 4 8]
B2 - 2B1 =>
[1 1 3 0 2 2]
[1 1 3
0 2 2]
B2/2 =>
[1 1 3 0 1 1]
[1 1 3
0 1 1]
B1 - B2 =>
[1 0 2 0 1 1]
Penyelesaian persamaan linier tersebut adalah sebagai berikut :
x1 = 2 dan x2 = 1
Berdasarkan algoritma metode eliminasi Gauss-Jordan maka dapat membangun sebuah fungsi menggunakan R. Fungsi tersebut adalah sebagai berikut :
gauss_jordan <- function (a){
m <- nrow (a)
n <- ncol (a)
piv <- 1
# cek elemen diagonal utama apakah bernilai nol
for(row_curr in 1:m){
if(piv <= n){
i <- row_curr
while(a[i, piv] == 0 && i < m){
i <- i + 1
if(i > m){
i <- row_curr
piv <- piv + 1
if(piv > n)
return (a)
}
}
# jika diagonal utama bernilai nol,lakukan row swapping
if(i != row_curr)
a <- swap_row(a, i, row_curr)
# proses pembentukan matriks reduced row echelon form
piv_val <- a[row_curr , piv]
a <- scale_row (a, row_curr , 1 / piv_val)
for(j in 1: m){
if(j != row_curr){
k <- a[j, piv]/a[row_curr, piv]
a <- replace_row (a, row_curr, j, -k)
}
}
piv <- piv + 1
}
}
return (a)
}
Dengan menggunakan fungsi gauss_jordan(), sistem persamaan linier pada Contoh soal 3.3 :
(m <- matrix(c(1,2,1,4,3,8), nrow=2))
## [,1] [,2] [,3]
## [1,] 1 1 3
## [2,] 2 4 8
gauss_jordan(m)
## [,1] [,2] [,3]
## [1,] 1 0 2
## [2,] 0 1 1
Contoh soal 3.4 Dengan menggunakan fungsi gauss_jordan(), carilah penyelesaian sistem persamaan linier Contoh soal 3.1 dan Contoh soal 3.2 :
Jawaban : Untuk Contoh soal 3.1
am <- c(1,1,2,
1,2,1,
1,-1,2,
6,2,10)
m <- matrix(am, nrow=3)
gauss_jordan(m)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 1
## [2,] 0 1 0 2
## [3,] 0 0 1 3
Untuk Contoh soal 3.2
m <- matrix(c(2,3,1,1,2,-5,
-1,-2,4,1,1,3),
nrow=3)
gauss_jordan(m)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 1
## [2,] 0 1 0 2
## [3,] 0 0 1 3
4. Matrik Tridiagonal
Metode eliminasi Gauss merupakan metode yang sederhana untuk digunakan khususnya ketika semua koefisien bukan nol berkumpul pada diagonal utama dan beberapa diagonal sekitarnya. Suatu sistem yang bersifat demikian disebut sebagai banded dan banyaknya diagonal yang memuat koefisien bukan nol disebut dengan bandwidth.
Proses eliminasi untuk matriks tridiagonal bersifat trivial karena dengan membentuk sebuah subdiagonal tambahan, proses substitusi mundur bisa dapat dilakukan.
Berdasarkan algoritma matrik tridiagonal maka dapat membangun sebuah fungsi pada R. Fungsi penyelesaian matriks tridiagonal sebagai berikut :
tridiagmatrix <- function (L, D, U, b){
n <- length (D)
L <- c(NA , L)
## forward sweep
U[1] <- U[1] / D[1]
b[1] <- b[1] / D[1]
for(i in 2:(n - 1)){
U[i] <- U[i] / (D[i] - L[i] * U[i - 1])
b[i] <- (b[i] - L[i] * b[i - 1]) /
(D[i] - L[i] * U[i - 1])
}
b[n] <- (b[n] - L[n] * b[n - 1])/(D[n] - L[n] * U[n - 1])
## backward sweep
x <- rep.int (0, n)
x[n] <- b[n]
for(i in (n - 1) :1)
x[i] <- b[i] - U[i] * x[i + 1]
return (x)
}
Contoh soal 4.1 Selesaikan sistem persamaan berikut menggunakan fungsi tridiagmatrix() dan fungsi gauss_jordan()!
3x1 + 4x2 = 20
4x1 + 5x2 - 2x3 = 28
2x2 + 5x3 - 3x4 = 18
3x3 + 5x4 = 18
Jawaban :
[3 4 0 0 4 5 2 0 0 2 5 3 0 0 5 3] x = [20 28 18 18]
Untuk menyelesaikan persamaan tersebut menggunakan fungsi tridiagmatrix() maka perlu membentuk vektor diagonal l, d, u, dan b.
l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5)
b <- c(20, 28, 18, 18)
Lalu, vektor tersebut langsung dimasukkan ke dalam fungsi tridiagmatrix().
tridiagmatrix(L=l, D=d, U=u, b=b)
## [1] 4 2 1 3
Untuk meneyelesaikannya menggunakan fungsi gauss_jordan() maka diperlukan membentuk augmented matrix-nya terlebih dahulu.
m <- matrix(c(3,4,0,0,4,5,2,0,
0,2,5,3,0,0,3,5,
20,28,18,18), nrow=4)
gauss_jordan(m)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 4
## [2,] 0 1 0 0 2
## [3,] 0 0 1 0 1
## [4,] 0 0 0 1 3
5. Penyelesaian Sistem Persamaan Linier Menggunakan Fungsi solve()
Dalam pemrograman R menyediakan dungsi bawaan solve() untuk menyelesaikan sistem persamaan linier. Format fungsi solve() adalah sebagai berikut :
**solve**(a, b)
Berikut adalah contoh penerapan fungsi solve() pada sistem persamaan linier yang telah dipaparkan pada Contoh soal 3.2 :
# memecah matriks m menjadi matriks koefisien dan vektor konstanta
a <- matrix(c(2,3,1,1,2,-5,-1,-2,4),nrow=3)
b <- c(1,1,3)
solve(a,b)
## [1] 1 2 3
Tetapi, jika hanya memasukkan matriks persegi, maka output yang akan dihasilkan adalah invers dari matriks yang dimasukkan.
solve(a)
## [,1] [,2] [,3]
## [1,] 2 -1 7.401487e-17
## [2,] 14 -9 -1.000000e+00
## [3,] 17 -11 -1.000000e+00
Apabila mengalikan invers dengan matriks semula, maka akan dihasilkan output berupa matriks identitas.
a%*%solve(a)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
6. Penyelesaian Sistem Persamaan Linier Menggunakan Fungsi ‘Solve.tridiag()’
Penyelesaian matriks tridiagonal bisa juga dengan cara menggunakan fungsi Solve.tridiag() dari paket limSolve.
Fungsi Solve.tridiag() memiliki format sebagai berikut :
**Solve.tridiag**(diam1, dia, diap1, B=rep(0, times=length(dia)))
Berikut adalah contoh penerapannya pada matriks Contoh soal 4.1 :
library(limSolve)
## Warning: package 'limSolve' was built under R version 4.2.2
l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5)
b <- c(20, 28, 18, 18)
Solve.tridiag(diam1=l, dia=d, diap1=u, B=b)
## [,1]
## [1,] 4
## [2,] 2
## [3,] 1
## [4,] 3
7. Dekomposisi LU
Pada metode dekomposisi LU, matriks A difaktorkan menjadi matriks L dan matriks U, dimana kedua matriks tersebut harus sama dengan ukuran matriks A atau bahwa hasil perkalian kedua matriks tersebut akan menghasilkan matriks A.
A = LU
Dekomposisi LU didasarkan pada operasi baris elementer. Pertama, perlu menemukan matriks segitiga atas yang sesuai dengan matriks A. Solusi untuk melakukan dekomposisi bisa jadi tak terhingga, tetapi solusi yang paling sederhana adalah mengubah matriks A menjadi matriks row echelon form. Kedua, L harus menjadi matriks segitiga bawah yang mereduksi ke-l dengan mengikuti operasi baris yang sama yang menghasilkan U. Dapat menggunakan algoritma Doolittle untuk menghasilkan L, dimana nilai tiap entri dalam matriks segitiga bawah merupakan pengali yang digunakan untuk menghilangkan entri yang sesuai untuk tiap proses row replacement.
Proses eliminasi Gauss unutk memperoleh matriks U kadang menghasilkan nol di kolom pivotnya. Kondisi tersebut mengharuskan untuk melakukan proses row swapping atau pertukaran baris unutk pivot bukan nol. Jika proses tersebut berhasil dilakukan bisa jadi matriks A mungkin setara dengan matriks LU, tetapi tidak sama dalam hal urutan nilai pada tiap barisnya. Agar kita dapat memperoleh hasil yang sama (matriks A sama dengan matriks LU), diperlukan matriks ketiga, P. Matriks ini merupakan matriks identitas dengan ukuran sama dengan matriks A. Jika pertukaran baris dilakukan selama proses pembentukan matriks U, maka pertukaran baris yang sama juga akan diimplemenntasikan pada matriks P. Oleh karena itu, dalam praktiknya matriks A = PLU dan perkalian dengan matriks P berfungsi untuk mengembalikan urutan baris.
Berikut adalah sintaks uang digunakan untuk menyusun algoritma faktorisasi LU menggunakan bahasa R :
lu_solve <- function(a, b=NULL){
m <- nrow(a)
n <- ncol(a)
piv <- 1
# membentuk matriks identitas P dan L
P <- L <- diag(n)
# cek elemen diagonal utama apakah bernilai nol
for(row_curr in 1:m){
if(piv <= n){
i <- row_curr
while(a[i, piv] == 0 && i < m){
i <- i + 1
if(i > m){
i <- row_curr
piv <- piv + 1
if(piv > n)
return(list(P = P, L = L, U = a))
}
}
# jika elemen diagonal utama bernilai nol,lakukan row swapping
if(i != row_curr){
a <- swap_row(a, i, row_curr)
P <- swap_row(P, i, row_curr)
}
# pembentukan matriks L dan U
for(j in row_curr:m)
if(j != row_curr){
k <- a[j, piv]/a[row_curr, piv]
# matriks U
a <- replace_row(a, row_curr, j, -k)
# pengisian elemen matriks L
L[j, piv] <- k
}
piv <- piv + 1
}
}
# penyelesaian persamaan linier
if(is.null(b)){
return(list(P = P, L = L, U = a))
}else{
# forward substitution
t <- forwardsolve(L, b)
# backward substitution
x <- backsolve(a, t)
return(list(P = P, L = L, U = a, result=x))
}
}
# membuat matriks a dan vektor b
a <- matrix(c(1,2,3,-1,1,1,-1,2,
0,-1,-1,3,3,1,2,-1),
nrow=4)
b <- c(4,1,-3,4)
# penyelesaian
decomp<-lu_solve(a,b)
decomp$L%*%decomp$U%*%decomp$P
## [,1] [,2] [,3] [,4]
## [1,] 1 1 0 3
## [2,] 2 1 -1 1
## [3,] 3 -1 -1 2
## [4,] -1 2 3 -1
Melakukan prose dekomposisi menggunakna fungsi lu_solve().
# membentuk matriks a
(A <- matrix(c(0, 1, 7, 1, 5, -1, -2, 9, -5), 3))
## [,1] [,2] [,3]
## [1,] 0 1 -2
## [2,] 1 5 9
## [3,] 7 -1 -5
# dekomposisi lu
decomp<-lu_solve(A)
Melakukan pengecekan apakah matriks hasil dekomposisi akan menghasilkan matriks A.
decomp$P %*% decomp$L %*% decomp$U
## [,1] [,2] [,3]
## [1,] 0 1 -2
## [2,] 1 5 9
## [3,] 7 -1 -5
Fungsi lu() pada Package Matrix dapat digunakan unutk melakukan dekomposisi LU. Untuk menggunakannya hanya perlu menginputkan matriks ke dalam fungsi tersebut. Berikut contoh penerapannya :
library(Matrix)
## Warning: package 'Matrix' was built under R version 4.2.2
# membuat matriks a
a <- Matrix::Matrix(round(rnorm(9),2), nrow=3)
# dekomposisi
lum <- Matrix::lu(a)
lum
## 'MatrixFactorization' of Formal class 'denseLU' [package "Matrix"] with 4 slots
## ..@ Dimnames:List of 2
## .. ..$ : NULL
## .. ..$ : NULL
## ..@ x : num [1:9] 1.37 0.54 0.299 -1.14 -1.884 ...
## ..@ perm : int [1:3] 3 2 3
## ..@ Dim : int [1:2] 3 3
Untuk menampilkan hasil dekomposisi, jalankan fungsi expand().
decomp <- Matrix::expand(lum)
decomp
## $L
## 3 x 3 Matrix of class "dtrMatrix" (unitriangular)
## [,1] [,2] [,3]
## [1,] 1.0000000 . .
## [2,] 0.5401460 1.0000000 .
## [3,] 0.2992701 -0.6056404 1.0000000
##
## $U
## 3 x 3 Matrix of class "dtrMatrix"
## [,1] [,2] [,3]
## [1,] 1.370000 -1.140000 -0.630000
## [2,] . -1.884234 -1.709708
## [3,] . . -1.046928
##
## $P
## 3 x 3 sparse Matrix of class "pMatrix"
##
## [1,] . . |
## [2,] . | .
## [3,] | . .
8. Dekomposisi QR
Dekomposisi QR merupakan dekomposisi yang menyelesaikan sistem persamaan linier. Dekomposisi ini juga berperan penting untuk menghitung koefisien regresi dan pengaplikasian algoritma Newton-Raphson.
Untuk memperoleh informasi terkait dekomposisi ini, maka perlu mengetik sintaks berikut pada R :
?qr
Berikut merupakan contoh penerapan fungsi qr() untuk menyelesaikan sistem persamaan linier :
# membuat matriks A dan B
set.seed(123)
A <- matrix((1:12)+rnorm(12), nrow=4)
b <- 2:5
# dekomposisi matriks A
qr(A)
## $qr
## [,1] [,2] [,3]
## [1,] -6.3777985 -12.1257372 -19.850120
## [2,] 0.2774974 -6.3105459 -7.939245
## [3,] 0.7147777 -0.6461294 2.351193
## [4,] 0.6382310 -0.5653624 0.276672
##
## $rank
## [1] 3
##
## $qraux
## [1] 1.068915 1.512720 1.960964
##
## $pivot
## [1] 1 2 3
##
## attr(,"class")
## [1] "qr"
# memperoleh penyelesaian SPL
qr.solve(A,b)
## [1] 0.3045952 -0.1111081 0.3236862
9. Dekomposisi Cholesky
Dekomposisi Cholesky memberikan faktorisasi matriks alternatif sehingga A = LLT, dimana LT merupakan transpose konjugat dari matriks L.
Seperti dekomposisi LU, dekomposisi Cholesky dapat digunakan untuk menyelesaikan sistem persamaan linier. Tetapi, kelebihannya adalah menemukan dekomposisi Cholesky jauh lebih cepat daripada dekomposisi LU. Namun, dekomposisi ini hanya terbatas pada matriks tertentu saja. Dekomposisi Cholesky hanya dapat digunakan pada matriks definit positif dan simetris. Matriks simetris adalah matriks yang nilai di atas dan di bawah diagonalnya simetris atau sama. Definit positif berarti bahwa setiap entri pivot (nilai elemen diagonal utama) selalu bernilai positif.
Beikut adalah sintaks yang digunakan untuk menyusun fungsi pada R untuk melakukan dekomposisi Cholesky :
cholesky_solve <- function(a, b=NULL){
m <- nrow(a)
# membentuk matriks L dengan elemen nol
L = diag(0,m)
# Perhitungan elemen matriks L
for(i in 1:m){
for(k in 1:i){
p_sum <- 0
for(j in 1:k)
p_sum <- p_sum + L[j,i]*L[j,k]
# Pehitungan elemen diagonal utama
if(i==k)
L[k,i]<-sqrt(a[i,i]-p_sum)
else
L[k,i]<-(a[k,i]-p_sum)/L[k,k]
}
}
# Perhitungan elemn matriks L*
tL <- t(L)
# penyelesaian persamaan linier
if(is.null(b)){
return(list(L = L, tL = tL, a = a))
}else{
# forward substitution
t <- forwardsolve(L, b)
# backward substitution
x <- backsolve(tL, t)
return(list(L = L, tL = tL, a = a, result=x))
}
}
Dekomposisi Cholesky menggunakan fungsi cholesky_solve(), ditampilkan dalam sintaks berikut :
a <- matrix(c(9,-3,6,-3,17,-10,6,-10,12),3)
# dekomposisi Cholesky
(decomp<-cholesky_solve(a))
## $L
## [,1] [,2] [,3]
## [1,] 3 -1 2
## [2,] 0 4 -2
## [3,] 0 0 2
##
## $tL
## [,1] [,2] [,3]
## [1,] 3 0 0
## [2,] -1 4 0
## [3,] 2 -2 2
##
## $a
## [,1] [,2] [,3]
## [1,] 9 -3 6
## [2,] -3 17 -10
## [3,] 6 -10 12
# mengecek hasil dekomposisi
decomp$tL %*% decomp$L
## [,1] [,2] [,3]
## [1,] 9 -3 6
## [2,] -3 17 -10
## [3,] 6 -10 12
Fungsi lain untuk melakukan dekomposisi Cholesky adalah menggunakan fungsi chol() pada Package Matrix. Pada fungsi tersebut, hanya perlu menginputkan objek matrik ke dalamnya. Berikut adalah contoh penerapan fungsi tersebut menggunakan matriks :
chol(a)
## [,1] [,2] [,3]
## [1,] 3 -1 2
## [2,] 0 4 -2
## [3,] 0 0 2
10. Dekomposisi Eigen
Fungsi eigen() pada Package base dapat digunakan untuk melakukan dekomposisi eigen. Berikut sintaks dari dekomposisi eigen :
?eigen
Contoh sintaks untuk melakukan dekomposisi eigen :
A <- matrix(c(2,-1,0,-1,2,-1,0,-1,2), nrow=3)
# dekomposisi matriks A
eigen(A)
## eigen() decomposition
## $values
## [1] 3.4142136 2.0000000 0.5857864
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.5000000 -7.071068e-01 0.5000000
## [2,] 0.7071068 1.099065e-15 0.7071068
## [3,] -0.5000000 7.071068e-01 0.5000000
11. Metode Iterasi
Metode iterasi dimulai dengan estimasi nilai akhir. Setelah menerapkan beberapa perlakuan pada nilai estimasi, hasil perlakuan selanjutnya menjadi nilai estimasi untuk iterasi berikutnya. Proses tersebut akan berlangsung secara terus-menerus hingga ambang batas dipenuhi. Nilai ambang batas dapat berupa jumlah iterasi maksimum atau selisih antara nilai estimasi baru dan estimasi semula lebih kecil dari suatu nilai toleransi yang ditetapkan.
Jumlah kuadrat merupakan metode yang sering digunakan untuk mengecek apakah selisih nilai estimasi baru terhadap estimasi lama lebih kecil dari nilai toleransi yang ditetapkan.
11.1 Iterasi Jacobi
Dengan menggunakan metode jacobi, pertama-tama mengamati bahwa terdapat matriks R dan D yang memiliki hubungan A = R + D.
Suatu persamaan linier yang hendak diselesaikan dengan menggunakan metode iterasi Jacobi harus memenuhi syarat nilai elemen diagonal utama matriks harus lebih dominan. Artinya adalah nilai absolut diagonal utama matriks harus lebih besar dari jumlah nilai absolut elemen matriks lainnya pada satu kolom.
Untuk mempermudah proses iterasi, maka akan menggunakan bantuan R untuk melakukan komputasi. Langkah pertama yang perlu dilakukan adalah menyiapkan matriks A, vektor b, dan vektor x (nilai taksiran awal).
(A <- matrix(c(5,2,1,2,7,3,3,4,8), 3))
## [,1] [,2] [,3]
## [1,] 5 2 3
## [2,] 2 7 4
## [3,] 1 3 8
(b <- c(40,39,55))
## [1] 40 39 55
(x <- rep(0,3))
## [1] 0 0 0
Langkah selanjutnya adalah memperoleh invers matriks D.
(Dinv <- diag(1/diag(A)))
## [,1] [,2] [,3]
## [1,] 0.2 0.0000000 0.000
## [2,] 0.0 0.1428571 0.000
## [3,] 0.0 0.0000000 0.125
Persiapan terakhir sebelum iterasi dilakukan adalah menyiapkan matriks R.
(R<-A-diag(diag(A)))
## [,1] [,2] [,3]
## [1,] 0 2 3
## [2,] 2 0 4
## [3,] 1 3 0
iterasi 1
(x1 <- Dinv %*% (b-R%*%x))
## [,1]
## [1,] 8.000000
## [2,] 5.571429
## [3,] 6.875000
iterasi 2
(x2 <- Dinv %*% (b-R%*%x1))
## [,1]
## [1,] 1.6464286
## [2,] -0.6428571
## [3,] 3.7857143
iterasi 3
(x3 <- Dinv %*% (b-R%*%x2))
## [,1]
## [1,] 5.985714
## [2,] 2.937755
## [3,] 6.910268
Selama proses iterasi, jumlah akar, jumlah kuadrat dihitung. Berikut contoh akar jumlah kuadrat pada iterasi ke-3 :
sqrt(sum(x3-x2)^2)
## [1] 11.04445
Berikut adalah sintaks yang digunakan untuk menyusun fungsi sebuah fungsi dalam melakukan iterasi Jacobi :
jacobi <- function(a, b, tol=1e-7, maxiter=100){
n <- length(b)
iter <- 0
Dinv <- diag(1/diag(a))
R <- a-diag(diag(a))
x <- rep(0,n)
x_new <- rep(tol, n)
while(sqrt(sum(x_new-x)^2)>tol){
if(iter>maxiter){
warning("iterasi maksimum tercapai")
break
}
x <- x_new
x_new <- Dinv %*% (b - R %*% x)
iter <- iter+1
}
return(list(X = x_new, iter=iter))
}
Berikut adalah penerapan fungsi jacobi() tersebut :
jacobi(A,b)
## $X
## [,1]
## [1,] 4
## [2,] 1
## [3,] 6
##
## $iter
## [1] 62
Matriks A (matriks koefisien) berdasarkan sistem persamaan linier tersebut telah memenuhi syarat dari algoritma Jacobi (nilai diagonal utama dominan dibanding nilai lainnya pada satu kolom). Penyelesaian sistem persamaan tersebut, sebagai berikut:
A <- matrix(c(27,6,1,6,15,1,-1,2,54), 3)
b <- c(85,72,110)
jacobi(A,b)
## $X
## [,1]
## [1,] 2.425476
## [2,] 3.573016
## [3,] 1.925954
##
## $iter
## [1] 17
Nilai vektor x sesungguhnya dapat diperoleh menggunakan fungsi solve().
solve(A,b)
## [1] 2.425476 3.573016 1.925954
Berdasarkan hasil perhitungan, vektor x hasil iterasi memiliki nilai identik dengan nilai penyelesaian yang sebenarnya.
Penggunaan fungsi jacobi() syarat utama matriks haruslah terpenuhi, seperti: nilai diagonal matriks A lebih besar dari nilai elemen lainnya pada satu kolom. Selain itu, nilai diagonal matriks D tidak boleh sama dengan nol agar inver matriks D dapat diperoleh. Jika syarat-syarat tersebut terpenuhi, maka metode Jacobi dapat diterapkan.
11.2 Iterasi Gauss-Seidel
Metode iterasi Gauss-Seidel melakukan dekomposisi pada matriks A menjadi matriks segitiga atas U dan matriks segitiga bawah L. Dekomposisi ini tidak sama dengan dekomposisi LU. Matriks U pada metode Gauss-Seidel merupakan elemen (entri) matriks A pada bagian atas diagonal utama, sedangkan matriks L merupakan elemen diagonal utama dan bagian bawah diagonal utama matriks A. Elemen selain yang disebutkan pad akedua matriks tersebut akan bernilai nol.
Syarat agar suatu sistem persamaan linier dapat diselesaikan menggunakan metode Gauss-Seidel adalah matriks harus memiliki nilai diagonal utama yang dominan. Maksudnya, nilai absolut diagonal utama lebih besar dari jumlah nilai absolut elemen lainnya dalam satu kolom. Jika syarat ini tidak terpenuhi maka metode ini tidak akan memperoleh penyelesaian yang konvergen.
# membentuk matriks U dan L dari matriks A
(L <- U <- A)
## [,1] [,2] [,3]
## [1,] 27 6 -1
## [2,] 6 15 2
## [3,] 1 1 54
# membentuk matriks L dari entri bagian bawah diagonal utama matriks A
L[upper.tri(A, diag=FALSE)]<-0
L
## [,1] [,2] [,3]
## [1,] 27 0 0
## [2,] 6 15 0
## [3,] 1 1 54
# membentuk matriks U dari entri bagian atas diagonal utama matriks A
U[lower.tri(A, diag=TRUE)]<-0
U
## [,1] [,2] [,3]
## [1,] 0 6 -1
## [2,] 0 0 2
## [3,] 0 0 0
Selanjutnya, lakukan invers terhadap matriks L menggunakan fungsi solve().
(Linv <- solve(L))
## [,1] [,2] [,3]
## [1,] 0.0370370370 0.000000000 0.00000000
## [2,] -0.0148148148 0.066666667 0.00000000
## [3,] -0.0004115226 -0.001234568 0.01851852
Tetapkan nilai estimasi awal dan nilai toleransi yang diinginkan.
# tebakan awal nilai x
(x <- rep(0, length(b)))
## [1] 0 0 0
iterasi 1
(x1 <- Linv %*% (b - U %*% x))
## [,1]
## [1,] 3.148148
## [2,] 3.540741
## [3,] 1.913169
# akar jumlah kuadrat
sqrt(sum(x1-x)^2)
## [1] 8.602058
iterasi 2
(x2 <- Linv %*% (b - U %*% x1))
## [,1]
## [1,] 2.432175
## [2,] 3.572041
## [3,] 1.925848
# akar jumlah kuadrat
sqrt(sum(x2-x1)^2)
## [1] 0.6719939
Berikut merupakan sintaks yang digunakan untuk menyusun fungsi sebuah fungsi dalam melakukan iterasi Gauss-Seidel :
gauss_seidel <- function(a, b, tol=1e-7, maxiter=100){
n <- length(b)
iter <- 0
L <- U <- a
L[upper.tri(a, diag=FALSE)] <- 0
U[lower.tri(a, diag=TRUE)] <- 0
Linv <- solve(L)
x <- rep(0,n)
x_new <- rep(tol, n)
while(sqrt(sum(x_new-x)^2)>tol){
if(iter>maxiter){
warning("iterasi maksimum tercapai")
break
}
x <- x_new
x_new <- Linv %*% (b - U %*% x)
iter <- iter+1
}
return(list(X = x_new, iter=iter))
}
Menggunakan fungsi gauss_seidel() yang ditampilkan pada sintaks berikut :
gauss_seidel(A,b)
## $X
## [,1]
## [1,] 2.425476
## [2,] 3.573016
## [3,] 1.925954
##
## $iter
## [1] 7
Metode Kuadrat Terkecil
Metode kuadrat terkecil merupakan salah satu aplikasi penerapan aljabar linier yang paling populer. Intuisi dibalik metode ini adalah bagaimana kita meminimalkan jarak antara sejumlah titik dengan garis regresi. Misalkan kita menggambarkan scatterplot antara dua buah variabel. Pola yang terbentuk dari plot tersebut adalah terjadi korelasi positif antara variabel pada sumbu x dan sumbu y. Garis regresi terbaik terjadi ketika jumlah kuadrat jarak antara titik observasi dan garis regresi yang terbentuk seminimal mungkin.
Mneggunakan dataset **trees* yang berisi data hasil pengukuran kayu dari pohon yang ditebang.
head(trees)
## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7
str(trees)
## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
Berikut adalah sintaks yang digunakan untuk membentuk matriks A :
# membentuk matriks A
pred <- cbind(intercept=1, Girth=trees$Girth, Height=trees$Height)
head(A)
## [,1] [,2] [,3]
## [1,] 27 6 -1
## [2,] 6 15 2
## [3,] 1 1 54
Langkah selanjutnya adalah membentuk matriks b. Berikut adalah sintaks yang digunakan :
resp<- trees$Volume
head(resp)
## [1] 10.3 10.3 10.2 16.4 18.8 19.7
A <- t(pred) %*% pred
b <- t(pred) %*% resp
Ab <- cbind(A,b)
(x <- gauss_jordan(Ab))
## intercept Girth Height
## intercept 1 0 0 -57.9876589
## Girth 0 1 0 4.7081605
## Height 0 0 1 0.3392512
Menggunakan fungsi lain untuk memperoleh nilai koefisien tersebut, seperti : lu_solve() dan solve(). Untuk fungsi jacobi() dan gauss_seidel maka harus memastikan syarat-syarat terkait metode tersebut. Berikut adalah contoh penyelesain menggunakan sintaks lainnya :
# metode LU
lu_solve(A,b)
## $P
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
##
## $L
## [,1] [,2] [,3]
## [1,] 1.00000 0.000000 0
## [2,] 13.24839 1.000000 0
## [3,] 76.00000 1.054369 1
##
## $U
## intercept Girth Height
## intercept 31 410.7000 2356.0000
## Girth 0 295.4374 311.5000
## Height 0 0.0000 889.5641
##
## $result
## [,1]
## [1,] -57.9876589
## [2,] 4.7081605
## [3,] 0.3392512
# fungsi solve()
solve(A,b)
## [,1]
## intercept -57.9876589
## Girth 4.7081605
## Height 0.3392512
Bahasa pemrograman R juga menyediakan fungsi untuk membentuk model regresi linier. Fungsi yang digunakan adalah lm(). Berikut sintaks yang digunakan untuk membentuk model linier menggunakan fungsi lm() :
lm(Volume~Girth+Height, data=trees)
##
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
##
## Coefficients:
## (Intercept) Girth Height
## -57.9877 4.7082 0.3393
Aliran Massa Dalam Reaktor
Menggunakan metode eliminasi Gauss-Jordan, dekomposisi LU, iterasi Jacobi, dan iterasi Gauss-Seidel.
Yang pertama dilakukan adalah membentuk matriksnya terlebih dahulu :
(A <- matrix(c(6,-3,0,0,-3,
0,3,-1,-1,-1,
-1,0,9,-8,0,
0,0,0,11,0,
0,0,0,-2,4),nrow=5))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 6 0 -1 0 0
## [2,] -3 3 0 0 0
## [3,] 0 -1 9 0 0
## [4,] 0 -1 -8 11 -2
## [5,] -3 -1 0 0 4
(b <- c(50,0,160,0,0))
## [1] 50 0 160 0 0
gauss_jordan(cbind(A,b))
## b
## [1,] 1 0 0 0 0 11.50943
## [2,] 0 1 0 0 0 11.50943
## [3,] 0 0 1 0 0 19.05660
## [4,] 0 0 0 1 0 16.99828
## [5,] 0 0 0 0 1 11.50943
lu_solve(A,b)
## $P
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1
##
## $L
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0 0.0000000 0.0000000 0 0
## [2,] -0.5 1.0000000 0.0000000 0 0
## [3,] 0.0 -0.3333333 1.0000000 0 0
## [4,] 0.0 -0.3333333 -0.9245283 1 0
## [5,] -0.5 -0.3333333 -0.0754717 0 1
##
## $U
## [,1] [,2] [,3] [,4] [,5]
## [1,] 6 0 -1.000000e+00 0 0
## [2,] 0 3 -5.000000e-01 0 0
## [3,] 0 0 8.833333e+00 0 0
## [4,] 0 0 0.000000e+00 11 -2
## [5,] 0 0 1.110223e-16 0 4
##
## $result
## [1] 11.50943 11.50943 19.05660 16.99828 11.50943
jacobi(A,b, maxiter=100)
## $X
## [,1]
## [1,] 11.50943
## [2,] 11.50943
## [3,] 19.05660
## [4,] 16.99828
## [5,] 11.50943
##
## $iter
## [1] 17
gauss_seidel(A,b, maxiter=200)
## $X
## [,1]
## [1,] 11.50943
## [2,] 11.50943
## [3,] 19.05660
## [4,] 16.99828
## [5,] 11.50943
##
## $iter
## [1] 7
DAFTAR PUSTAKA
Rosidi, M. (2019). Metode Numerik Menggunakan R Untuk Teknik Lingkungan.
Primartha, R. (2018). Belajar Machine Learning Teori dan Praktik. Penerbit Informatika : Bandung.
Sanjaya, M. (2015). Metode Numerik Berbasis Phython. Penerbit Gava Media: Yogyakarta.
Howard, J.P. (2017). Computational Methods for Numerical Analysis with R. CRC Press.