namespace Gauss public class GaussSolutionNotFoun Exception public Gau

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
namespace Gauss {
public class GaussSolutionNotFound : Exception {
public GaussSolutionNotFound(string msg)
: base("Решение не может быть найдено: \r\n" + msg) {
}
}
public class LinearSystem {
private double[,] matrix_A2;
private double[,] matrix_A; // матрица A
private double[] mas_X; // вектор неизвестных x
private double[] mas_B2;
private double[] mas_B; // вектор b
private double[] mas_U; // вектор невязки U
private double eps; // порядок точности для сравнения вещественных чисел
private int size; // размерность задачи
private double[,] treg_Matrix;
private double det_A;
public LinearSystem(double[,] matrix_A2, double[] mas_B)
: this(matrix_A2, mas_B, 0.0001) {
}//наследуемый конструктор, если точность не задана
public LinearSystem(double[,] matrix_A, double[] mas_B, double eps) {
if (matrix_A == null || mas_B == null)
throw new ArgumentNullException("Один из параметров равен null.");
int b_length = mas_B.Length;
int a_length = matrix_A.Length;
if (a_length != b_length * b_length)
throw new ArgumentException("Количество строк и столбцов в матрице A должно совпадать с количеством элементров в векторе B.");
this.matrix_A2 = matrix_A; // запоминаем исходную матрицу для невязки
this.matrix_A = (double[,])matrix_A2.Clone(); //дублируем матрицу для
this.mas_B2 = mas_B;
this.mas_B = (double[])mas_B2.Clone();
this.mas_X = new double[b_length];
this.mas_U = new double[b_length];
this.size = b_length;
this.eps = eps;
this.treg_Matrix = new double[a_length, a_length];
GaussSolve();
}
public double[] XVector {
get {return mas_X;}
}
public double[] UVector {
get {return mas_U;}
}
public double[,] Treg_Matrix {
get{return treg_Matrix;}
}
public double Determ_A {
get{return det_A;}
}
// инициализация массива индексов столбцов
private int[] NumIndex() {
int[] index = new int[size];
for (int i = 0; i < index.Length; ++i)
index[i] = i;
return index;
}
// поиск главного элемента в матрице
public double FindBig(int row, int[] index) {
int max_i = row;
double max = matrix_A[row, index[max_i]];
double max_abs = Math.Abs(max);
//if(row < size - 1)
for (int index_now = row + 1; index_now < size; ++index_now) {
double now = matrix_A[row, index[index_now]];
double now_abs = Math.Abs(now);
if (now_abs > max_abs) {
max_i = index_now;
max = now;
max_abs = now_abs;
}
}
if (max_abs < eps) {
if (Math.Abs(mas_B[row]) > eps)
throw new GaussSolutionNotFound("Система уравнений несовместна.");
else
throw new GaussSolutionNotFound("Система уравнений имеет множество решений.");
}
// меняем местами
for (int i = 0; i < size; ++i) {
double tmp = matrix_A[i, index[row]];
matrix_A[i, index[row]] = matrix_A[i, index[max_i]];
matrix_A[i, index[max_i]] = tmp;
}
return max;
}
// Нахождение решения СЛАУ методом Гаусса
private void GaussSolve() {
int[] index = NumIndex();
Forward(index);
Backward(index);
Find_U();
}
// Прямой ход метода Гаусса (треугольный вид)
private void Forward(int[] index) {
for (int i = 0; i < size; ++i) {
double r = FindBig(i, index);
for (int j = 0; j < size; ++j)
matrix_A[i, j] /= r;
mas_B[i] /= r;
for (int k = i + 1; k < size; ++k) {
double p = matrix_A[k, index[i]];
for (int j = i; j < size; ++j)
matrix_A[k, index[j]] -= matrix_A[i, index[j]] * p;
mas_B[k] -= mas_B[i] * p;
//matrix_A[k, index[i]] = 0.0;
}
}
treg_Matrix = matrix_A;
}
// Обратный ход метода Гаусса
private void Backward(int[] index) {
for (int i = size - 1; i >= 0; --i) {
double x_i = mas_B[i];
for (int j = i + 1; j < size; ++j)
x_i -= mas_X[index[j]] * matrix_A[i, index[j]];
mas_X[index[i]] = x_i;
}
}
// Вычисление невязки решения
// U = b - x * A
private void Find_U() {
for (int i = 0; i < size; ++i) {
double actual_b_i = 0.0; // результат перемножения i-строки
// исходной матрицы на вектор x
for (int j = 0; j < size; ++j)
actual_b_i += matrix_A2[i, j] * mas_X[j];
// i-й элемент вектора невязки
mas_U[i] = mas_B2[i] - actual_b_i;
}
}
//Получаем минор
public double[,] get_minor(double[,] matrix, int row, int column) {
if (matrix.GetLength(0) != matrix.GetLength(1)) throw new Exception(" Матрица не квадратная");
double[,] buffer = new double[matrix.GetLength(0) - 1, matrix.GetLength(0) - 1];
for (int i = 0; i < matrix.GetLength(0); i++)
for (int j = 0; j < matrix.GetLength(1); j++)
{
if ((i != row) || (j != column))
{
if (i > row && j < column) buffer[i - 1, j] = matrix[i, j];
if (i < row && j > column) buffer[i, j - 1] = matrix[i, j];
if (i > row && j > column) buffer[i - 1, j - 1] = matrix[i, j];
if (i < row && j < column) buffer[i, j] = matrix[i, j];
}
}
return buffer;
}
//Рассчитываем определитель
public double Determinant(double[,] matrix) {
if (matrix.GetLength(0) != matrix.GetLength(1)) throw new Exception(" Число строк в матрице не совпадает с числом столбцов");
double det_A = 0;
int Rang = matrix.GetLength(0);
if (Rang == 1) det_A = matrix[0, 0];
if (Rang == 2) det_A = matrix[0, 0] * matrix[1, 1] - matrix[0, 1] * matrix[1, 0];
if (Rang > 2)
{
for (int j = 0; j < matrix.GetLength(1); j++)
{
det_A += Math.Pow(-1, 0 + j) * matrix[0, j] * Determinant(get_minor(matrix, 0, j));
}
}
return det_A;
}
}
}