Метод Гаусса

Наиболее известным, универсальным и весьма эффективным точным методом решения систем линейных алгебраических уравнений является метод Гаусса или метод исключения.

В методе Гаусса, общая схема которого рассмотрена в разделе 8.5.1, матрица СЛАУ с помощью элементарных алгебраических операций преобразуется в верхнюю (нижнюю) треугольную матрицу, получающуюся в результате так называемого прямого хода. В процедуре обратного хода определяются неизвестные.

Рассмотрим подробнее одну из возможных реализаций этого метода на примере решения системы из n уравнений с n неизвестными.

Обозначим неизвестные через x ₁, x ₂, …, x _n и запишем СЛАУ в следующем виде:

(11.1)

Напомним, что элементарные преобразования, к которым относится перестановка двух уравнений системы, умножение уравнения системы на любое отличное от нуля число, а также сложение (вычитание) одного уравнения, умноженного на любое отличное от нуля число, с любым другим уравнением системы, переводят данную систему уравнений в эквивалентную. Таким образом, решение системы при элементарных преобразованиях не изменяется.

Первый шаг прямого хода метода Гаусса организуем так. Введем коэффициенты m ₂, m ₃, …, m_n, каждый из которых определяется следующим образом

Заметим, что обязательным условием данной операции является неравенство нулю элемента матрицы a ₁₁.

Начиная со второго уравнения, вычтем из каждого i-го уравнения первое, умноженное на соответствующий коэффициент m _i. Новые значения элементов матрицы системы и свободные члены будут равны

Индекс i в этих формулах обозначает номер уравнения, а индекс j номер столбца матрицы. Верхний индекс ⁽¹⁾ у элементов означает, что это преобразованные элементы, полученные на первом шаге прямого хода. После сделанного таким образом преобразования, во всех уравнениях, начиная со второго, коэффициенты при x ₁ будут равны нулю:

Преобразованная система уравнений (эквивалентная исходной) будет иметь вид

(11.2)

На этом первый шаг прямого хода метода Гаусса заканчивается.

Полученные уравнения, начиная со второго до n-го, представляют систему из n-1 уравнения с n-1 неизвестным. Повторяя описанные действия к этой системе на втором шаге прямого хода можно исключить x₂ из всех уравнений, начиная с третьего. На третьем шаге прямого хода исключается x₃ из всех уравнений, начиная с четвертого и т.д.

На k -м шаге прямого хода исключаем x _k из всех уравнений, начиная с (k +1)‑го с помощью множителей

, (11.3)

где верхний индекс (k -1) у элементов матрицы системы означает, что значения этих элементов рассчитаны на (k -1)-м шаге прямого хода метода Гаусса.

Также как и на первом шаге, обязательным условием данной операции является неравенство нулю элемента матрицы .

Повторяя действия, описанные выше для первого шага, вычтем из каждого уравнения, начиная с (k+1), k-е уравнение, умноженное на соответствующий коэффициент .

Новые значения элементов матрицы системы и свободные члены после k-го шага будут равны

(11.4)

На (n –1)-м шаге прямого хода произойдет исключение неизвестного
x_{n -1} из последнего уравнения. В итоге после (n –1)-го шага система примет треугольный вид:

(11.5)

На этом шаге прямой ход метода Гаусса заканчивается.

Описанный выше словесный алгоритм прямого хода метода Гаусса показан на рис. 9 в виде блок схемы и приведен в листинге 1 в виде фрагмента программы на Visual Basic (VB).

Рис.9. Блок схема алгоритма прямого хода метода Гаусса решения СЛАУ

Листинг 1. Прямой ход метода Гаусса без выбора главного элемента

For k = 1 To N - 1

For i = k + 1 To N

m = A(i, k) / A(k, k)

A(i, k) = 0

For j = k + 1 To N

A(i, j) = A(i, j) - m * A(k, j)

Next

B(i) = B(i) - m * B(k)

Next

Next

Для лучшего понимания алгоритма, выделим еще раз смысл переменных, использованных в блок схеме и фрагменте программы. Переменная n обозначает число уравнений в системе, индекс i обозначает номер уравнения, из которого исключается неизвестное, k – номер шага прямого хода и соответственно номер уравнения, которое вычитается из лежащих ниже n–k уравнений и одновременно номер исключаемого на данном шаге неизвестного, j – номер столбца. Двумерный массив A(1..n,1..n) и одномерный массив B(1..n) в программе на VB имеют тип, соответствующий числам с плавающей точкой (типы single или double в VB) и к моменту выполнения приведенного в листинге фрагмента должны быть инициализированы значениями матрицы системы и вектором правых частей системы.

Особенностями алгоритма, приведенного в виде блок схемы и программы, является то, что весь массив множителей не запоминается. В переменной m, имеющей тип single, или double, хранится только текущее значение множителя на k -м шаге для i -го уравнения. В профессиональных пакетах подобных программ все множители сохраняются обычно в самой матрице системы, где их размещают на месте нулевых элементов. Это позволяет рационально построить расчеты в случае, когда решение надо получить для нескольких правых частей системы.

Блок схема алгоритма, приведенная на рис. 9, дополнена, по сравнению со словесным алгоритмом, также специальной процедурой выбора главного элемента, которая требует отдельного пояснения. Относительно этой процедуры можно сказать следующее.

Во всех рассуждениях при построении алгоритма мы предполагали, что диагональные элементы матрицы системы, на которые производится деление при нахождении коэффициентов m, не равны нулю. Если это условие не выполняется, построенный нами алгоритм не позволит нам найти решение. Например, решение не будет найдено для систем, подобных приведенной ниже, так как элемент матрицы a ₁₁ равен нулю.

