《算法》BEYOND 部分程序 part 1

Posted cuancuancuanhao

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了《算法》BEYOND 部分程序 part 1相关的知识,希望对你有一定的参考价值。

? 书中第六章部分程序,加上自己补充的代码,包括高斯消元法求解线性方程组,高斯 - 约旦消元法求解线性方程组

● 高斯消元法求解线性方程组,将原方程转化为上三角矩阵,然后从最后一个方程开始求解

  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 }

 

以上是关于《算法》BEYOND 部分程序 part 1的主要内容,如果未能解决你的问题,请参考以下文章

《算法》BEYOND 部分程序 part 2

《算法》第五章部分程序 part 1

《算法》第六章部分程序 part 7

《算法》第三章部分程序 part 4

《算法》第三章部分程序 part 5

《算法》第五章部分程序 part 2