Г Л А В А   3

 

решение  систем  линейных
алгебраических  уравнений

 

 

 

§ 1. Основные понятия

 

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

Линейные системы. Запишем систему n линейных уравнений с n неизвестными:

 

                        (3.1)

 

Систему (3.1) можно записать более компактно используя матричную форму записи

,                                           (3.2)

 

где A – квадратная матрица, составленная из коэффициентов при неизвестных; x и с – вектор-столбец неизвестных и вектор-столбец правых частей (свободных членов), соответственно:

 

,     ,      .             (3.3)

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

Все эти случаи легко проиллюстрировать геометрически для системы двух уравнений

                                      (3.4)

 

Каждое уравнение описывает прямую на плоскости; координаты точки пересечения указанных прямых являются решением системы (3.4).

Возможны три случая взаимного расположения двух прямых на плоскости:

1)     прямые пересекаются, т.е. коэффициенты системы (3.4) не пропорциональны:

;                                             (3.5)

 

2)     прямые параллельны, т.е. коэффициенты системы (3.4) подчиняются условиям

;                                        (3.6)

 

3)     прямые совпадают, т.е. все коэффициенты пропорциональны

 

.                                      (3.7)

 

Запишем определитель системы (3.4)

 

.

 

Лишь при выполнении условия (3.5) определитель  не равен нулю и система имеет единственное решение (прямые пересекаются в одной точке). В случаях отсутствия решений (прямые параллельны) или при бесконечном множестве решений (прямые совпадают) выполняются соответственно соотношения (3.6) и (3.7), из которых получаем .

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

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

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

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

Важное достоинство итерационных методов состоит в том, что погрешности окончательных результатов не накапливаются, поскольку точность вычислений в каждой итерации определяется лишь результатами предыдущей итерации и практически не зависит от ранее выполненных вычислений.

 

 

§ 2. Прямые методы

 

Вводные замечания. Одним из способов решения систем линейных алгебраических уравнений является правило Крамера, согласно которому каждое неизвестное представляется в виде отношения определителей:

 

,

 

где  – определитель матрицы, получающейся из матрицы А заменой i-го столбца столбцом правых частей системы (3.1). Определители при этом предлагается вычислять по формулам, рассматриваемым в курсах линейной алгебры, например

 

.                          (3.8)

 

Здесь индексы a, b, …, w  пробегают все возможные  перестановок номеров ; k – число инверсий в данной перестановке.

Однако в качестве конкретного метода решения системы (3.1) данные формулы совершенно неприменимы, так как при подсчете каждого определителя по приведенной выше формуле надо вычислить  слагаемых, что нереально при весьма умеренных значениях n. Например, уже при  имеем . Если одно слагаемое вычисляется, скажем, за  с, что вполне допустимо для современных машин, то время расчета составит совершенно фантастическую цифру:

 

!

 

Фактически же в настоящее время с использованием подходящих методов решаются системы гораздо более высокого порядка  (до ).

Очевидно, что сложность системы (3.1) определяется структурой ее матрицы А. Существуют два случая, когда система имеет простые решения. Если А – диагональная матрица

 

,

 

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

Просто решается задача и в случае, когда матрица А является треугольной

.

 

В этом случае из последнего уравнения следует

 

,

и далее

,

 

для .

Оценим объем вычислений, связанный с решением системы с треугольной матрицей. Для того, чтобы вычислить  требуется одна операция, для вычисления   – три,  – пять и т.д. Как можно подсчитать, общее число операций при этом равно .

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

Метод Гаусса. Метод исключения Гаусса содержит два этапа: прямой ход – приведение исходной матрицы к треугольному виду путем последовательного исключения неизвестных из уравнений системы и обратный ход – решение системы с треугольной матрицей. Рассмотрим одну из возможных реализаций прямого хода.

Пусть . Тогда этот элемент называется ведущим или главным. Введем обозначение , .

Вычтем последовательно из второго, третьего, …, n-го уравнения системы (3.1) первое уравнение, умноженное соответственно на , , …, . Это позволит обратить в нуль  коэффициенты при  во всех уравнениях, кроме первого. В результате получим эквивалентную систему в виде

 

                        (3.9)

 

где ,  ;  .

Теперь исключим неизвестное  из уравнений, начиная с третьего. Пусть  – новый ведущий элемент; положим снова ,  и вычтем из третьего, четвертого, …, n-го уравнения второе уравнение, умноженное соответственно на , , …, . В результате получим

                      (3.10)

 

где ,  ;  .

Проделывая описанные действия  раз получим систему уравнений с треугольной матрицей

 

                    (3.11)

 

