▶ 书中第六章部分程序，加上自己补充的代码，包括高斯消元法求解线性方程组，高斯 - 约旦消元法求解线性方程组
● 高斯消元法求解线性方程组，将原方程转化为上三角矩阵，然后从最后一个方程开始求解
1 package package01; 2 3 import edu.princeton.cs.algs4.StdOut; 4 import edu.princeton.cs.algs4.StdRandom; 5 6 public class class01 7 { 8 private static final double EPSILON = 1e-8; 9 private final int m; // 方程个数 10 private final int n; // 变量个数 11 private double[][] a; // 增广矩阵 12 private double[] x; // 方程的解 13 private boolean haveSolution; // 方程是否有解 14 15 public class01(double[][] A, double[] b) 16 { 17 m = A.length; 18 n = A[0].length; 19 if (b.length != m) 20 throw new IllegalArgumentException("Dimensions disagree"); 21 a = new double[m][n + 1]; 22 x = new double[n]; 23 24 for (int i = 0; i < m; i++) // 搬运 A 和 b 到 a 中 25 { 26 for (int j = 0; j < n; j++) 27 a[i][j] = A[i][j]; 28 } 29 for (int i = 0; i < m; i++) 30 a[i][n] = b[i]; 31 32 for (int p = 0; p < Math.min(m, n); p++) // 消元计算，a[p][p] 为当前主元 33 { 34 int max = p; // 从第 p 行往下寻找主元最大的行 35 for (int i = p + 1; i < m; i++) 36 max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max; 37 38 double[] temp = a[p]; // 将最大主元行换到第 p 行 39 a[p] = a[max]; 40 a[max] = temp; 41 42 if (Math.abs(a[p][p]) <= EPSILON) // 主元为 0，退化 43 continue; 44 45 for (int i = p + 1; i < m; i++) // 用主元消去后面的所有行 46 { 47 double alpha = a[i][p] / a[p][p]; 48 for (int j = p; j <= n; j++) 49 a[i][j] -= alpha * a[p][j]; 50 } 51 } 52 53 for (int i = Math.min(n - 1, m - 1); i >= 0; i--) // 从最后一个方程开始求解 x[i] 54 { 55 double sum = 0.0; 56 for (int j = i + 1; j < n; j++) 57 sum += a[i][j] * x[j]; 58 if (Math.abs(a[i][i]) > EPSILON) 59 x[i] = (a[i][n] - sum) / a[i][i]; // 解 x[i] 60 else if (Math.abs(a[i][n] - sum) > EPSILON) // x[i] 系数为 0，但是方程右边非 0，无解 61 { 62 StdOut.println("No solution"); 63 return; 64 } 65 } 66 for (int i = n; i < m; i++) // m > n 的情形，检查剩余的方程是否有矛盾 67 { 68 double sum = 0.0; 69 for (int j = 0; j < n; j++) 70 sum += a[i][j] * x[j]; 71 if (Math.abs(a[i][n] - sum) > EPSILON) 72 { 73 StdOut.println("No solution"); 74 return; 75 } 76 } 77 78 for (int i = 0; i < m; i++) // 检验 A*x = b 79 { 80 double sum = 0.0; 81 for (int j = 0; j < n; j++) 82 sum += A[i][j] * x[j]; 83 if (Math.abs(sum - b[i]) > EPSILON) 84 StdOut.println("not feasible b[" + i + "] = " + b[i] + ", sum = " + sum); 85 } 86 haveSolution = true; 87 } 88 89 public double[] solution() 90 { 91 return haveSolution ? x : null; 92 } 93 94 private static void test(String name, double[][] A, double[] b) 95 { 96 StdOut.println("---------------------------------------------------- " + name + " ----------------------------------------------------"); 97 class01 gaussian = new class01(A, b); 98 double[] x = gaussian.solution(); 99 if (x != null) 100 { 101 for (int i = 0; i < x.length; i++) 102 StdOut.printf("%.6f ", x[i]); 103 } 104 StdOut.println(); 105 } 106 107 private static void test1() // 3 × 3 非退化矩阵，x = [-1, 2, 2] 108 { 109 double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 } }; 110 double[] b = { 4, 2, 36 }; 111 test("test 1 (3-by-3 system, nonsingular)", A, b); 112 } 113 114 private static void test2() // 3 × 3 无解 115 { 116 double[][] A = { { 1, -3, 1 },{ 2, -8, 8 },{ 3, -11, 9 } }; 117 double[] b = { 4, -2, 9 }; 118 test("test 2 (3-by-3 system, nonsingular)", A, b); 119 } 120 121 private static void test3() // 5 × 5 无解 122 { 123 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 } }; 124 double[] b = { 4, 4, 9, -6, 5 }; 125 test("test 3 (5-by-5 system, no solutions)", A, b); 126 } 127 128 private static void test4() // 5 × 5 无穷多组解，x = [-6.25, -4.5, 0, 0, 1] 129 { 130 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 } }; 131 double[] b = { 4, 4, 9, -5, 5 }; 132 test("test 4 (5-by-5 system, infinitely many solutions)", A, b); 133 } 134 135 private static void test5() // 4 × 3 列满秩多解，x = [-1, 2, 2, ?] 136 { 137 double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } }; 138 double[] b = { 4, 2, 36, 42 }; 139 test("test 7 (4-by-3 system, full rank)", A, b); 140 } 141 142 private static void test6() // 4 × 3 列满秩无解 143 { 144 double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } }; 145 double[] b = { 4, 2, 36, 40 }; 146 test("test 8 (4-by-3 system, no solution)", A, b); 147 } 148 149 private static void test7() // 3×4 满秩，x = [3, -1, -2, 0] 150 { 151 double[][] A = { { 1, -3, 1, 1 },{ 2, -8, 8, 2 },{ -6, 3, -15, 3 } }; 152 double[] b = { 4, -2, 9 }; 153 test("test 9 (3-by-4 system, full rank)", A, b); 154 } 155 156 public static void main(String[] args) 157 { 158 test1(); 159 test2(); 160 test3(); 161 test4(); 162 test5(); 163 test6(); 164 test7(); 165 int n = Integer.parseInt(args[0]); // 随机测试 166 double[][] A = new double[n][n]; 167 for (int i = 0; i < n; i++) 168 { 169 for (int j = 0; j < n; j++) 170 A[i][j] = StdRandom.uniform(1000); 171 } 172 double[] b = new double[n]; 173 for (int i = 0; i < n; i++) 174 b[i] = StdRandom.uniform(1000); 175 test(n + "-by-" + n + " (probably nonsingular)", A, b); 176 } 177 }
● 高斯 - 约旦消元法求解线性方程组，将原方程转化为约旦标准型，无解时产生一个 yA == 0，yb != 0 的解
1 package package01; 2 3 import edu.princeton.cs.algs4.StdOut; 4 import edu.princeton.cs.algs4.StdRandom; 5 6 public class class01 7 { 8 private static final double EPSILON = 1e-8; 9 private final int n; // 方程个数等于未知数个数 10 private double[][] a; 11 private double[] x; 12 private boolean haveSolution = true; // 方程是否有解 13 14 public class01(double[][] A, double[] b) 15 { 16 n = b.length; 17 a = new double[n][n + n + 1]; 18 x = new double[n]; 19 20 for (int i = 0; i < n; i++) // 搬运 21 { 22 for (int j = 0; j < n; j++) 23 a[i][j] = A[i][j]; 24 } 25 for (int i = 0; i < n; i++) // 有解时 a[0:n-1][n:n*2] 最终为 A 的逆，否则某行为扩展解 26 a[i][n + i] = 1.0; 27 for (int i = 0; i < n; i++) 28 a[i][n + n] = b[i]; 29 30 for (int p = 0; p < n; p++) // 消元计算 31 { 32 int max = p; 33 for (int i = p + 1; i < n; i++) 34 max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max; 35 36 double[] temp = a[p]; 37 a[p] = a[max]; 38 a[max] = temp; 39 40 if (Math.abs(a[p][p]) <= EPSILON) 41 continue; 42 43 for (int i = 0; i < n; i++) // 高斯-约旦消元 44 { 45 double alpha = a[i][p] / a[p][p]; 46 for (int j = 0; j <= n + n; j++) 47 a[i][j] -= (i != p) ? alpha * a[p][j] : 0; 48 } 49 for (int j = 0; j <= n + n; j++) // 主行归一化 50 a[p][j] /= (j != p) ? a[p][p] : 1; 51 for (int i = 0; i < n; i++) // 被消去行的主元所在列元素强制归 0，主元归 1 52 a[i][p] = 0.0; 53 a[p][p] = 1.0; 54 } 55 56 for (int i = 0; i < n; i++) // 求解 x 57 { 58 if (Math.abs(a[i][i]) > EPSILON) 59 x[i] = a[i][n + n] / a[i][i]; 60 else if (Math.abs(a[i][n + n]) > EPSILON) 61 { 62 StdOut.println("No solution"); 63 haveSolution = false; 64 break; 65 } 66 } 67 if (!haveSolution) // 无解，输出约旦阵中第 i 行的结果 68 { 69 for (int i = 0; i < n; i++) 70 { 71 if ((Math.abs(a[i][i]) <= EPSILON) && (Math.abs(a[i][n + n]) > EPSILON)) 72 { 73 for (int j = 0; j < n; j++) 74 x[j] = a[i][n + j]; 75 } 76 } 77 } 78 79 if (haveSolution) // 检验结果 80 { 81 for (int i = 0; i < n; i++) // 有解则检验 A*x = b 82 { 83 double sum = 0.0; 84 for (int j = 0; j < n; j++) 85 sum += A[i][j] * x[j]; 86 if (Math.abs(sum - b[i]) > EPSILON) 87 StdOut.printf("not feasible b[%d] = %8.3f, sum = %8.3f ", i, b[i], sum); 88 } 89 } 90 else 91 { 92 for (int j = 0; j < n; j++) // 无解则检验 yA = 0, yb != 0 93 { 94 double sum = 0.0; 95 for (int i = 0; i < n; i++) 96 sum += A[i][j] * x[i]; 97 if (Math.abs(sum) > EPSILON) 98 StdOut.printf("invalid certificate of infeasibility sum = %8.3f ", sum); 99 } 100 double sum = 0.0; 101 for (int i = 0; i < n; i++) 102 sum += x[i] * b[i]; 103 if (Math.abs(sum) < EPSILON) 104 StdOut.printf("invalid certificate of infeasibility yb = %8.3f ", sum); 105 } 106 } 107 108 public double[] solution() 109 { 110 return x; 111 } 112 113 public boolean haveActualSolution() 114 { 115 return haveSolution; 116 } 117 118 private static void test(String name, double[][] A, double[] b) 119 { 120 StdOut.println("---------------------------------------------------- " + name + " ----------------------------------------------------"); 121 class01 gaussian = new class01(A, b); 122 double[] x = gaussian.solution(); 123 StdOut.println(gaussian.haveActualSolution() ? "Solution:" : "Certificate of infeasibility"); 124 for (int i = 0; i < x.length; i++) 125 StdOut.printf("%10.6f ", x[i]); 126 StdOut.println(); 127 } 128 129 private static void test1() // 3×3 非退化，x = [ -1, 2, 2 ] 130 { 131 double[][] A = { { 0, 1, 1 },{ 2, 4, -2 },{ 0, 3, 15 } }; 132 double[] b = { 4, 2, 36 }; 133 test("test 1", A, b); 134 } 135 136 private static void test2() // 3×3 无解，x = [ 1, 0, 1/3 ] 137 { 138 double[][] A = { { 2, -1, 1 },{ 3, 2, -4 },{ -6, 3, -3 } }; 139 double[] b = { 1, 4, 2 }; 140 test("test 5", A, b); 141 } 142 143 private static void test3() // 5×5 非退化，x = [ -6.25, -4.5, 0, 0, 1 ] 144 { 145 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 } }; 146 double[] b = { 4, 4, 9, -5, 5 }; 147 test("test 4", A, b); 148 } 149 150 private static void test4() // 5×5 无解，x = [ -1, 0, 1, 1, 0 ] 151 { 152 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 } }; 153 double[] b = { 4, 4, 9, -6, 5 }; 154 test("test 3", A, b); 155 } 156 157 private static void test5() // 3×3 无穷多组解，x = [ -1.375, 1.625, 0 ] 158 { 159 double[][] A = { { 1, -1, 2 },{ 4, 4, -2 },{ -2, 2, -4 } }; 160 double[] b = { -3, 1, 6 }; 161 test("test 6 (infinitely many solutions)", A, b); 162 } 163 164 public static void main(String[] args) 165 { 166 test1(); 167 test2(); 168 test3(); 169 test4(); 170 test5(); 171 int n = Integer.parseInt(args[0]); 172 double[][] A = new double[n][n]; 173 double[] b = new double[n]; 174 for (int i = 0; i < n; i++) 175 { 176 for (int j = 0; j < n; j++) 177 A[i][j] = StdRandom.uniform(1000); 178 } 179 for (int i = 0; i < n; i++) 180 b[i] = StdRandom.uniform(1000); 181 test("random " + n + "-by-" + n + " (likely full rank)", A, b); 182 183 for (int i = 0; i < n - 1; i++) 184 { 185 for (int j = 0; j < n; j++) 186 A[i][j] = StdRandom.uniform(1000); 187 } 188 for (int i = 0; i < n - 1; i++) 189 { 190 double alpha = StdRandom.uniform(11) - 5.0; 191 for (int j = 0; j < n; j++) 192 A[n - 1][j] += alpha * A[i][j]; 193 } 194 for (int i = 0; i < n; i++) 195 b[i] = StdRandom.uniform(1000); 196 test("random " + n + "-by-" + n + " (likely infeasible)", A, b); 197 } 198 }