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.
library(nnet)## Warning: package 'nnet' was built under R version 3.2.3#A is 3x3 coefficient matrix, B are the contraints
#function solves for vector X
upper_tri_solver <- function(A, b) {
n <- nrow(A)
for(k in 1:(n-1)) {
max_index <- which.is.max(abs(A[,k]))
if (max_index != k) {
A1 <- A[k,]
A2 <- A[max_index,]
b1 <- b[k]
b2 <- b[max_index]
A[k,] <- A2
A[max_index,] <- A1
b[k] <- b2
b[max_index] <- b1
}
for (i in (k+1):n) {
multiplier <- A[i,k]/A[k,k]
A[i,k] <- 0
for (j in (k+1):n) {
A[i,j] <- A[i,j] - multiplier * A[k,j]
}
b[i] <- b[i] - multiplier * b[k]
}
}
x <- vector(mode="numeric",length=n)
x[n] <- b[n]/A[n,n]
x[n-1] <- (b[n-1] - (A[n-1,n] * x[n])) / A[n-1,n-1]
x[n-2] <- (b[n-2] - (A[n-2,n-1] * x[n-1]) - (A[n-2,n] * x[n])) / A[n-2,n-2]
return(x)
}Running the function on the test system of equations produces the correct result:
matrix_elements <- c(1,2,-1,1,-1,-2,3,5,4)
A <- matrix(matrix_elements, nrow=3, ncol=3)
b <- c(1,2,6)
upper_tri_solver(A, b)## [1] -1.5454545 -0.3181818 0.9545455