▶ 书中第六章部分程序，加上自己补充的代码，包括高斯消元法求解线性方程组，高斯 - 约旦消元法求解线性方程组
● 高斯消元法求解线性方程组，将原方程转化为上三角矩阵，然后从最后一个方程开始求解
package package01; import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom; public class class01
{
private static final double EPSILON = 1e-8;
private final int m; // 方程个数
private final int n; // 变量个数
private double[][] a; // 增广矩阵
private double[] x; // 方程的解
private boolean haveSolution; // 方程是否有解 public class01(double[][] A, double[] b)
{
m = A.length;
n = A[0].length;
if (b.length != m)
throw new IllegalArgumentException("Dimensions disagree");
a = new double[m][n + 1];
x = new double[n]; for (int i = 0; i < m; i++) // 搬运 A 和 b 到 a 中
{
for (int j = 0; j < n; j++)
a[i][j] = A[i][j];
}
for (int i = 0; i < m; i++)
a[i][n] = b[i]; for (int p = 0; p < Math.min(m, n); p++) // 消元计算，a[p][p] 为当前主元
{
int max = p; // 从第 p 行往下寻找主元最大的行
for (int i = p + 1; i < m; i++)
max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max; double[] temp = a[p]; // 将最大主元行换到第 p 行
a[p] = a[max];
a[max] = temp; if (Math.abs(a[p][p]) <= EPSILON) // 主元为 0，退化
continue; for (int i = p + 1; i < m; i++) // 用主元消去后面的所有行
{
double alpha = a[i][p] / a[p][p];
for (int j = p; j <= n; j++)
a[i][j] -= alpha * a[p][j];
}
} for (int i = Math.min(n - 1, m - 1); i >= 0; i--) // 从最后一个方程开始求解 x[i]
{
double sum = 0.0;
for (int j = i + 1; j < n; j++)
sum += a[i][j] * x[j];
if (Math.abs(a[i][i]) > EPSILON)
x[i] = (a[i][n] - sum) / a[i][i]; // 解 x[i]
else if (Math.abs(a[i][n] - sum) > EPSILON) // x[i] 系数为 0，但是方程右边非 0，无解
{
StdOut.println("No solution");
return;
}
}
for (int i = n; i < m; i++) // m > n 的情形，检查剩余的方程是否有矛盾
{
double sum = 0.0;
for (int j = 0; j < n; j++)
sum += a[i][j] * x[j];
if (Math.abs(a[i][n] - sum) > EPSILON)
{
StdOut.println("No solution");
return;
}
} for (int i = 0; i < m; i++) // 检验 A*x = b
{
double sum = 0.0;
for (int j = 0; j < n; j++)
sum += A[i][j] * x[j];
if (Math.abs(sum - b[i]) > EPSILON)
StdOut.println("not feasible\n b[" + i + "] = " + b[i] + ", sum = " + sum);
}
haveSolution = true;
} public double[] solution()
{
return haveSolution ? x : null;
} private static void test(String name, double[][] A, double[] b)
{
StdOut.println("----------------------------------------------------\n" + name + "\n----------------------------------------------------");
class01 gaussian = new class01(A, b);
double[] x = gaussian.solution();
if (x != null)
{
for (int i = 0; i < x.length; i++)
StdOut.printf("%.6f\n", x[i]);
}
StdOut.println();
} private static void test1() // 3 × 3 非退化矩阵，x = [-1, 2, 2]
{
double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 } };
double[] b = { 4, 2, 36 };
test("test 1 (3-by-3 system, nonsingular)", A, b);
} private static void test2() // 3 × 3 无解
{
double[][] A = { { 1, -3, 1 },{ 2, -8, 8 },{ 3, -11, 9 } };
double[] b = { 4, -2, 9 };
test("test 2 (3-by-3 system, nonsingular)", A, b);
} private static void test3() // 5 × 5 无解
{
double[][] A = { { 2, -3, -1, 2, 3 },{ 4, -4, -1, 4, 11 },{ 2, -5, -2, 2, -1 },{ 0, 2, 1, 0, 4 },{ -4, 6, 0, 0, 7 } };
double[] b = { 4, 4, 9, -6, 5 };
test("test 3 (5-by-5 system, no solutions)", A, b);
} private static void test4() // 5 × 5 无穷多组解，x = [-6.25, -4.5, 0, 0, 1]
{
double[][] A = { { 2, -3, -1, 2, 3 },{ 4, -4, -1, 4, 11 },{ 2, -5, -2, 2, -1 },{ 0, 2, 1, 0, 4 },{ -4, 6, 0, 0, 7 } };
double[] b = { 4, 4, 9, -5, 5 };
test("test 4 (5-by-5 system, infinitely many solutions)", A, b);
} private static void test5() // 4 × 3 列满秩多解，x = [-1, 2, 2, ?]
{
double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } };
double[] b = { 4, 2, 36, 42 };
test("test 7 (4-by-3 system, full rank)", A, b);
} private static void test6() // 4 × 3 列满秩无解
{
double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } };
double[] b = { 4, 2, 36, 40 };
test("test 8 (4-by-3 system, no solution)", A, b);
} private static void test7() // 3×4 满秩，x = [3, -1, -2, 0]
{
double[][] A = { { 1, -3, 1, 1 },{ 2, -8, 8, 2 },{ -6, 3, -15, 3 } };
double[] b = { 4, -2, 9 };
test("test 9 (3-by-4 system, full rank)", A, b);
} public static void main(String[] args)
{
test1();
test2();
test3();
test4();
test5();
test6();
test7();
int n = Integer.parseInt(args[0]); // 随机测试
double[][] A = new double[n][n];
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
A[i][j] = StdRandom.uniform(1000);
}
double[] b = new double[n];
for (int i = 0; i < n; i++)
b[i] = StdRandom.uniform(1000);
test(n + "-by-" + n + " (probably nonsingular)", A, b);
}
}
● 高斯 - 约旦消元法求解线性方程组，将原方程转化为约旦标准型，无解时产生一个 yA == 0，yb != 0 的解
package package01; import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom; public class class01
{
private static final double EPSILON = 1e-8;
private final int n; // 方程个数等于未知数个数
private double[][] a;
private double[] x;
private boolean haveSolution = true; // 方程是否有解 public class01(double[][] A, double[] b)
{
n = b.length;
a = new double[n][n + n + 1];
x = new double[n]; for (int i = 0; i < n; i++) // 搬运
{
for (int j = 0; j < n; j++)
a[i][j] = A[i][j];
}
for (int i = 0; i < n; i++) // 有解时 a[0:n-1][n:n*2] 最终为 A 的逆，否则某行为扩展解
a[i][n + i] = 1.0;
for (int i = 0; i < n; i++)
a[i][n + n] = b[i]; for (int p = 0; p < n; p++) // 消元计算
{
int max = p;
for (int i = p + 1; i < n; i++)
max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max; double[] temp = a[p];
a[p] = a[max];
a[max] = temp; if (Math.abs(a[p][p]) <= EPSILON)
continue; for (int i = 0; i < n; i++) // 高斯-约旦消元
{
double alpha = a[i][p] / a[p][p];
for (int j = 0; j <= n + n; j++)
a[i][j] -= (i != p) ? alpha * a[p][j] : 0;
}
for (int j = 0; j <= n + n; j++) // 主行归一化
a[p][j] /= (j != p) ? a[p][p] : 1;
for (int i = 0; i < n; i++) // 被消去行的主元所在列元素强制归 0，主元归 1
a[i][p] = 0.0;
a[p][p] = 1.0;
} for (int i = 0; i < n; i++) // 求解 x
{
if (Math.abs(a[i][i]) > EPSILON)
x[i] = a[i][n + n] / a[i][i];
else if (Math.abs(a[i][n + n]) > EPSILON)
{
StdOut.println("No solution");
haveSolution = false;
break;
}
}
if (!haveSolution) // 无解，输出约旦阵中第 i 行的结果
{
for (int i = 0; i < n; i++)
{
if ((Math.abs(a[i][i]) <= EPSILON) && (Math.abs(a[i][n + n]) > EPSILON))
{
for (int j = 0; j < n; j++)
x[j] = a[i][n + j];
}
}
} if (haveSolution) // 检验结果
{
for (int i = 0; i < n; i++) // 有解则检验 A*x = b
{
double sum = 0.0;
for (int j = 0; j < n; j++)
sum += A[i][j] * x[j];
if (Math.abs(sum - b[i]) > EPSILON)
StdOut.printf("not feasible\nb[%d] = %8.3f, sum = %8.3f\n", i, b[i], sum);
}
}
else
{
for (int j = 0; j < n; j++) // 无解则检验 yA = 0, yb != 0
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i][j] * x[i];
if (Math.abs(sum) > EPSILON)
StdOut.printf("invalid certificate of infeasibility\nsum = %8.3f\n", sum);
}
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += x[i] * b[i];
if (Math.abs(sum) < EPSILON)
StdOut.printf("invalid certificate of infeasibility\nyb = %8.3f\n", sum);
}
} public double[] solution()
{
return x;
} public boolean haveActualSolution()
{
return haveSolution;
} private static void test(String name, double[][] A, double[] b)
{
StdOut.println("----------------------------------------------------\n" + name + "\n----------------------------------------------------");
class01 gaussian = new class01(A, b);
double[] x = gaussian.solution();
StdOut.println(gaussian.haveActualSolution() ? "Solution:" : "Certificate of infeasibility");
for (int i = 0; i < x.length; i++)
StdOut.printf("%10.6f\n", x[i]);
StdOut.println();
} private static void test1() // 3×3 非退化，x = [ -1, 2, 2 ]
{
double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 } };
double[] b = { 4, 2, 36 };
test("test 1", A, b);
} private static void test2() // 3×3 无解，x = [ 1, 0, 1/3 ]
{
double[][] A = { { 2, -1, 1 },{ 3, 2, -4 },{ -6, 3, -3 } };
double[] b = { 1, 4, 2 };
test("test 5", A, b);
} private static void test3() // 5×5 非退化，x = [ -6.25, -4.5, 0, 0, 1 ]
{
double[][] A = { { 2, -3, -1, 2, 3 },{ 4, -4, -1, 4, 11 },{ 2, -5, -2, 2, -1 },{ 0, 2, 1, 0, 4 },{ -4, 6, 0, 0, 7 } };
double[] b = { 4, 4, 9, -5, 5 };
test("test 4", A, b);
} private static void test4() // 5×5 无解，x = [ -1, 0, 1, 1, 0 ]
{
double[][] A = { { 2, -3, -1, 2, 3 },{ 4, -4, -1, 4, 11 },{ 2, -5, -2, 2, -1 },{ 0, 2, 1, 0, 4 },{ -4, 6, 0, 0, 7 } };
double[] b = { 4, 4, 9, -6, 5 };
test("test 3", A, b);
} private static void test5() // 3×3 无穷多组解，x = [ -1.375, 1.625, 0 ]
{
double[][] A = { { 1, -1, 2 },{ 4, 4, -2 },{ -2, 2, -4 } };
double[] b = { -3, 1, 6 };
test("test 6 (infinitely many solutions)", A, b);
} public static void main(String[] args)
{
test1();
test2();
test3();
test4();
test5();
int n = Integer.parseInt(args[0]);
double[][] A = new double[n][n];
double[] b = new double[n];
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
A[i][j] = StdRandom.uniform(1000);
}
for (int i = 0; i < n; i++)
b[i] = StdRandom.uniform(1000);
test("random " + n + "-by-" + n + " (likely full rank)", A, b); for (int i = 0; i < n - 1; i++)
{
for (int j = 0; j < n; j++)
A[i][j] = StdRandom.uniform(1000);
}
for (int i = 0; i < n - 1; i++)
{
double alpha = StdRandom.uniform(11) - 5.0;
for (int j = 0; j < n; j++)
A[n - 1][j] += alpha * A[i][j];
}
for (int i = 0; i < n; i++)
b[i] = StdRandom.uniform(1000);
test("random " + n + "-by-" + n + " (likely infeasible)", A, b);
}
}