Gaussian elimination to solve AX=b

2020-11-03

library(matlib)

## Warning: package 'matlib' was built under R version 4.0.3

A<-matrix(c(2,3,-2,1,2,2,4,2,3),3,3)
b<-c(1,2,3)

gaussian_elimiation<-function(A,b){
  for(J in  1 : nrow(A)){
    C = 1 / A[J,J]
    for(K in 1:nrow(A)){
      A[J,K] = A[J,K] * C;
    }
    b[J] = b[J] * C
    
    for(I in 1 : nrow(A)){
      if(I != J){ C = A[I,J]
      for(K in 1:nrow(A)){
        A[I,K] = A[I,K] - A[J,K] * C;
      }
      b[I] = b[I] - b[J] * C;
      }
      
    }
  } 
  return (list("A"=A,"b"=b))
}

gaussian_elimiation(A,b)

## $A
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
## 
## $b
## [1] -0.19354839  1.25806452  0.03225806

# check with other methods
solve(A)%*%b

##             [,1]
## [1,] -0.19354839
## [2,]  1.25806452
## [3,]  0.03225806

gaussianElimination(A, b)

##      [,1] [,2] [,3]        [,4]
## [1,]    1    0    0 -0.19354839
## [2,]    0    1    0  1.25806452
## [3,]    0    0    1  0.03225806