Eliminasi Gauss

# Membuat matriks A dan vektor b 
A = matrix(c(3, 1, -1,
              4, 7, -3,
              2, -2, 5), nrow = 3, byrow = TRUE)
b = c(5, 20, 10)

# Fungsi Eliminasi Gauss
gauss_elim = function(A,b) {
  n = nrow(A)
  A = cbind(A, b)
  for (i in 1:(n-1)){
    for (j in (i+1):n){
    factor = A[j, i]/ A[i, i]
    A[j, ] = A[j, ] - factor*A[i, ]
    }
  }
  return(A)
}

# Lakukan eliminasi 
A_gauss = gauss_elim(A,b)
A_gauss
##                                  b
## [1,] 3 1.000000 -1.000000  5.00000
## [2,] 0 5.666667 -1.666667 13.33333
## [3,] 0 0.000000  4.882353 12.94118

Terbentuk matriks A_gauss dengan segitiga atas, untuk solusinya dilakukan subtitusi balik

Subtitusi Balik

# Fungsi subtitusi balik 
subs_balik = function(A){
  n = nrow(A)
  x = numeric(n)
  x[n] =A[n, n+1]/A[n, n]
  
  for (i in (n-1):1){
    x[i] =(A[i, n+1] - sum(A[i, (i+1):n]*x[(i+1):n]))/A[i,i]
  }
  return(x)
}

# Solusi 
solusi_gauss = subs_balik(A_gauss)
solusi_gauss
## [1] 1.506024 3.132530 2.650602

Kelebihan:

  • Mudah dipahami dan dipraktikan untuk SPL kecil.

  • Sapat digunakan pada sistem persamaan dengan koefisien real.

Kekurangan:

  • Rentan terhadap error pembulatan pada matriks besar.

  • Membutuhkan operasi komputasi yang signifikan untuk matriks berukuran besar.

Eliminasi Gauss Jordan

# Eliminasi Gauss Jordan 
gauss_j_elim = function(A,b){
n = nrow(A)
  A = cbind(A, b)
  for (k in 1:n) {
    A[k, ] = A[k, ]/A[k, k]
    for (i in 1:n) {
      if (i != k) {
        factor = A[i, k]/A[k, k]
        A[i, ] = A[i, ] - factor * A[k, ]
      }
    }
  }
  return(A)
}

# Lakukan eliminasi 
A_gauss_jordan = gauss_j_elim(A, b)
A_gauss_jordan
##                   b
## [1,] 1 0 0 1.506024
## [2,] 0 1 0 3.132530
## [3,] 0 0 1 2.650602

sehingga bentuk x, y, dan z adalah

Solusi

A_gauss_jordan[, 4]
## [1] 1.506024 3.132530 2.650602

Kelebihan:

  • Memberikan solusi langsung tanpa subtitusi balik.

  • Berguna untuk analisis data atau aljabar linier simbolik.

Kekurangan:

  • Lebih kompleks dibangkan Eliminasi Gauss.

  • Membutuhkan lebih banyak langkah komputasi.

Iterasi Jacobi

Membutuhkan matriks A dengan dominansi diagonal untuk menjamin konvergensi.

# Iterasi Jacobi 
jacobi = function(A, b, tol = 1e-6, maxiter = 1000){
  n = length(b)
  x = rep(0, n) # inialisasi solusi 
  for(iter in 1:maxiter){
    x_new = rep(0, n)
    for(i in 1:n){
      s = sum(A[i, - i]*x[-i])
      x_new[i] = (b[i] - s)/A[i, i]
    }
    if(all(abs(x - x_new) < tol)){
      return(list(x = x_new, iter = iter))
    }
    x = x_new
  }
 return(list(x = x, iter = maxiter)) 
}

# Jalankan iterasi Jacobi 
jacobi_result = jacobi(A, b)
jacobi_result
## $x
## [1] 1.506024 3.132530 2.650602
## 
## $iter
## [1] 22

Kelebihan:

  • Cocok untuk matriks besar yang jarang/sparse.

  • Proses iterasi paralel lebih mudah diimplementasikan.

Kekurangan:

  • Tidak selalu konvergen, terutama jika tidak ada dominasi diagonal.

  • Membutuhkan iterasi lebih banyak dibandingkan Gauss-Seidel.

Iterasi Gauss-Seidel

Membutuhkan matriks A dengan dominansi diagonal untuk menjamin konvergensi.

# Metode Gauss-Seidel
gaussSeidel = function(A, b, epsilon, maxIterations) {
    # Extract the size of the system
    n = length(b)
    # Initialize the solution vector x
    x = numeric(n)
    # Iterate until either epsilon is reached or maxIterations is reached
    for (iteration in 1:maxIterations) {
        x_old = x
        # Update each variable using the Gauss-Seidel formula
        for (i in 1:n) {
            x[i] = (b[i] - sum(A[i, -i] * x[-i])) / A[i, i]
        }
        # Check if epsilon is reached
        if (max(abs(x - x_old)) < epsilon) {
            break
        }
    }
    # Return the solution and the number of iterations
    list(x = x, iterations = iteration)
}

