Tridiagonal

Metode eliminasi Gauss adalah suatu metode yang mudah digunakan, terutama ketika semua koefisien yang bukan nol terletak pada diagonal utama dan beberapa diagonal sekitarnya. Sistem yang memiliki karakteristik ini disebut sebagai “banded”, dan jumlah diagonal yang mengandung koefisien bukan nol disebut sebagai “bandwidth”. Salah satu contoh khusus yang sering ditemui adalah matriks tridiagonal yang memiliki bandwidth tiga. Proses eliminasi pada matriks tridiagonal dapat dianggap sebagai proses yang sederhana karena dengan menambahkan sebuah subdiagonal tambahan, proses substitusi mundur dapat dilakukan langsung. Berikut ini adalah fungsi penyelesaian untuk matriks tridiagonal:

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)
}

Dibawah ini adalah contoh soal menyelesaikan sistem persamaan dengan tridiagmatrix() Untuk menyelesaikannya, langkah pertama yang harus dilakukan adalah kita harus merubah persamaan tersebut kedalam bentuk matriks dengan membuat 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)

Masukkan vektor ke dalam fungsi tridiagmatrix().

tridiagmatrix(L=l, D=d, U=u, b=b)

## [1] 4 2 1 3

Terlebih dahulu buatlah augmented matriks dan gunakan sintaks gauss_jordan() untuk langkah terakhir.

scale_row <- function(m, row, k){
 m[row, ] <- m[row, ]*k
 return(m)
}

replace_row <- function(m, row1, row2, k){
  m[row2, ] <- m[row2, ] + m[row1, ]*k
  return(m)
}

gauss_jordan <- function (a){
    m <- nrow (a)
    n <- ncol (a)
    piv <- 1
    
    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)
                }
            }

            if(i != row_curr)
                a <- swap_row(a, i, row_curr)
            
            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)
}

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