Данная система имеет тривиальное решение (x ₁=1, x ₂=1) и могла бы быть решена с помощью построенного алгоритма, если бы мы поменяли уравнения местами. На самом деле для того, чтобы алгоритм работал более эффективно, необходимо выполнение еще более жесткого условия, наложенного на значения диагональных элементов, а именно, необходимо, чтобы диагональные элементы матрицы были наибольшими по модулю, по сравнению с другими элементами матрицы. В противном случае, если диагональные элементы не равны нулю, но малы по сравнению с другими, значительно возрастают ошибки округления при решении системы. (Анализ показывает, что абсолютная ошибка округления в определении коэффициентов a_ij на каждом шаге прямого хода пропорциональна модулю коэффициентов m, вычисляемых по формуле (11.3).)

Для выполнения указанных условий в алгоритм необходимо ввести дополнительную процедуру, идея которой заключается в следующем. Будем просматривать в начале каждого шага прямого хода значения коэффициентов при исключаемой (k -й) переменной во всех уравнениях, лежащих ниже k -го уравнения, которое должно вычитаться. Найдя максимальный по модулю элемент матрицы (в k -м столбце), поменяем местами уравнение с максимальным элементом и k -е уравнение. После такой перестановки на диагонали матрицы будет оказываться каждый раз максимальный по модулю элемент столбца. Данная процедура называется выбором главного элемента.

Из теории матриц известно, что каждая перестановка уравнений не влияет на решение системы, но изменяет знак определителя матрицы. Поэтому, если мы хотим одновременно с решением системы вычислять определитель матрицы, который равен по модулю произведению диагональных элементов преобразованной к треугольному виду матрицы (уравнения (11.5)), необходимо предусмотреть в процедуре выбора главного элемента подсчет общего количества перестановок уравнений.

На листинге 2 приведен алгоритм прямого хода метода Гаусса на языке VB, дополненный процедурой выбора главного элемента.

Листинг 2. Прямой ход метода Гаусса с процедурой выбора главного элемента

For k = 1 To N – 1

‘ Начало процедуры выбора главного элемента

‘ Переменная l целого типа предназначена для запоминания номера уравнения,

‘ содержащего наибольший по модулю элемент матрицы в k-м столбце

l = k

For i = k + 1 To N

If Abs(A(i, k)) > Abs(A(l, k)) Then

l = i

End If

Next

‘ если в k-м столбце найден элемент, больший по модулю диагонального,

‘ то меняем местами уравнение с номером k и уравнение с номером l

If l <> k Then

For j = k To N

Z = A(k, j) ‘(Z- вспомогательная переменная, ‘используемая только в данном фрагменте для перестановки)

A(k, j) = A(l, j)

A(l, j) = Z

Next

Z = B(k)

B(k) = B(l)

B(l) = Z

End If

‘ окончание процедуры выбора главного элемента и продолжение прямого ‘хода

For i = k + 1 To N

m = A(i, k) / A(k, k)

A(i, k) = 0

For j = k + 1 To N

A(i, j) = A(i, j) - m * A(k, j)

Next

B(i) = B(i) - m * B(k)

Next

Next

Обратный ход метода Гаусса в нашей реализации, когда система уравнений приведена к верхнему треугольному виду (11.5), начинается с определения значения x_n из последнего уравнения в (11.5):

.

Подставляя найденное значение в предпоследнее уравнение в (11.5), можно найти x_{n -1}:

.

В свою очередь известные значения x_n и x_{n -1} можно подставить в уравнение с номером n –2 и найти x_{n -2} и т.д.

Общую формулу нахождения x _k на k-м шаге обратного хода можно записать следующим образом:

. (11.6)

Алгоритм обратного хода метода Гаусса, записанный на языке VB приведен на листинге 3.

Листинг 3. Обратный ход метода Гаусса

‘ Предполагается, что в массиве a(1..n,1..n) содержится верхняя треугольная матрица из (5)

‘Одномерный массив b(1..n) содержит правые части преобразованных уравнений (5)

‘Решение системы уравнений сохраняется в одномерном массиве X(1..n)

X(n) = B(n) / A(n, n)

For k = n - 1 To 1 Step -1

s = 0 ‘s- вспомогательная переменная, используемая только

‘ в данном фрагменте для накопления разности

For j = i + 1 To n

s = s - A(i, j) * X(j)

Next

X(i) = (B(i) + s) / A(i, i)

Next

Подсчитаем число арифметических операций, требуемых для выполнения описанного алгоритма. На первом шаге прямого хода требуется n ² умножений и делений, на втором шаге (n –1)² умножений и делений и т.д. Всего на n шагов прямого хода требуется

умножений и делений.

Для n шагов обратного хода требуется сделать n (n– 1)/2 умножений и делений. Суммируя число операций прямого и обратного хода, получим, что для решения системы (11.1) методом Гаусса без использования процедуры выбора главного элемента требуется n (n ²+3 n –1)/3 операций умножения и деления. Примерно столько же требуется операций сложения. Следовательно, для больших n можно считать, что метод Гаусса требует примерно арифметических действий и n ² ячеек памяти.

Если бы отсутствовали ошибки округления, то после выполнения алгоритма обратного хода было бы найдено точное решение системы. Наличие ошибок округления, от которых невозможно избавиться при решении задачи на ЭВМ вследствие конечной разрядности представления чисел, приводит к тому, что найденное решение всегда будет приближенным. (В методе Гаусса ошибки округления накапливаются в ходе выполнения алгоритма и, чем больше размерность системы, тем сильнее они могут повлиять на получаемое решение.)