Problem Set 1

2. Problem set 2

Set up a system of equations with 3 variables and 3 constraints and solve for x. Please write a function in R that will take two variables (matrix A & constraint vector b) and solve using elimination. Your function should produce the right answer for the system of equations for any 3-variable, 3-equation system. You don’t have to worry about degenerate cases and can safely assume that the function will only be tested with a system of equations that has a solution. Please note that you do have to worry about zero pivots, though. Please note that you should not use the built-in function solve to solve this system or use matrix inverses. The approach that you should employ is to construct an Upper Triangular Matrix and then back-substitute to get the solution. Alternatively, you can augment the matrix A with vector b and jointly apply the Gauss Jordan elimination procedure.

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)
}
Please test it with the system below and it should produce a solution \(x = [−1.55,−0.32,0.95]\)

\[\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