פתרון באמצעות אלימינציה

function pluSolve(matrix, rhs, relativeTolerance=1e-12) {
  const n = matrix.length;
  if (!n || rhs.length !== n || matrix.some((row) => row.length !== n)) {
    throw new Error("A must be square and match b");
  }

  const U = matrix.map((row) => row.map(Number));
  const rowScales = U.map((row) => Math.max(...row.map(Math.abs)));
  if (rowScales.some((scale) => scale === 0)) {
    throw new Error("no unique solution");
  }
  const b = rhs.map(Number);
  const L = Array.from({ length: n }, () => Array(n).fill(0));
  const P = Array.from(
    { length: n },
    (_, i) => Array.from({ length: n }, (_, j) => Number(i === j)),
  );

  for (let k = 0; k < n; k += 1) {
    let pivot = k;
    for (let i = k + 1; i < n; i += 1) {
      const candidateRatio = Math.abs(U[i][k]) / rowScales[i];
      const pivotRatio = Math.abs(U[pivot][k]) / rowScales[pivot];
      if (candidateRatio > pivotRatio) pivot = i;
    }
    if (U[pivot][k] === 0) throw new Error("no unique solution");
    if (Math.abs(U[pivot][k]) <= relativeTolerance*rowScales[pivot]) {
      throw new Error("nonzero pivot is numerically unsafe at this tolerance");
    }
    if (pivot !== k) {
      [U[k], U[pivot]] = [U[pivot], U[k]];
      [b[k], b[pivot]] = [b[pivot], b[k]];
      [P[k], P[pivot]] = [P[pivot], P[k]];
      [rowScales[k], rowScales[pivot]] = [rowScales[pivot], rowScales[k]];
      for (let j = 0; j < k; j += 1) {
        [L[k][j], L[pivot][j]] = [L[pivot][j], L[k][j]];
      }
    }

    L[k][k] = 1;
    for (let i = k + 1; i < n; i += 1) {
      const multiplier = U[i][k] / U[k][k];
      L[i][k] = multiplier;
      for (let j = k; j < n; j += 1) {
        U[i][j] -= multiplier * U[k][j];
      }
      b[i] -= multiplier * b[k];
    }
  }

  const x = Array(n).fill(0);
  for (let i = n - 1; i >= 0; i -= 1) {
    const known = U[i].slice(i + 1)
      .reduce((sum, value, j) => sum + value*x[i + 1 + j], 0);
    x[i] = (b[i] - known) / U[i][i];
  }
  return { x, P, L, U };
}

function multiplyMatrices(A, B) {
  return A.map((row) =>
    B[0].map((_, j) =>
      row.reduce((sum, value, k) => sum + value*B[k][j], 0),
    ),
  );
}

const A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]];
const result = pluSolve(A, [5, -2, 9]);
const PA = multiplyMatrices(result.P, A);
const LU = multiplyMatrices(result.L, result.U);
const factorizationError = Math.max(
  ...PA.flatMap((row, i) => row.map((value, j) => Math.abs(value - LU[i][j]))),
);
console.log({ ...result, factorizationError });
console.log(pluSolve([[1, 0], [0, 1e-13]], [1, 2e-13]).x);

def plu_solve(matrix, rhs, relative_tolerance=1e-12):
    n = len(matrix)
    if not n or len(rhs) != n or any(len(row) != n for row in matrix):
        raise ValueError("A must be square and match b")

    U = [[float(value) for value in row] for row in matrix]
    row_scales = [max(abs(value) for value in row) for row in U]
    if any(scale == 0 for scale in row_scales):
        raise ValueError("no unique solution")
    b = [float(value) for value in rhs]
    L = [[0.0]*n for _ in range(n)]
    P = [[float(i == j) for j in range(n)] for i in range(n)]

    for k in range(n):
        pivot = max(
            range(k, n),
            key=lambda i: abs(U[i][k])/row_scales[i],
        )
        if U[pivot][k] == 0:
            raise ValueError("no unique solution")
        if abs(U[pivot][k]) <= relative_tolerance*row_scales[pivot]:
            raise ArithmeticError(
                "nonzero pivot is numerically unsafe at this tolerance"
            )
        if pivot != k:
            U[k], U[pivot] = U[pivot], U[k]
            b[k], b[pivot] = b[pivot], b[k]
            P[k], P[pivot] = P[pivot], P[k]
            row_scales[k], row_scales[pivot] = row_scales[pivot], row_scales[k]
            for j in range(k):
                L[k][j], L[pivot][j] = L[pivot][j], L[k][j]

        L[k][k] = 1.0
        for i in range(k+1, n):
            multiplier = U[i][k] / U[k][k]
            L[i][k] = multiplier
            for j in range(k, n):
                U[i][j] -= multiplier*U[k][j]
            b[i] -= multiplier*b[k]

    x = [0.0]*n
    for i in range(n-1, -1, -1):
        known = sum(U[i][j]*x[j] for j in range(i+1, n))
        x[i] = (b[i] - known) / U[i][i]
    return {"x": x, "P": P, "L": L, "U": U}

def multiply_matrices(A, B):
    return [
        [
            sum(a*b for a, b in zip(row, column))
            for column in zip(*B)
        ]
        for row in A
    ]

A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]]
result = plu_solve(A, [5, -2, 9])
PA = multiply_matrices(result["P"], A)
LU = multiply_matrices(result["L"], result["U"])
result["factorization_error"] = max(
    abs(value - LU[i][j])
    for i, row in enumerate(PA)
    for j, value in enumerate(row)
)
print(result)
print(plu_solve([[1, 0], [0, 1e-13]], [1, 2e-13])["x"])