How to do Gaussian elimination in R (do not use "solve")

Question

A <- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
U.pls <- cbind(A,b)
for (i in 1:p){
    for (j in (i+1):(p+1)) U.pls[i,j] <- U.pls[i,j]/U.pls[i,i]
    U.pls[i,i] <- 1
    if (i < p) {
       for (k in (i+1):p) U.pls[k,] <- 
                   U.pls[k,] - U.pls[k,i]/U.pls[i,i]*U.pls[i,]
    }
}
U.pls
x <- rep(0,p)
for (i in p:1){
  if (i < p){
     temp <- 0
     for (j in (i+1):p) temp <- temp + U.pls[i,j]*x[j]
     x[i] <- U.pls[i,p+1] - temp
  }
  else x[i] <- U.pls[i,p+1]
}
x

output

> U.pls
     [,1] [,2] [,3] [,4]
[1,]    1 -2.5    2 -1.5
[2,]    0  1.0 -Inf  Inf
[3,]    0  0.0    1  NaN
> x
[1] NaN NaN NaN

In this way, I cannot solve solution in some case. Even though I know the reason why the error happens mathematically, I cannot fix the error in R. Help me with some fix. Thank you in advance.

Solution 1

I had to reorder you matrix A, don't know how do a generic code:

A <- matrix(c(2,-5,4,1,-4,6,1,-2.5,1),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

for (i in 2:p){
 for (j in i:p) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
}

for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls

EDIT:

A <- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

i <- 2
while (i < p+1) {
 j <- i
 while (j < p+1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
  j <- j+1
 }
 while (U.pls[i,i] == 0) {
  U.pls <- rbind(U.pls[-i,],U.pls[i,])
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
 i <- i+1
}
for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls

Solution 2

Sorry if this answer is too late.
I've wrote a function to make the gaussian elimination

Note:
i) This function only works with real numbers and not with variables.
ii) In the return at the end I have used the function as.fractions().This is only available in the MASS package and you need to have at least R version 3.2.5 or newer to use it. You can ommit this is if you want but the elements of the matrixes will be as decimals.
iii) This function also returns the matrixes used in the decomposition PA = LDU If you want only the decomposition PA = LU use as a second argument FALSE.

Hope this help!

gauss_elimination = function(matrix,diagonal = T){
  #This function performs the gaussian elimination for one column
  for_one_column = function(matrix, starting_row = 1, column = 1){
    for (i in (starting_row + 1):nrow(matrix)){
      L_table[i,column] <<-  matrix[i,column]/matrix[starting_row,column]
      matrix[i,] = matrix[i,] - matrix[i,column]/
        matrix[starting_row,column]*
        matrix[starting_row,]
    }
    return(matrix)
  }

  #This function changes a row of a matrix with another one
  change_lines = function(matrix,line_1,line_2,column_1 = 1,
                          column_2 =  ncol(matrix)){
    scapegoat = matrix[line_1,column_1:column_2]
    matrix[line_1,column_1:column_2] = matrix[line_2,column_1:column_2]
    matrix[line_2,column_1:column_2] = scapegoat
    return(matrix)
  }

  #This function checks if alternations need to be made
  checking_zeroes = function(matrix,starting_row,column){
    #If pilot is not zero
    if (matrix[starting_row,column] != 0){
      return ("No alternation needed")
      #If pilot is zero
    }else{
      row_to_change = 0
      for (i in starting_row:nrow(matrix)){
        if (matrix[i,column] != 0){
          row_to_change = i
          break
        }
      }
      #If both the pilot and all the elements below it are zero
      if (row_to_change == 0){
        return("Skip this column")
        #If the pilot is zero but at least one below it is not
      }else{
        return(c("Alternation needed",row_to_change))
      }
    }
  }

  #The main program
  row_to_work = 1 ; alternation_table = diag(nrow(matrix))
  L_table = diag(nrow(matrix))
  for (column_to_work in 1:ncol(matrix)){
    a = checking_zeroes(matrix,row_to_work,column_to_work)
    if (a[1] == "Alternation needed"){
      matrix = change_lines(matrix,as.numeric(a[2]),row_to_work)
      alternation_table = change_lines(alternation_table,row_to_work,
                                       as.numeric(a[2]))
      if (as.numeric(a[2]) != 1){
        L_table = change_lines(L_table,row_to_work,
                               as.numeric(a[2]),1,column_to_work - 1)
      }
    }
    if (a[1] == "Skip this column"){
      next()
    }
    matrix = for_one_column(matrix,row_to_work,column_to_work)
    if(row_to_work + 1 == nrow(matrix)){
      break
    }
    row_to_work = row_to_work + 1
  }
  if (diagonal == FALSE){
    return(list("P" = as.fractions(alternation_table),
                "L" = as.fractions(L_table),
                "U : The matrix after the elimination"=as.fractions(matrix)))
  }else{
    D = diag(nrow(matrix))
    diag(D) = diag(matrix)
    for (i in 1:nrow(D)){
      matrix[i,] = matrix[i,]/diag(D)[i]
    }
    return(list("P" = as.fractions(alternation_table),
                "L" = as.fractions(L_table),
                "D" = as.fractions(D),
                "U : The matrix after the elimination"=as.fractions(matrix)))
  }
}

Solution 3

Gauss <- function () { A <- matrix ( c ( 2 , -5 , 4 , 1 , 4 , 1 , 1 , -4 , 6 ), byrow = T , nrow = 3 , ncol = 3 ) b <- matrix ( c ( -3 , 5 , 10 ), nrow = 3 , ncol = 1 ) U.pls <- cbind (A,b) p <- nrow ( A )

  r <- ncol ( U.pls )
  X <- matrix ( c(rep (0,p)))
  for ( i in 1 :(p-1))  {
    cons <-    U.pls[i+1,i]/ U.pls[i,i]
    for ( j in (i+1):(p))   {
        U.pls[j,] <-  U.pls[j,] -   U.pls[i,]  *   cons
            }
     }

    X [p] <- U.pls[ p,r]/ U.pls[p,p]

  for (k in (p-1):1)      {
       suma = 0
       for (l in (k+1):p) {
            suma = suma + U.pls[k,l]*  X[l]  }
          X[k] <- (U.pls[k,r]-   suma   )/ U.pls[k,k]
          }
  return (X)

}