请求修改R代码以支持多方程增广系数矩阵并实现消元算法

需求说明 #

原增广系数矩阵生成代码仅支持3元方程组，无法适配8元方程组，需修改该代码；同时实现高斯消元和高斯-若尔当消元两个函数，每个函数返回包含以下内容的命名列表：

• variables：未知变量组成的向量
• augcoeffmatrix：最终变换后的增广系数矩阵
• solution：解向量（方程组无解或有无穷多解时为NA）

目标8元方程组 #

# 定义8元线性方程组
E1 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 8000 * x1 + 4500 * x2 + 4000 * x3 + 3000 * x4 + 2000 * x5 + 1000 * x6 + 900 * x7 + 250 * x8 + -143145000
E2 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 7800 * x1 + 6500 * x2 + 5800 * x3 + 0 * x4 + 3100 * x5 + 1600 * x6 + 1000 * x7 + 300 * x8 + -58870000
E3 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 10000 * x1 + 0 * x2 + 3100 * x3 + 0 * x4 + 2600 * x5 + 1300 * x6 + 850 * x7 + 150 * x8 + -108440000
E4 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 5200 * x1 + 3700 * x2 + 3100 * x3 + 2700 * x4 + 2400 * x5 + 1800 * x6 + 1200 * x7 + 450 * x8 + -143805000
E5 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 7700 * x1 + 7100 * x2 + 0 * x3 + 5700 * x4 + 5100 * x5 + 1300 * x6 + 950 * x7 + 95 * x8 + -181390500
E6 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 9300 * x1 + 8700 * x2 + 6100 * x3 + 5100 * x4 + 4000 * x5 + 1000 * x6 + 700 * x7 + 70 * x8 + -209273000
E7 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 6000 * x1 + 0 * x2 + 5000 * x3 + 4300 * x4 + 3000 * x5 + 1900 * x6 + 1400 * x7 + 920 * x8 + -174388000
E8 <- function(x1, x2, x3, x4, x5, x6, x7, x8) 8500 * x1 + 3700 * x2 + 4200 * x3 + 3900 * x4 + 3500 * x5 + 2400 * x6 + 1000 * x7 + 250 * x8 + -183065000

# 打包为方程组列表
system_8var <- list(E1, E2, E3, E4, E5, E6, E7, E8)

修改后的增广系数矩阵生成函数 #

原代码核心问题为矩阵列数定义错误（用方程数而非变量数定义列数），同时未处理负号分割的边界情况，修复后版本如下：

AugCoeffMatrix <- function(system) {
  # 检查所有方程的变量数是否一致
  numOfVar <- sapply(system, function(eq) length(formals(eq)))
  if (length(unique(numOfVar)) != 1) {
    stop("所有方程的未知变量数量必须一致")
  }
  n_vars <- numOfVar[1]
  n_eq <- length(system)
  
  # 提取变量名（从第一个方程提取，已验证所有方程变量数一致）
  variables <- names(formals(system[[1]]))
  
  # 初始化增广矩阵：行数=方程数，列数=变量数+1（右侧常数项）
  augcoeffmatrix <- matrix(0, nrow = n_eq, ncol = n_vars + 1)
  colnames(augcoeffmatrix) <- c(variables, "RHS")
  rownames(augcoeffmatrix) <- paste0("E", 1:n_eq)
  
  # 遍历每个方程提取系数
  for (i in 1:n_eq) {
    # 将方程转为字符串，把负号替换为"+ -"以便统一分割
    eq_str <- deparse(system[[i]])[2]
    eq_str <- gsub("-", "+ -", eq_str)
    # 按"+"分割项，去除空项和多余空格
    terms <- strsplit(eq_str, "\\+")[[1]]
    terms <- trimws(terms)
    terms <- terms[terms != ""]
    
    for (term in terms) {
      if (grepl("\\*", term, fixed = TRUE)) {
        # 分割系数与变量
        parts <- strsplit(term, "\\*", fixed = TRUE)[[1]]
        coeff <- as.numeric(trimws(parts[1]))
        var <- trimws(parts[2])
        augcoeffmatrix[i, var] <- coeff
      } else {
        # 处理常数项：移到等号右侧取反
        augcoeffmatrix[i, "RHS"] <- -as.numeric(term)
      }
    }
  }
  
  return(list(variables = variables, augcoeffmatrix = augcoeffmatrix))
}

高斯消元函数 #