# Menjalankan metode Gauss-Seidel
gauss_seidel_result = gaussSeidel(A, b, 1e-6, 1000)
cat("Gauss-Seidel Method: Solusi x y z =", gauss_seidel_result$x, "setelah", gauss_seidel_result$iterations, "iterasi.\n")
## Gauss-Seidel Method: Solusi x y z = 1.506024 3.13253 2.650602 setelah 15 iterasi.

Kelebihan:

  • Konvergensi lebih cepat dari Jacobi.

  • Efisien untuk matriks besar yang jarang.

Kekurangan:

  • Membutuhkan dominansi diagonal untuk menjadi konvergen.

Kasus Regresi Linear

data trees digunakan untuk memperkirakan volume pohon (Y) berdasarkan diameter dan tinggi pohon (X). Metode Eliminasi dan Iterasi dipakai untuk menghitung beta \(\beta\) dengan validasifungsi lm. Hasilnya diverifikasi dengan berbagai metode untuk memastikan kesesuaian solusi

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

Visualisasi

# Hubungan Diameter dan Volume Pohon
library(ggplot2)
ggplot(trees, aes(x = Girth, y = Volume)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Hubungan Diameter dan Volume Pohon",
       x = "Diameter Pohon (Girth)",
       y = "Volume Pohon")
## `geom_smooth()` using formula = 'y ~ x'

GGally::ggpairs(trees)

Terlihat terdapat hubungan linear antar peubah pada data trees. Dalam hal ini, Volume akan dilakukan sebagai variabel respons (Y).

# Bentuk data dalam sebuah matiks A
pred = cbind(itercept = 1, Girth = trees$Girth, Height = trees$Height)
head(pred)
##      itercept Girth Height
## [1,]        1   8.3     70
## [2,]        1   8.6     65
## [3,]        1   8.8     63
## [4,]        1  10.5     72
## [5,]        1  10.7     81
## [6,]        1  10.8     83
# Bentuk vektor B
resp = trees$Volume
head(resp)
## [1] 10.3 10.3 10.2 16.4 18.8 19.7
# Transformasi 
A = t(pred) %*% pred #(X'X)
b = t(pred) %*% resp #(X'y)
ab = cbind(A, b)

Regresi Linear Sederhana sebagai validasi untuk setiap metode yang akan digunakan. Base function dari regresi linear sebagai acuan.

reglin = lm(Volume ~., data = trees)
coef_reglin = reglin$coefficients 
coef_reglin
## (Intercept)       Girth      Height 
## -57.9876589   4.7081605   0.3392512

Nilai ini sebagai validasi untuk metode selanjutnya

Eliminasi Gauss

A_Gauss = gauss_elim(A,b)
A_Gauss
##          itercept    Girth    Height          
## itercept       31 410.7000 2356.0000  935.3000
## Girth           0 295.4374  311.5000 1496.6435
## Height          0   0.0000  889.5641  301.7857
coef_gauss = subs_balik(A_Gauss)
rbind(coef_gauss, coef_reglin)
##             (Intercept)    Girth    Height
## coef_gauss    -57.98766 4.708161 0.3392512
## coef_reglin   -57.98766 4.708161 0.3392512

coef_gauss dan coef_reglin memiliki output yang sama.

Eliminasi Gauss-Jordan

A_Gauss_Jordan = gauss_j_elim(A, b)
A_Gauss_Jordan
##          itercept Girth Height            
## itercept        1     0      0 -57.9876589
## Girth           0     1      0   4.7081605
## Height          0     0      1   0.3392512
coef_gauss_jordan = A_Gauss_Jordan[,4]
rbind(coef_gauss_jordan, coef_reglin)
##                    itercept    Girth    Height
## coef_gauss_jordan -57.98766 4.708161 0.3392512
## coef_reglin       -57.98766 4.708161 0.3392512

coef_gaus_jordan dan coef_ reglin memiliki output yang sama.

Gauss-Seidel

gauss_seidel_result = gaussSeidel(A, b, 1e-6, 1000)
gauss_seidel_result
## $x
## [1] -57.7752385   4.7123396   0.3357443
## 
## $iterations
## [1] 1000
coef_gauss_seidel = gauss_seidel_result$x
rbind(coef_gauss_jordan, coef_reglin)
##                    itercept    Girth    Height
## coef_gauss_jordan -57.98766 4.708161 0.3392512
## coef_reglin       -57.98766 4.708161 0.3392512

coef_gauss_jordan dan coef_reglin juga memiliki output yang sama.