На этом завершается прямой ход метода Гаусса. Общее количество арифметических операций, которые необходимо совершить для приведения системы к треугольному виду описанным методом, оценивается величиной  ~. Это вполне приемлемая величина (при n ~  и быстродействии ЭВМ порядка  операций в секунду требуемое для расчета время порядка одного часа).

Ø Замечание. В процессе исключения неизвестных приходится выполнять операции деления на коэффициенты ,  и т.д. Поэтому они должны быть отличны от нуля. В противном случае необходимо соответственным образом переставить уравнения системы. Перестановка уравнений должна быть предусмотрена в вычислительном алгоритме при его реализации на компьютере.<

Обратный ход начинается с решения последнего уравнения системы (3.11):

,

 

затем из предпоследнего уравнения находится :

 

,

 

аналогично определяются и все остальные неизвестные.

Ниже приведено описание вычислительного алгоритма метода Гаусса.

 

1.     Ввод n, ,

2.     Для i от 1 до  :

2.1.Если , тогда Перестановка уравнений

2.2.Для  k  от  до n :

2.2.1.  

2.2.2.   Для  j  от  до  n :  

2.2.3.  

3.     Для i  от  n  до 1 с шагом –1 :

3.1.

3.2.Для j от  до n : 

3.3.

4.     Вывод

 

В приведенном алгоритме первый цикл с параметром i реализует прямой ход, а второй – обратный ход метода Гаусса.

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

Алгоритм выбора наибольшего элемента по столбцу достаточно прост и состоит из двух этапов: отыскание номера наибольшего элемента и перестановки элементов двух строк. Описание алгоритма приведено ниже.

 

1.    

2.     Для m от  до n :

2.1.Если , тогда

3.     Если , тогда

3.1.Для  j от i  до n :

3.2.

 

Здесь k – номер наибольшего по абсолютной величине элемента матрицы в столбце с номером i. Этот алгоритм встраивается в приведенный выше алгоритм метода Гаусса вместо условной конструкции, выполняющей перестановку уравнений в случае равенства нулю элемента  (шаг 2.1.).

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

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

Метод прогонки специально создан для решения систем линейных алгебраических уравнений с трехдиагональной матрицей, т.е. для систем вида

                (3.12)

Системы такого вида часто возникают при решении различных задач вычислительной математики (например, при численном решении дифференциальных уравнений и интерполировании сплайнами).

Выразим из первого уравнения :

 

.                    (3.13)

 

Выражение для  подставим во второе уравнение системы (3.12):

 

,

 

Выражая из последнего соотношения  получим:

 

.             (3.14)

 

Выражение для  подставим в третье уравнение системы (3.12) и так далее. На i-м шаге описанного процесса () i-е уравнение системы приводится к такому же виду:

 

,    где

(3.15)

, 

 

На последнем n-м шаге подстановка получаем n-е уравнение в виде . Отсюда можно выразить неизвестное :

 

.                                  (3.16)

 

Значения остальных неизвестных вычисляются в процессе обратного хода по формулам (3.15).

Таким образом, прямой ход состоит в вычислении прогоночных  коэффициентов  и  по формулам:

 

,   ,   ,        (3.17)

,   ,                   (3.17а)

 

Обратный ход метода прогонки дает значения неизвестных по формулам:

,   .                      (3.18)

 

Отметим, что метод прогонки относится к классу экономичных методов. Экономичными называются методы, для которых число требуемых арифметических операций пропорционально числу неизвестных. Для реализации вычислений по описанному алгоритму требуется примерно 8n арифметических операций (тогда как в методе исключения Гаусса эта величина ~). Экономичность в данном случае достигается за счет того, что при реализации метода исключения, приводящего к формулам (3.17) и (3.17а), не выполнялись операции над нулевыми элементами исходной матрицы.

Вычисление определителей. Как было отмечено ранее, непосредственное вычисление определителя матрицы по формулам линейной алгебры требует достаточно большого количества операций, что при больших значениях n (порядка матрицы) приводит к нереальным затратам машинного времени. Вместе с тем известно, что легко вычисляется определитель треугольной матрицы: он равен произведению ее диагональных элементов.

Для приведения матрицы к треугольному виду может быть использован метод исключения, т.е. прямой ход метода Гаусса. В процессе исключения неизвестных величина определителя не меняется, а знак определителя меняется на противоположный при каждой перестановке его столбцов или строк. Следовательно, значение определителя матрицы А, приведенной к треугольному виду,  вычисляется по формуле

 

.

 

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

Нахождение обратных матриц. Матрица  называется обратной по отношению к матрице А, если их произведение равно единичной матрице: . Согласно формулам линейной алгебры каждый элемент  обратной матрицы  равен отношению алгебраического дополнения  элемента  исходной матрицы к значению её определителя. Алгебраическое дополнение , в свою очередь, вычисляется как определитель -го порядка, образованный из определителя матрицы А вычеркиванием i-й строки и j-го столбца и взятый со знаком плюс, если сумма  является четной, и со знаком минус в противном случае. Таким образом, задача нахождения обратной матрицы сводится к вычислению определителей, что без использования специальных методов может потребовать огромных затрат времени.