gaussian_elimination <- function(system) {
  # 生成增广矩阵
  aug_result <- AugCoeffMatrix(system)
  variables <- aug_result$variables
  aug_mat <- aug_result$augcoeffmatrix
  n_eq <- nrow(aug_mat)
  n_vars <- length(variables)
  tol <- 1e-9 # 浮点精度阈值，避免浮点误差干扰判断
  
  # 前向消元：生成上三角矩阵
  for (col in 1:min(n_eq, n_vars)) {
    # 寻找当前列绝对值最大的主元行（避免除数过小）
    pivot_row <- which.max(abs(aug_mat[col:n_eq, col])) + col - 1
    # 交换主元行与当前行
    if (pivot_row != col) {
      aug_mat[c(col, pivot_row), ] <- aug_mat[c(pivot_row, col), ]
    }
    # 主元接近0则跳过（列线性相关）
    if (abs(aug_mat[col, col]) < tol) {
      next
    }
    # 消去下方所有行的当前列
    for (row in (col+1):n_eq) {
      factor <- aug_mat[row, col] / aug_mat[col, col]
      aug_mat[row, col:(n_vars+1)] <- aug_mat[row, col:(n_vars+1)] - factor * aug_mat[col, col:(n_vars+1)]
    }
  }
  
  # 判断解的情况：计算系数矩阵与增广矩阵的秩
  rank_coeff <- sum(apply(aug_mat[, 1:n_vars], 1, function(row) any(abs(row) > tol)))
  rank_aug <- sum(apply(aug_mat, 1, function(row) any(abs(row) > tol)))
  
  solution <- NA
  if (rank_coeff < rank_aug) {
    # 无解：增广矩阵秩大于系数矩阵秩
    solution <- NA
  } else if (rank_coeff < n_vars) {
    # 无穷多解：系数矩阵秩小于变量数
    solution <- NA
  } else {
    # 唯一解：回代求解
    solution <- numeric(n_vars)
    names(solution) <- variables
    for (row in n_eq:1) {
      solution[row] <- (aug_mat[row, n_vars+1] - sum(aug_mat[row, 1:n_vars] * solution)) / aug_mat[row, row]
    }
  }
  
  return(list(variables = variables, augcoeffmatrix = aug_mat, solution = solution))
}

高斯-若尔当消元函数 #

gauss_jordan_elimination <- function(system) {
  # 生成增广矩阵
  aug_result <- AugCoeffMatrix(system)
  variables <- aug_result$variables
  aug_mat <- aug_result$augcoeffmatrix
  n_eq <- nrow(aug_mat)
  n_vars <- length(variables)
  tol <- 1e-9 # 浮点精度阈值
  
  # 行变换：生成行最简形矩阵
  for (col in 1:min(n_eq, n_vars)) {
    # 寻找主元行
    pivot_row <- which.max(abs(aug_mat[col:n_eq, col])) + col - 1
    if (pivot_row != col) {
      aug_mat[c(col, pivot_row), ] <- aug_mat[c(pivot_row, col), ]
    }
    if (abs(aug_mat[col, col]) < tol) {
      next
    }
    # 主元归一化
    aug_mat[col, ] <- aug_mat[col, ] / aug_mat[col, col]
    # 消去所有其他行的当前列
    for (row in 1:n_eq) {
      if (row != col && abs(aug_mat[row, col]) > tol) {
        factor <- aug_mat[row, col]
        aug_mat[row, ] <- aug_mat[row, ] - factor * aug_mat[col, ]
      }
    }
  }
  
  # 判断解的情况
  rank_coeff <- sum(apply(aug_mat[, 1:n_vars], 1, function(row) any(abs(row) > tol)))
  rank_aug <- sum(apply(aug_mat, 1, function(row) any(abs(row) > tol)))
  
  solution <- NA
  if (rank_coeff < rank_aug) {
    solution <- NA
  } else if (rank_coeff < n_vars) {
    solution <- NA
  } else {
    # 直接提取行最简形中的解
    solution <- aug_mat[1:n_vars, n_vars+1]
    names(solution) <- variables
  }
  
  return(list(variables = variables, augcoeffmatrix = aug_mat, solution = solution))
}

测试代码 #

# 测试高斯消元
gauss_result <- gaussian_elimination(system_8var)
cat("高斯消元解：\n")
print(gauss_result$solution)

# 测试高斯-若尔当消元
gauss_jordan_result <- gauss_jordan_elimination(system_8var)
cat("\n高斯-若尔当消元解：\n")
print(gauss_jordan_result$solution)

# 验证解的正确性（代入原方程，残差接近0则正确）
verify_solution <- function(solution, system) {
  vars_list <- as.list(solution)
  residuals <- sapply(system, function(eq) do.call(eq, vars_list))
  return(residuals)
}

cat("\n解的残差：\n")
print(verify_solution(gauss_result$solution, system_8var))

内容的提问来源于stack exchange，提问作者Zack Lorton