如何在R中进行高斯消除(不要使用"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

<标题> 输出

> 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

在某些情况下，我无法解出解。即使我知道数学上错误发生的原因，我也无法在R中修正错误。帮我修理一下。

Answer 1

我不得不重新排序你的矩阵A，不知道怎么做一个通用代码:

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

编辑:

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

Answer 2

很抱歉，这个回答太迟了。
我写了一个函数来进行高斯消去

注意:
i)这个函数只能处理实数，不能处理变量。
ii)在最后的返回中，我使用了函数as.fractions()。这只在MASS包中可用，您需要至少有R 3.2.5或更新的版本才能使用它。如果你愿意，你可以提交这个，但是矩阵的元素将作为小数。
iii)此函数还返回分解PA = LDU中使用的矩阵如果只需要分解PA = LU，则使用第二个参数FALSE。

希望这有帮助!

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)))
  }
}

Answer 3

高斯<- function (){A <-矩阵(c (2， - 5,4,1,4,1,1， - 4,6)， byrow = T, nrow = 3, ncol = 3)-矩阵(c (- 3,5,10)， nrow = 3, ncol = 1)-中国(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)

}

Answer 4

可以通过pracma库进行高斯消去。

要安装它，你可以输入:

install.packages("pracma")

然后，按如下方式使用:

library(pracma)
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)
rref(cbind(A, b)

结果是:

     [,1] [,2] [,3]  [,4]
[1,]    1    0    0 -51.0
[2,]    0    1    0 -25.0
[3,]    0    0    1  -6.5