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;
}
}
}