К прямымметодам относятся такие методы, которые в предположении, что вычисления ведутся без округлений, позволяют получить точные значения неизвестных. Они просты, универсальны и используются для широкого класса систем. Однако они не применимы к системам больших порядков (n < 200) и к плохо обусловленным системам из-за возникновения больших погрешностей. К ним можно отнести: правило Крамера, методы обратных матриц, Гаусса, прогонки, квадратного корня и др.
К приближенным относятся методы, которые даже в предположении, что вычисления ведутся без округлений, позволяют получить решение системы лишь с заданной точностью. Это итерационные методы, т. е. методы последовательных приближений. К ним относятся методы простой итерации, Зейделя.
Правило Крамера. Рассмотрим систему (2.1). Как отмечалось в подразд. 2.1, если определитель этой системы не равен нулю, то будет иметь место единственное решение, - это необходимое и достаточное условие. Тогда по правилу Крамера
, (2.4)
где Dk – определитель, получающийся из D при замене элементов a1k, a2k, ..., ank k-го столбца (соответствующими) свободными членами b1, b2, ..., bn из (2.1):
,
где Аik – алгебраическое дополнение элемента aik в определителе D.
Стоит существенная проблема вычисления определителей высоких порядков.
Метод обратных матриц. Дана система . Умножим левую и правую части этого выражения на А–1:
; .
При реализации данного метода стоит проблема нахождения обратной матрицы А–1, выбора экономичной схемы ее получения и достижения приемлемой точности.
Метод Гаусса. Этот метод является наиболее распространенным методом решения СЛАУ. В его основе лежит идея последовательного исключения неизвестных, в основном приводящая исходную систему к треугольному виду, в котором все коэффициенты ниже главной диагонали равны нулю.
Существуют различные вычислительные схемы, реализующие этот метод. Наибольшее распространение получили схемы с выбором главного элемента либо по строке, либо по столбцу, либо по всей матрице. С точки зрения простоты реализации, хотя и с потерей точности, перед применением этих схем целесообразно использовать так называемую схему единственного деления. Рассмотрим ее суть.
Посредством первого уравнения системы (2.1) исключается х1 из последующих уравнений. Далее посредством второго уравнения исключается х2 из последующих уравнений и т. д. Этот процесс называется прямым ходом Гаусса. Исключение неизвестных повторяется до тех пор, пока в левой части последнего n-го уравнения не останется одно неизвестное хn:
a¢nnxn = b¢, (2.5)
где a¢nn и b¢ – коэффициенты, полученные в результате линейных (эквивалентных) преобразований.
Прямой ход реализуется по формулам
(2.6)
где m – номер уравнения, из которого исключается xk; i – номер столбца исходной матрицы; k – номер неизвестного, которое исключается из оставшихся (n – k) уравнений и обозначает номер уравнения, с помощью которого исключается xk; akk – главный (ведущий) элемент матрицы.
Во время счета необходимо следить, чтобы akk ¹ 0. В противном случае прибегают к перестановке строк матрицы.
Обратный ход метода Гаусса состоит в последовательном вычислении xn, xn–1, ..., x1, начиная с (2.5) по алгоритму:
xn = b¢ / a¢nn; . (2.7)
Точность полученного решения оценивается посредством невязки (2.3). В векторе невязки (r1, r2, ..., rn)Т отыскивается максимальный элемент и сравнивается с заданной точностью e. Приемлемым решение будет, если rmax < e. В противном случае следует применить схему уточнения решения.
Полученные методом Гаусса приближенные значения корней можноуточнить. Пусть для системы найдено приближенное решение , не удовлетворяющее по «невязке». Положим тогда. Для получения поправки d = (d1, d2, ..., dn)Т корня следует рассмотреть новую систему
или ,
где – невязка для исходной системы.
Таким образом, решив линейную систему с прежней матрицей А и новым свободным членом = (e1, e2, ..., en)Т, получим поправки (d1, d2, ..., dn).
Пример решения СЛАУ по методу Гаусса (с точностью до трех знаков). Нужно уточнить корни до 10–4:
В результате = 4,67; = 7,62; = 9,05. Невязки равны = −0,02; = 0; = −0,01.
Модифицированный метод Гаусса. В данном случае, помимо соблюдения требования akk ¹ 0, при реализации формул (2.6) накладываются дополнительные требования: ведущий (главный) элемент в текущем столбце в процессе преобразований исходной матрицы должен иметь максимальное по модулю значение, что достигается перестановкой строк матрицы.
Пример. В качестве иллюстрации преимущества модифицированного метода Гаусса рассмотрим систему третьего порядка:
(2.8)
Прямой ход метода Гаусса. Исключаем х1 из второго и третьего уравнений (2.8). Для этого первое уравнение умножаем на 0,3 и складываем со вторым, а затем умножаем первое уравнение на (–0,5) и складываем с третьим. В результате получаем
(2.8а)
Замена второго уравнения третьим не производится, так как вычисления выполняются в рамках точной арифметики.
Умножив второе уравнение на 25 и сложив с третьим, получим
(2.8б)
Обратный ход метода Гаусса. Выполняем вычисления начиная с последнего уравнения в полученной системе:
Подставляя полученное решение [0; –1; 1] в исходную систему (2.8), убеждаемся в его истинности.
Теперь изменим коэффициенты системы таким образом, чтобы сохранить прежнее решение, но при вычислении будем использовать округления в рамках арифметики с плавающей точкой, сохраняя пять разрядов. Изменению вычислений будет соответствовать следующая система:
(2.8в)
Прямой ход метода для системы (2.8в) повторим по аналогичной технологии с исходной системой (2.8):
(2.8г)
После исключения х2 третье уравнение примет вид (остальные – без изменения):
15005х3 = 15004. (2.8д)
Выполнив обратный ход, получим
Очевидно, что полученные решения [0; –1; 1] и [–0,35; –1,4; 0,99993] различны. Причиной этого является малая величина ведущего элемента во втором уравнении преобразования (2.8г). Чтобы это исключить, переставим в (2.8г) вторую и третью строки:
º (2.8е)
Система (2.8е) после исключения х2 из третьего уравнения примет следующий вид:
6,002 х3 = 6,002. (2.8ж)
В данном случае, выполнив обратный ход:
получим решение системы (2.8в) [0; –1; 1], которое в точности совпадает с решением исходной системы.
Решая систему (2.8г), мы использовали модифицированный метод Гаусса, при выполнении которого на диагонали должен был находиться максимальный в текущем столбце элемент.
Рассмотрим блок-схему модифицированного метода Гаусса (рис. 2.1). Проведем анализ предложенной схемы на примере системы n = 3 (e = 0,001):
(2.9)
; . (2.9а)
Блок 1. Ввод исходных данных: n – порядок системы, A – матрица коэффициентов при неизвестных, – вектор свободных членов.
Блок 2. Цикл I прямого хода (для k, изменяющегося от единицы до предпоследнего значения, т. е. до n – 1) обеспечивает исключение из главной диагонали матрицы А элемента akk= 0 благодаря поиску максимального элемента akk в текущем столбце, осуществляемому в блоках 3 - 6 с помощью цикла II.
Далее с помощью цикла III в блоках 7 - 13 выполняется перестановка текущей строки и строки с максимальным элементом в k-м столбце (ее номер р).
Затем реализуются расчеты по формулам (2.6) прямого хода Гаусса в блоках циклов IV и V.
Проведем поблочный анализ в среде рассмотренных циклов I - V на примере (2.9).
Рис. 2.1
Блок 3 p = k = 1;
Вход в цикл II
Блок 4 m = k + 1 = 2; до n = 3;
Блок 5 a11 = 2 < a21 = 4; из (2.9а);
Блок 6 p = 2;
Блок 4 m = 2 + 1 = 3;
Блок 5 a21 = 4 < a31 = 6; из (2.9а);
Блок 6 p = 3.
Выход из цикла II и вход в цикл III, блоки 7 - 10 выполняют перестановку строк матрицы А поэлементно:
Блок 7 j = 1 (j от 1 до 3);
Блок 8 r = a11 = 2; из (2.9а);
Блок 9 a11 = a31 = 6;
Блок 10 a31 = r;
Блок 7 j = 2;
Блок 8 r = a12 = 1;
Блок 9 a12 = a32 = 5;
Блок 10 a32 = r = 1;
Блок 7 j = 3 и по аналогии r = a13; a13 = a33; a33 = r = −1.
Выход из цикла III и вход в блок 11. Блоки 11 - 13 выполняют аналогичную перестановку значений свободных членов:
r = b1 = 1; b1 = b3 = 14; b3 = r = 1.
Вход в цикл IV с измененной системой:
; . (2.9б)
В цикле IV выполняем пересчет значения b2 вектора :
m = k + 1 = 1 + 1 = 2; до n = 3;
c = amk / akk = a21 / a11 = 4/6; из (2.9б);
b2 = b2 – cb1 = 6 – 4/6 × 14 = −20/6; из (2.9б).
Вход во вложенный цикл V для пересчета второй строки:
i = 1 (i от 1 до 3); a21 = a21 – с × a11 = 4 – 4/6 × 6 = 0;
i = 2; a22 = a22 – с × a12 = 6 – 4/6 × 5 = 16/6;
i = 3; a23 = a23 – с × a13 = 2 – 4/6 × 8 = −20/6.
Выход из цикла V и возврат в цикл IV:
m = 3; c = a31 / a11 = 2/6.
Вход в блок 16:
b3 = b3 – c × b1 = 1 – 2/6 × 14 = −22/6.
Вход во вложенный цикл V (блок 17):
i = 1 (i от 1 до 3); a31 = a31 – с × a11 = 2 – 2/6 × 6 = 0;
i = 2; a32 = a32 – с × a12 = 1 – 2/6 × 5 = −4/6;
i = 3; a33 = a33 – с × a13 = −1 – 2/6 × 8 = −22/6.
Выход из циклов V и IV с преобразованной системой:
; ; (2.9в)
и возврат по линии А в цикл I:
k = 2; p = k = 2; m = k+1 = 3; вход в блок 5;
| a22 | < | a32 | = | 16/6 | > | 4/6 | из (2.9в).
Выход из цикла II и вход в цикл III:
j = 2 (j от 2 до 3);
r = akj = a22 = 16/6; a22 = a22 ; a22 = r = 16/6; из (2.9в);
j = 3;
r = a23 = −20/6; a23 = a23 ; a23 = r = −20/6; из (2.9в).
В данном случае на диагонали оказался максимальный элемент, поэтому перестановка 2-й и 3-й строк не выполняется.
Выход из цикла III и выполнение блоков 11 – 13:
r = b2; b2 = b2; b2 = r = −20/6.
Свободный член b2 остается на своем месте.
Вход в цикл IV:
m = k + 1 = 2 + 1 = 3;
c = amk / akk = a32 / a22 = (–4/6) / (16/6); из (2.9в);
b3 = b3 – c b2 = −22/6 – (–1/4) × (–20/6) = −27/6; из (2.9в).
Вход во вложенный цикл V:
i = 2 (i от 2 до 3); a32 = a32 – с × a22 = −4/6 – (–1/4) × 16/6 = 0;
i = 3; a33 = a33 – с × a23 = −22/6 – (–1/4) × (–20/6) = −27/6.
Выход из циклов V и IV и завершение (выход) цикла I.
Выполнение обратного хода метода Гаусса:
В блоках 19 - 24 реализуются формулы (2.7). В блоке 19 из последнего уравнения (2.7) находится значение xn (n = 3):
x3 = bn / ann = b3 / a33 = (–27/6) / (–27/6) = 1.
Вход в цикл VI (блок 20), в котором значение переменной цикла k изменяется от n – 1 до 1 с шагом (–1):
В блоке 25 (цикл опущен) выполняется вывод на экран полученного решения СЛАУ – вектора , т. е. xi, i = 1, ..., n. В данном случае (1; 0; 1).
Метод прогонки. Данный метод является модификацией метода Гаусса для частного случая разреженных систем – систем с матрицей трехдиагонального типа (краевая задача дифференциальных уравнений).
Метод реализуется в два этапа – прямой и обратный ходы.
Прямой ход. Каждое неизвестное xi выражается через xi+1:
xi= Ai× xi+1+ Bi для i = 1, 2, ..., n – 1, (2.12)
посредством прогоночных коэффициентов Ai и Bi. Определим алгоритм их вычисления.
Из первого уравнения системы (2.11) находим x1:
.
Из уравнения (2.12) при i = 1: x1= A1× x2+ B1, следовательно,
. (2.13)
Из второго уравнения системы (2.11) определяем x2 через x3, подставляя найденное значение x1:
а2( A1x2 + B1) + b2 x2 + c2 x3 = d2,
откуда
. (2.13а)
Согласно (2.12) при i = 2: x2= A2× x3+ B2, следовательно,
, где е2= а2× А1 + b2.
Ориентируясь на соотношения индексов при коэффициентах (2.12) и (2.12а), можно получить эти соотношения для общего случая:
, где еi= аi× Аi–1 + bi (i = 2, 3, ..., n – 1). (2.14)
Обратный ход.Из последнего уравнения системы (2.11) с использованием (2.12) при i = n – 1:
. (2.15)
Далее посредством (2.12) и прогоночных коэффициентов (2.13), (2.14) последовательно вычисляем xn– 1, xn– 2, ..., x1.
При реализации метода прогонки нужно учитывать, что при выполнении условия
| bi | ³ | ai | + | ci | (2.16)
или хотя бы для одного bi имеет место строгое неравенство (2.16), деление на ноль исключается и система имеет единственное решение.
Заметим, что условие (2.16) является достаточным, но не необходимым. В ряде случаев для хорошо обусловленных систем (2.11) метод прогонки может быть устойчивым и при несоблюдении условия (2.16).
Схема алгоритма метода прогонки может иметь вид, представленный на рис. 2.2.
Рис. 2.2
Метод квадратного корня. Данный метод используется для решения линейной системы
, (2.17)
у которой матрица А симметрическая, т. е. АТ = А, aij = aji (i = j = 1, ..., n).
Решение системы (2.17) осуществляется в два этапа.
Прямой ход. Преобразование матрицы А и представление ее в виде произведения двух взаимно транспонированных треугольных матриц:
А = SТ × S, (2.18)
где
Перемножив SТ и S и приравняв матрице А, получим следующие формулы для определения sij:
(2.19)
После нахождения матрицы S систему (2.17) заменяем двумя ей эквивалентными системами с треугольными матрицами (2.18):
. (2.20)
Обратный ход. Записываем системы (2.20) в развернутом виде:
(2.21)
(2.22)
Используя (2.21) и (2.22), последовательно находим:
(2.23)
(2.24)
Метод квадратного корня дает большой выигрыш во времени по сравнению с рассмотренными ранее прямыми методами, так как, во-первых, существенно уменьшает число умножений и делений (почти в два раза), во-вторых, позволяет накапливать сумму произведений без записи промежуточных результатов. Числовой пример ручного счета рассмотрен в [1].
Машинная реализация метода квадратного корня предусматривает его следующую трактовку.
Исходная матрица А системы (2.17) представляется в виде произведения трех матриц:
A = S Т D S,
где ST – транспонированная нижняя треугольная матрица; D – диагональная матрица с элементами dii = ±1; S – верхняя треугольная (sik = 0, если i > k , причем sii > 0).
Выполнение условия sii > 0 необходимо для полной определенности разложения. Это и определяет необходимость введения диагональной матрицы D.
Рассмотрим алгоритм разложения матрицы А с использованием матрицы D на примере матрицы второго порядка.
Пусть А – действительная симметричная матрица:
.
Будем искать S и D в виде
, , dii = ±1.
Тогда
.
Из условия равенства A = SТ D S получим три уравнения:
Из первого уравнения находим
d11 = sign a11; s11 = .
Далее, если а11 ¹ 0, то s12 = а12 / (s11 d11), и, наконец
,
т. е. d22 = sign (a22); s22 = .
Здесь, как и для общего случая, матрицу S можно по аналогии с числами трактовать как корень квадратный из матрицы А, отсюда и название метода.
Итак, если SТDS известно, то решение исходной системы сводится к последовательному решению систем:
. (2.24а)
Нахождение элементов матрицы S (извлечение корня из А) выполняем по рекуррентным формулам, что исключает использование комплексных чисел:
dk = sign;
skk = ; (2.25)
skj = ;
k = 1, 2, ..., n; j = k + 1, k + 2, ..., n.
В (2.25) сначала полагаем k = 1 и последовательно вычисляем
d1 = sign (a11); s11 =
и все элементы первой строки матрицы S (s1j, j > 1), затем полагаем k = 2, вычисляем s22 и вторую строку матрицы s1j для j > 2 и т. д.
Решение систем (2.24а) ввиду треугольности матрицы S осуществляется по формулам, аналогичным обратному ходу метода Гаусса:
Метод квадратного корня почти вдвое эффективнее метода Гаусса, так как использует симметричность матрицы.
Схема алгоритма метода квадратного корня представлена на рис. 2.3. Значение функции sign(x) равно +1 для всех х > 0 и –1 для всех x < 0. Алгоритм реализован в [6].
Проиллюстрируем метод квадратного корня, решив систему трех уравнений:
, .
Нетрудно проверить, что матрица А есть произведение двух треугольных матриц (здесь dii = 1):
, (2.4)
,
. Умножим левую и правую части этого выражения на А–1:
; 
.
(2.6)
. (2.7)
, не удовлетворяющее по «невязке». Положим тогда
. Для получения поправки d = (d1, d2, ..., dn)Т корня 
или 
,
– невязка для исходной системы.
= (e1, e2, ..., en)Т, получим поправки (d1, d2, ..., dn).
= 4,67; 
= 7,62; 
= 9,05. Невязки равны 
= −0,02; 
= 0; 
= −0,01.
= −0,0039; 
= −0,0011; 
= −0,0025. Следовательно, х1 = 4,6661; х2 = 7,6189; х3 = 9,0475.
(2.8)
(2.8а)
(2.8б)
(2.8в)
(2.8г)
º 
(2.8е)
(2.9)
; 
. (2.9а)
– вектор свободных членов.
; 
. (2.9б)
; 
; (2.9в)
, т. е. xi, i = 1, ..., n. В данном случае (1; 0; 1).
; a1 = cn = 0, (2.10)
или в развернутом виде
.
. (2.13)
. (2.13а)
, где е2= а2× А1 + b2.
, где еi = аi × Аi–1 + bi (i = 2, 3, ..., n – 1). (2.14)
. (2.15)
, (2.17)
(2.19)
. (2.20)
(2.21)
(2.22)
(2.23)
(2.24)
.
, 
, dii = ±1.
.
.
,
); s22 = 
.
сводится к последовательному решению систем:
. (2.24а)
dk = sign
;
; (2.25)
;
, 
.
.
.
.
получим систему 
:
, откуда y1 = 3; y2 = 2; y3 = 1.
:
, откуда х3 = 1; х2 = 1; х1 = 1.