solve static method

List<double?>? solve(List<List<double>> a, List<double> b)

Solves a system of equations using Gaussian Elimination. In order to avoid overhead the algorithm runs in-place on A - if A should not be modified the client must supply a copy.

@param a an nxn matrix in row/column order )modified by this method) @param b a vector of length n

@return a vector containing the solution (if any) or null if the system has no or no unique solution

@throws IllegalArgumentException if the matrix is the wrong size

Implementation

static List<double?>? solve(List<List<double>> a, List<double> b) {
  int n = b.length;
  if (a.length != n || a[0].length != n)
    throw ArgumentError("Matrix A is incorrectly sized");

  // Use Gaussian Elimination with partial pivoting.
  // Iterate over each row
  for (int i = 0; i < n; i++) {
    // Find the largest pivot in the rows below the current one.
    int maxElementRow = i;
    for (int j = i + 1; j < n; j++)
      if (a[j][i].abs() > a[maxElementRow][i].abs()) maxElementRow = j;

    if (a[maxElementRow][i] == 0.0) return null;

    // Exchange current row and maxElementRow in A and b.
    swapRows(a, i, maxElementRow);
    swapRowsList(b, i, maxElementRow);

    // Eliminate using row i
    for (int j = i + 1; j < n; j++) {
      double rowFactor = a[j][i] / a[i][i];
      for (int k = n - 1; k >= i; k--) a[j][k] -= a[i][k] * rowFactor;
      b[j] -= b[i] * rowFactor;
    }
  }

  /**
   * A is now (virtually) in upper-triangular form.
   * The solution vector is determined by back-substitution.
   */
  List<double> solution = List.filled(n, 0.0);
  for (int j = n - 1; j >= 0; j--) {
    double t = 0.0;
    for (int k = j + 1; k < n; k++) t += a[j][k] * solution[k];
    solution[j] = (b[j] - t) / a[j][j];
  }
  return solution;
}