For assistance Rosetta Code and Stackoverflow was used.
Gaussian Elimination
gauss_elim <- function(mat) {
row_count <- nrow(mat)
col_count <- ncol(mat)
pivot <- 1
for(curr_row in 1:row_count) {
if(col_count > pivot) {
i <- curr_row
while(mat[i, pivot] == 0) {
i <- i + 1
if(row_count == i) {
i <- curr_row
pivot <- pivot + 1
if(col_count == pivot)
return(mat)
}
}
mat <- swapRows(mat, curr_row, i)
for(j in curr_row:row_count)
if(j != curr_row) {
p <- mat[j, pivot] / mat[curr_row, pivot]
mat <- replaceRow(mat, curr_row, j, p)
}
pivot <- pivot + 1
}
}
return(mat)
}
Function to swap rows
swapRows <- function(mat, r1, r2) {
tmpr <- mat[r1,]
mat[r1,] <- mat[r2,]
mat[r2,] <- tmpr
return(mat)
}
Function to replace rows after elimination
replaceRow <- function(mat, r1, r2, p) {
mat[r2,] <- mat[r2,] - mat[r1,] * p
return(mat)
}
\[\left[\begin{array} {rrr} 1 & 1 & 3\\ 2 & -1 & 5\\ -1 & -2 & 4 \end{array}\right] = \left[\begin{array} {rrr} x_{1}\\ x_{2}\\ x_{3} \end{array}\right] = \left[\begin{array} {rrr} 1\\ 2\\ 6 \end{array}\right] \]
A <- matrix(c(1, 1, 3,
2, -1, 5,
-1, -2, 4), 3, 3, byrow=T)
b <- matrix(c(1, 2, 6), 3, 1, byrow = T)
aug_A <- cbind(A, b)
final_A <- gauss_elim(aug_A)
final_A # final product
## [,1] [,2] [,3] [,4]
## [1,] 1 1 3.000000 1
## [2,] 0 -3 -1.000000 0
## [3,] 0 0 7.333333 7
#Substitutions
x3 <- final_A[3,4] / final_A[3,3]
x2 <- (final_A[2,4] - (final_A[2,3] * x3)) / final_A[2, 2]
x1 <- final_A[1,4] - ((final_A[1,2] * x2) + (final_A[1,3] * x3)) / final_A[1,1]
cat("Solutions:", "\n x1 -> ", round(x1, 2), "\t\t\n", "x2 -> ", round(x2, 2), "\t\t\n", "x3 -> ", round(x3, 2))
## Solutions:
## x1 -> -1.55
## x2 -> -0.32
## x3 -> 0.95