Рассмотрим сейчас алгоритм Гаусса, позволяющий найти решение всех интересующих нас задач - вычислить определитель матрицы, решить m систем линейных уравнений, найти обратную матрицу. Построим вначале расширенную матрицу , состоящую из двух клеток:
Матрица , дополняющая матрицу , зависит от того, какую задачу предполагается решить. Если нужно вычислить только определитель матрицы , то расширенная матрица совпадает с исходной и матрица в этом случае отсутствует. Если нужно решить одну систему линейных уравнений, то матрица состоит из одного столбца - правых частей системы уравнений. Если нужно решить m систем уравнений, то матрица состоит из векторов, каждый из которых задает правые части своей системы уравнений. Если нужно найти обратную матрицу, то матрица задается единичной матрицей .
После того, как построена расширенная матрица, вся специфика конкретной задачи теряется - над расширенной матрицей выполняются одни и те же действия с параллельным вычислением определителя матрицы . В чем суть этих действий? Над матрицей последовательно выполняются элементарные преобразования - деление элементов строки на число, что изменяет величину определителя, и вычитание из одной строки матрицы другой строки, умноженной на некоторое число. Цель наших действий состоит в том, чтобы в расширенной матрице клетку преобразовать в единичную матрицу . Поскольку каждое элементарное действие можно рассматривать, как умножение слева на некоторую матрицу, совокупность преобразований, переводящая в , эквивалентна умножению слева на матрицу . Но это означает также, что эти преобразования переводят клетку в матрицу , что и дает решение исходных задач. Поскольку в результате преобразования переходит в единичную матрицу, определитель которой известен и равен 1, а для каждого преобразования известно, как меняется величина определителя, параллельно вычисляется и величина определителя исходной матрицы .
Рассмотрим на простом примере матричный вид элементарных операций. Пусть элементарная операция состоит в том, что к первой строке прибавляется вторая строка, умноженная на число . Это действие эквивалентно умножению почти единичной матрицы на исходную матрицу:
Матрица, задающая элементарную операцию, отличается от единичной матрицы тем, что у нее в первой строке на втором месте стоит число q, а не ноль. Если бы к первой строке прибавлялась не вторая строка, а строка с номером j, то число q стояло бы не на втором месте, а в позиции j. Если строка j прибавляется не к первой строке, а к строке с номером i, то число q появлялось бы в i-ой строке матрицы.
Рассмотрим теперь возможную реализацию алгоритма Гаусса:
public void Gauss(double[,] M) { det = 1; int n = M.GetLength(0); int m = M.GetLength(1); double d =0,r=0; for (int i = 0; i < n; i++) { //Приведение столбца i к единичному вектору d = M[i, i]; det *= d; //деление на диагональный элемент: M[i,i]теперь = 1; for (int k = 0; k < m; k++) M[i, k] /= d; //Элементарная операция: сложение строк for (int j=0; j<n; j++) { //К строке j прибавляется строка i, умноженная на r //В результате M[j,i]=0 if(j!=i) { r=-M[j,i]; for (int k = 0; k < m; k++) M[j, k] += r * M[i, k]; } } }
Аргументом метода является расширенная матрица . В результате работы метода матрица приобретает вид: . В зависимости от того, как задана матрица B, находится решение одной системы уравнений, нескольких систем или вычисляется значение обратной матрицы. Параллельно в переменной det формируется значение определителя матрицы A.
Алгоритм Гаусса в том виде, как он выше рассмотрен, не всегда обеспечивает получение результата. Действительно, пусть в матрице А элемент a[1,1] равен нулю. Тогда при выполнении элементарных операций в процессе преобразования матрицы А к единичной матрице Е возникнет ошибка уже на первом шаге при делении первой строки на элемент a[1,1]. Однако равенство нулю диагонального элемента вовсе не означает, что определитель матрицы равен нулю (если речь не идет о диагональной матрице) или что для нее не существует обратной матрицы.
Возможны различные модификации рассматриваемого алгоритма, исправляющие ситуацию.
Алгоритм с выбором первого ненулевого элемента
В случае, когда а[i, i] равно нулю, алгоритм ищет первую строку ниже i-й, в которой элемент a[i, j] не равен нулю. Эта строка добавляется к строке i, что гарантирует возможность деления на а[i, i].
Алгоритм с выбором главного элемента в столбце
Прежде чем приводить столбец к единичному виду, алгоритм ищет в столбце максимальный по модулю элемент и меняет местами строку i и строку j, в которой находится максимальный элемент. При обмене строк может измениться знак определителя матрицы.
Алгоритм с выбором главного элемента во всей матрице
На каждом шаге приведения очередного столбца к диагональному виду в еще не приведенной матрице отыскивается максимальный элемент и меняются местами не только строки, но и столбцы матрицы, ставя максимальный элемент в позицию а[i, i]. Этот прием гарантирует отсутствие переполнения при выполнении операции деления. Гарантируется также, что при умножениях не будет получено слишком большое число, поскольку деление на максимальный элемент с последующим умножением на один из элементов приводит к тому, что элементы преобразованной матрицы не увеличиваются по модулю. Однако ничто не дается даром. Выбор главного элемента, перестановка строк и столбцов, необходимость обратной перестановки в конце вычислений - все это усложняет алгоритм. Как правило, страдает и точность вычислений, особенно для плохо обусловленных матриц. Все модификации алгоритма стоит применять тогда, когда в основной схеме возникла исключительная ситуация, требующая корректировки алгоритма. Обработчик исключительной ситуации при делении на ноль, возникновении переполнения, потере значащих цифр может вызывать модифицированный вариант алгоритма в надежде получить решение, когда отказывается работать основная схема.
Замечу, что никакая модификация не может помочь найти обратную матрицу, если она не существует и определитель матрицы действительно равен нулю. В этом случае, например, все элементы в столбце, начиная с диагонального и ниже его, будут равны нулю. Это и будет означать, что определитель матрицы равен нулю, обратная матрица и решение системы уравнений не существует.