4 
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 中的错​​误。帮我解决一些问题。先感谢您。

4

4 回答 4

3

我不得不重新排序你的矩阵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
于 2013-04-16T19:15:54.250 回答
1

抱歉，如果这个答案为时已晚。
我写了一个函数来进行高斯消除

注意：
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)))
  }
}
于 2016-05-15T08:21:04.220 回答
1

高斯 <- 函数 () { A <- 矩阵 ( c ( 2 , -5 , 4 , 1 , 4 , 1 , 1 , -4 , 6 ), byrow = T , nrow = 3 , ncol = 3 ) b <-矩阵 ( 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)

}

于 2016-06-06T05:38:32.643 回答
0

您可以通过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
于 2017-12-16T12:42:29.780 回答