В тоже время, нахождение обратной матрицы можно свести к решению последовательности систем линейный алгебраических уравнений. Запишем определение обратной матрицы  в виде

 

,

 

откуда следует, что

 

,   ,                          (3.19)

 

где  и  j-е столбцы матриц   и Е, соответственно. Таким образом, для нахождения j-го столбца обратной матрицы нужно решить систему уравнений (3.19). Решив n таких систем для  мы найдем все столбцы  и, следовательно, саму обратную матрицу.

Поскольку при разных j матрица А системы (3.19) не меняется, исключение неизвестных при использовании метода Гаусса (прямой ход) проводится только один раз, причем сразу для всех правых частей – столбцов . Затем для каждой из систем (3.19) выполняется обратный ход с соответствующей правой частью.

 

 

§ 3. Итерационные методы

 

Метод простых итераций. Мы переходим к обсуждению итерационных методов (т.е. методов последовательных приближений) решения системы линейных уравнений:

 

.                                           (3.20)

 

Как и в предыдущем параграфе, будем считать, что решение (3.20) существует и единственно.

Различные варианты метода простых итераций связаны с переходом от системы (3.20) к эквивалентной системе

 

.                                          (3.21)

 

Операция приведения исходной системы (3.20) к виду (3.21) не является очевидной и будет подробно рассмотрена в дальнейшем.

Опираясь на выражение (3.21) можно построить простой итерационный процесс следующим образом:

 

,      – задано.                 (3.22)

 

Здесь  – номер приближения.

Если оказывается, что с увеличением числа итераций приближенное решение стремится к точному:

 

,

 

то итерационный процесс называют сходящимся.

Доказано (см., например, [8]), что для сходимости итераций (3.22) к точному решению системы независимо от выбора начального приближения  достаточно, чтобы  в какой-либо норме[1] выполнялось условие

 

.                                       (3.23)

 

При практической реализации метода простых итераций достижение требуемой точности оценивается, как правило, путем  сравнения двух последовательных приближений  и , и при малом отличии  от  процесс прекращается.

Для определения малости отличия х на двух последовательных итерациях можно использовать разные подходы. Например, если задана допустимая погрешность , то критерием окончания итерационного процесса можно считать выполнение одного из двух неравенств:

 

,                (3.24)

 

.                                 (3.25)

 

Здесь в первом случае отличие векторов  и  “на величину e” понимается в смысле малости модуля их разности, а во втором случае в смысле малости разностей всех соответствующих компонент векторов.

Метод итераций Якоби. Самый простой способ приведения системы (3.20) к виду (3.21) состоит в записи в явном виде неизвестного  из первого уравнения,  из второго и т.д. Метод итераций в такой реализации называется методом Якоби.

Проиллюстрируем сначала этот метод на примере системы трех уравнений:

                            (3.26)

 

Предположим, что диагональные элементы ,  и  отличны от нуля. Выразим неизвестные ,  и  соответственно из первого, второго и третьего уравнений системы (3.26):

 

                            (3.27)

 

Тем самым, мы получили систему уравнений, записанную в виде (3.21), причем матрица В и вектор g имеют здесь следующий вид:

 

,   .           (3.28)

 

С учетом выражений (3.27) построим итерационный процесс метода Якоби, в котором приближение с номером  вычисляется на основе приближения с номером k следующим образом:

 

                        (3.29)

 

В качестве вектора начальных приближений можно взять нулевой вектор или вектор g, т.е.

 

  или  .                                       (3.30)

 

Итерационный процесс (3.29) продолжается до тех пор, пока не выполнится одно из условий (3.24), (3.25).

Рассмотренный метод легко обобщается на систему из n уравнений. Элементы матрицы В и вектора g будут определятся в этом случае следующим образом

 

    ,    .       (3.31)

 

А система (3.29) принимает вид

 

,    .                     (3.32)

 

Для сходимости метода Якоби достаточно, чтобы для исходной матрицы А имело бы место диагональное преобладание, т.е. чтобы коэффициенты исходного уравнения удовлетворяли условиям

 

  для всех i.

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

Метод Зейделя.  Этот метод отличается от метода Якоби только тем, что при вычислении -го приближения i-й компоненты вектора х используются уже вычисленные -е приближения предыдущих (1-й, 2-й, …, -й) компонент вектора х.

Для системы (3.27) итерационный процесс Зейделя будет иметь следующий вид:

 

                   (3.33)

 

Если ввести нижнюю и верхнюю треугольные матрицы

 

   и  ,

 

