|
| 1 | +using System; |
| 2 | +using System.Collections; |
| 3 | +using System.Data; |
| 4 | + |
| 5 | +namespace Gauss { |
| 6 | + public class GaussSolutionNotFound : Exception { |
| 7 | + public GaussSolutionNotFound(string msg) |
| 8 | + : base("Ðåøåíèå íå ìîæåò áûòü íàéäåíî: \r\n" + msg) { |
| 9 | + } |
| 10 | + } |
| 11 | + |
| 12 | + public class LinearSystem { |
| 13 | + private double[,] initial_a_matrix; |
| 14 | + private double[,] a_matrix; // ìàòðèöà A |
| 15 | + private double[] x_vector; // âåêòîð íåèçâåñòíûõ x |
| 16 | + private double[] initial_b_vector; |
| 17 | + private double[] b_vector; // âåêòîð b |
| 18 | + private double eps; // ïîðÿäîê òî÷íîñòè äëÿ ñðàâíåíèÿ âåùåñòâåííûõ ÷èñåë |
| 19 | + private int size; // ðàçìåðíîñòü çàäà÷è |
| 20 | + |
| 21 | + |
| 22 | + public LinearSystem(double[,] a_matrix, double[] b_vector) |
| 23 | + : this(a_matrix, b_vector, 0.0001) { |
| 24 | + } |
| 25 | + public LinearSystem(double[,] a_matrix, double[] b_vector, double eps) { |
| 26 | + if (a_matrix == null || b_vector == null) |
| 27 | + throw new ArgumentNullException("Îäèí èç ïàðàìåòðîâ ðàâåí null."); |
| 28 | + |
| 29 | + int b_length = b_vector.Length; |
| 30 | + int a_length = a_matrix.Length; |
| 31 | + if (a_length != b_length * b_length) |
| 32 | + throw new ArgumentException(@"Êîëè÷åñòâî ñòðîê è ñòîëáöîâ â ìàòðèöå A äîëæíî ñîâïàäàòü ñ êîëè÷åñòâîì ýëåìåíòðîâ â âåêòîðå B."); |
| 33 | + |
| 34 | + this.initial_a_matrix = a_matrix; // çàïîìèíàåì èñõîäíóþ ìàòðèöó |
| 35 | + this.a_matrix = (double[,])a_matrix.Clone(); // ñ å¸ êîïèåé áóäåì ïðîèçâîäèòü âû÷èñëåíèÿ |
| 36 | + this.initial_b_vector = b_vector; // çàïîìèíàåì èñõîäíûé âåêòîð |
| 37 | + this.b_vector = (double[])b_vector.Clone(); // ñ åãî êîïèåé áóäåì ïðîèçâîäèòü âû÷èñëåíèÿ |
| 38 | + this.x_vector = new double[b_length]; |
| 39 | + this.size = b_length; |
| 40 | + this.eps = eps; |
| 41 | + |
| 42 | + GaussSolve(); |
| 43 | + } |
| 44 | + |
| 45 | + public double[] XVector { |
| 46 | + get { |
| 47 | + return x_vector; |
| 48 | + } |
| 49 | + } |
| 50 | + |
| 51 | + // èíèöèàëèçàöèÿ ìàññèâà èíäåêñîâ ñòîëáöîâ |
| 52 | + private int[] InitIndex() { |
| 53 | + int[] index = new int[size]; |
| 54 | + for (int i = 0; i < index.Length; ++i) |
| 55 | + index[i] = i; |
| 56 | + return index; |
| 57 | + } |
| 58 | + |
| 59 | + // ïîèñê ãëàâíîãî ýëåìåíòà â ìàòðèöå |
| 60 | + private double FindR(int row, int[] index) { |
| 61 | + int max_index = row; |
| 62 | + double max = a_matrix[row, index[max_index]]; |
| 63 | + double max_abs = Math.Abs(max); |
| 64 | + for (int cur_index = row + 1; cur_index < size; ++cur_index) { |
| 65 | + double cur = a_matrix[row, index[cur_index]]; |
| 66 | + double cur_abs = Math.Abs(cur); |
| 67 | + if (cur_abs > max_abs) { |
| 68 | + max_index = cur_index; |
| 69 | + max = cur; |
| 70 | + max_abs = cur_abs; |
| 71 | + } |
| 72 | + } |
| 73 | + |
| 74 | + if (max_abs < eps) { |
| 75 | + if (Math.Abs(b_vector[row]) > eps) |
| 76 | + throw new GaussSolutionNotFound("Ñèñòåìà óðàâíåíèé íåñîâìåñòíà."); |
| 77 | + else |
| 78 | + throw new GaussSolutionNotFound("Ñèñòåìà óðàâíåíèé èìååò ìíîæåñòâî ðåøåíèé."); |
| 79 | + } |
| 80 | + |
| 81 | + // ìåíÿåì ìåñòàìè èíäåêñû ñòîëáöîâ |
| 82 | + int temp = index[row]; |
| 83 | + index[row] = index[max_index]; |
| 84 | + index[max_index] = temp; |
| 85 | + |
| 86 | + return max; |
| 87 | + } |
| 88 | + |
| 89 | + // ðåøåíèå ìåòîäîì Ãàóññà |
| 90 | + private void GaussSolve() { |
| 91 | + int[] index = InitIndex(); |
| 92 | + GaussForwardStroke(index); |
| 93 | + GaussBackwardStroke(index); |
| 94 | + GaussDiscrepancy(); |
| 95 | + } |
| 96 | + |
| 97 | + // Ïðÿìîé õîä ìåòîäà Ãàóññà |
| 98 | + private void GaussForwardStroke(int[] index) { |
| 99 | + for (int i = 0; i < size; ++i) { |
| 100 | + double r = FindR(i, index); |
| 101 | + for (int j = 0; j < size; ++j) |
| 102 | + a_matrix[i, j] /= r; |
| 103 | + b_vector[i] /= r; |
| 104 | + for (int k = i + 1; k < size; ++k) { |
| 105 | + double p = a_matrix[k, index[i]]; |
| 106 | + for (int j = i; j < size; ++j) |
| 107 | + a_matrix[k, index[j]] -= a_matrix[i, index[j]] * p; |
| 108 | + b_vector[k] -= b_vector[i] * p; |
| 109 | + a_matrix[k, index[i]] = 0.0; |
| 110 | + } |
| 111 | + } |
| 112 | + } |
| 113 | + |
| 114 | + // Îáðàòíûé õîä ìåòîäà Ãàóññà |
| 115 | + private void GaussBackwardStroke(int[] index) { |
| 116 | + for (int i = size - 1; i >= 0; --i) { |
| 117 | + double x_i = b_vector[i]; |
| 118 | + for (int j = i + 1; j < size; ++j) |
| 119 | + x_i -= x_vector[index[j]] * a_matrix[i, index[j]]; |
| 120 | + x_vector[index[i]] = x_i; |
| 121 | + } |
| 122 | + } |
| 123 | + |
| 124 | + // x - ðåøåíèå óðàâíåíèÿ, ïîëó÷åííîå ìåòîäîì Ãàóññà |
| 125 | + private void GaussDiscrepancy() { |
| 126 | + for (int i = 0; i < size; ++i) { |
| 127 | + double actual_b_i = 0.0; |
| 128 | + for (int j = 0; j < size; ++j) |
| 129 | + actual_b_i += initial_a_matrix[i, j] * x_vector[j]; |
| 130 | + } |
| 131 | + } |
| 132 | + } |
| 133 | +} |
0 commit comments