то система (3.33) может быть компактно записана в матричном виде

 

.                       (3.34)

 

Здесь g – вектор, определенный выше выражением (3.31)

Итерационный процесс (3.26) непосредственно обобщается на систему n уравнений, при этом элементы матриц  и  определяются следующим образом

 

,    .

.                                      (3.35)

 

При программной реализации метода Зейделя его алгоритм можно организовать таким образом, что необходимость в введении двух треугольных матриц отпадает, и расчет ведется с одной матрицей В, элементы которой вычисляются по формуле (3.31). Опишем этот алгоритм:

 

1.     Ввод n, , , e

2.     Для i от 1 до n :

2.1.Для j от 1 до n :

2.1.1.   Если , тогда , иначе

2.2.;

3.     Для i от 1 до n : 

4.     Для i от 1 до n :

4.1.

4.2.Для  j от 1 до n : 

4.3.

5.     Если , тогда перейти к шагу 3

6.     Вывод

 

Поясним смысл некоторых переменных и шагов алгоритма. На первом шаге организуется ввод элементов матрицы А, вектора с и заданной точности e. На втором шаге вычисляются элементы матрицы В, вектора g, а также задается вектор начальных приближений x, в качестве которого берется вектор g. На третьем шаге запоминаются текущие значения элементов вектора x в векторе y, его смысл – хранить значения неизвестных до выполнения очередной итерации с тем, чтобы можно было сравнивать их значения на -й и k-й итерациях. На четвертом шаге совершается одна итерации в соответствии с формулой (3.34). На пятом шаге происходит сравнение значение элементов вектора x на двух последовательных итерациях и, если вектора отличаются больше чем на величину e (вычисления модуля разности векторов выполняется по формуле (3.24)), то шаги, начиная с третьего, повторяются еще раз.

Доказано, что метод Зейделя гарантировано сходится, если:

       для матрицы А выполнено условие диагонального преобладания

или

       матрица А является симметричной и положительно определенной, т.е. такой что  и .

Системы уравнений с симметричными и положительно определенными матрицами называются нормальными системами. Для приведения любой системы к нормальному виду, достаточно выполнить простые преобразования, а именно умножить левую и правую части системы слева на транспонированную матрицу[2] :

 

.

 

Если ввести новую матрицу  и вектор , то получившаяся эквивалентная система

 

 

будет обладать необходимыми свойствами нормальной системы.

Для систем с нормальной матрицей метод Зейделя сходится примерно в два раза быстрее метода Якоби. Однако в общем случае возможны ситуации, когда метод Якоби сходится, а метод Зейделя или сходится медленнее, или не сходится вовсе. Возможны также и противоположные ситуации.

 

 

Дополнение к § 3.

 

Нормы векторов и матриц.

Нормой вектора X называется число, обозначаемое  и удовлетворяющее условиям:

1)     ;

2)     ;

3)     ;

Наиболее употребительными являются следующие три нормы:

,

 – евклидова норма,

 – равномерная норма.

Векторное пространство с введенной в нем нормой называют нормированным. Одновременно оно является метрическим, так как норма определяет метрику – расстояние между точками пространства:

 

.

 

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

 

.

 

Введение понятия нормы вектора позволяет, в частности,  определять количественные характеристики погрешности приближенных значений векторов.

Пусть, например,   – точное решение системы линейных алгебраических уравнений, а  – приближенное решение. Тогда абсолютную и относительную погрешности приближенного решения можно определить следующим образом:

 

   и  .

 

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

 

.

 

При этом в условии (3.24) используется евклидова норма, а в условии (3.25) – равномерная.

Нормой матрицы А называется величина . Введенная таким образом норма обладает свойствами, аналогичными свойствам нормы вектора:

1)     , причем  тогда и только тогда, когда ;

2)      для любой матрицы А и любого числа a ;

3)      для любых матриц А и В;

4)     .

Каждой из векторных норм соответствует своя подчиненная норма матрицы:

,   ,   .

 

 

Вопросы для самоконтроля

1.     В чем разница прямых и итерационных методов решения систем линейных алгебраических уравнений?

2.     В чем заключается прямой и обратный ход метода исключения Гаусса?

3.     Почему в прямом ходе метода Гаусса целесообразно так переставить уравнения, чтобы ведущий элемент имел наибольшее значение?

4.     Для решения каких систем применяется метод прогонки?

5.     В чем отличие итерационного процесса метода Зейделя от аналогичного процесса метода Якоби?

6.     Каковы достаточные условия сходимости методов Якоби и Зейделя?

 



[1] Подробнее о нормах векторов и матриц см. дополнение к параграфу 3.

[2] Транспонированной называется матрица, в которой строки и столбцы поменяны местами.