Решение за конечное число арифметических действий. Сюда относятся известное правило Крамера нахождения решения при помощи определите

В.Б. Андреев
ЧИСЛЕННЫЕ МЕТОДЫ
Часть I

2

Глава I
Вычислительные методы линейной алгебры
3

5
С вычислительной точки зрения в линейной алгебре имеются, если понимать их достаточно широко, две основные задачи:
1
?
решение систем линейных уравнений,
2
?
вычисление собственных значений и собственных векторов матрицы.
Основное внимание в лекциях будет уделено решению первой задачи, да и то при весьма ограничительных предположениях. Вторая задача более трудная  ее мы коснемся менее подробно.
В силу теоремы Кронекера - Капелли система линейных алгебраических уравнений
Ax = b разрешима тогда и только тогда, когда ранг матрицы A равен рангу расширенной матрицы [Ab]. Это заведомо так, если матрица A квадратная и невырожденная, т.е.
det A = 0
. В этом случае система не только разрешима при любых b, но и имеет единственное решение. (Разрешима однозначно).
Именно этот случай мы и будем изучать.
Методы решения систем линейных алгебраических уравнений делятся на две груп- пы. К первой группе принадлежат так называемые прямые методы  алгоритмы,
позволяющие получить решение за конечное число арифметических действий. Сюда относятся известное правило Крамера нахождения решения при помощи определите- лей, метод исключения Гаусса, метод прогонки  метод решения систем с трехдиа- гональными матрицами. Существуют и другие методы, из которых отметим метод
Холецкого (метод квадратных корней), применяемый к системам с симметричными положительно определенными матрицами, метод вращений и метод отражений.
Вторую группу составляют приближенные методы, в частности, итерационные.
В итерационных методах решение системы получается как предел при стремлении числа итераций n к бесконечности. При конечных n, как правило, получаются лишь приближенные решения.
Прямые и итерационные методы имеют свою область применения: если размер- ность системы не слишком велика, то часто предпочтительнее использовать прямые методы. Итерационные методы выгодны для систем большого порядка. Особенно в случае матриц специального вида.
В настоящем курсе основное внимание будет уделено прямым методам, а итераци- онных методов коснемся лишь кратко. Более подробно итерационные методы будут изложены во второй части курса "Численные методы", которая читается на четвертом курсе.

6

џ 1
Метод исключения Гаусса и треугольное (LU) разложение матрицы
1.1 Метод исключения Гаусса
Система Ax = b в развернутой форме имеет вид n
j=1
a ij x
j
= b i
,
i = 1, . . . , n.
(1.1)
Как известно, метод Гаусса или метод последовательного исключения неизвестных состоит в том, что неизвестные x j
, j = 1, . . . , n ? 1 последовательно исключаются из соответствующих уравнений системы (1.1) , в результате чего она преобразуется к эквивалентной системе с треугольной матрицей a
(0)
11
x
1
+ a
(0)
12
x
2
+ a
(0)
13
x
3
+ . . . + a
(0)
1n x
n
= b
(0)
1
,
a
(1)
22
x
2
+ a
(1)
23
x
3
+ . . . + a
(1)
2n x
n
= b
(1)
2
,
a
(i?1)
ii x
i
+ . . . + a
(i?1)
in x
n
= b
(i?1)
i
,
a
(n?1)
nn x
n
= b
(n?1)
n
,
(1.2)
коэффициенты a
(k)
ij которой и компоненты ее правой части b
(k)
i вычисляются по фор- мулам a
(k)
ij
= a
(k?1)
ij
? l ik a
(k?1)
kj
,
i, j = k + 1, . . . , n;
k = 1, . . . , n ? 1;
(1.3)
7

8
џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ
b
(k)
i
= b
(k?1)
i
? l ik b
(k?1)
k
,
i, j = k + 1, . . . , n;
k = 1, . . . , n ? 1;
(1.4)
а l
ik
= a
(k?1)
ik
/a
(k?1)
kk
,
i = k + 1, . . . , n;
k = 1, . . . , n ? 1;
(1.5)
причем a
(0)
ij
= a ij
, b
(0)
i
= b i
Вычисления по формулам (1.3)-(1.5)называются прямым ходом метода Гаусса.
После этого неизвестные x k
последовательно, начиная с x n
, находятся из (1.2) по формулам x
i
= b
(i?1)
i
?
n j=i+1
a
(i?1)
ij x
j a
(i?1)
ii
,
i = n, . . . , 1.
(1.6)
Вычисления по этим формулам называют обратным ходом метода Гаусса.
Замечание 1.1. В формулах (1.2) при пребразованиях системы (1.1) первое уравне- ние осталось без изменения. С равным успехом может быть использован и другой вариант исключения, когда первое уравнение (1.1) делится на a
11
, а вместо (1.2)
получается система с единичными коэффициентами при x j
в j-ом уравнении.
Замечание 1.2. Вычисления по формулам (1.5), (1.6) , а, следовательно, и по фор- мулам (1.3), (1.4) возможны лишь тогда, когда все числа a
(i?1)
ii
= 0,
i = 1, . . . , n.
(1.7)
Необходимые и достаточные условия выполнения (1.7) устанавливаются в доказы- ваемой чуть позже теореме 1.2 1.2 LU разложение матрицы.
Покажем, что метод Гаусса эквивалентен разложению матрицы A системы (1.1) в произведение нижней L и верхней U треугольных матриц с последующим решением вспомогательных систем с этими матрицами. В самом деле, из (1.4) находим, что b
(k?1)
i
= b
(k?2)
i
? l i k?1
b
(k?2)
k?1
Подставляя это соотношение в (1.4), получим b
(k)
i
= b
(k?2)
i
? l i k?1
b
(k?2)
k?1
? l ik b
(k?1)
k
Точно так же подставляя сюда выражения для b
(k?2)
i
, а затем для b
(k?3)
i и т.д., будем иметь b
(k)
i
= b
(0)
i
? l i1
b
(0)
1
? l i2
b
(1)
2
? · · · ? l ik b
(k?1)
k

1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ.
9
Полагая здесь k = i ? 1 и выражая b
(0)
i
= b i
через b
(j)
i
, получим b
i
=
i?1
j=1
l ij b
(j?1)
j
+ b
(i?1)
i
(1.8)
Обозначим столбец правой части системы (1.2) через y = [y
1
. . . y n
]
T
, полагая y
i
= b
(i?1)
i
(1.9)
В этих обозначениях (1.8) перепишется так b
i
=
i?1
j=1
l ij y
j
+ y i
,
i = 1, . . . , n.
(1.10)
Обозначим через L нижнюю треугольную матрицу с коэффициентами l ij
, вычисляе- мым по формулам (1.5), и единичной главной диагональю
L =
?
?
?
?
?
?
1 0
0
. . . 0
l
21 1
0
. . . 0
l
31
l
32 1
. . . 0
l n1
l n2
l n3
. . . 1
?
?
?
?
?
?
(1.11)
Тогда (1.10) можно записать в матричном виде b = Ly.
(1.12)
Если верхнюю треугольную матрицу системы (1.2) обозначить через U и переписать
(1.2) в матричном виде, то будем иметь
Ux = y.
(1.13)
Действуя теперь на левую и правую часть (1.13) невырожденной матрицей L и при- нимая во внимание (1.12), получим
LUx = Ly = b
?
A = LU.
(1.14)
Итак, мы показали, что реализация вычислений по формулам (1.3) и (1.5) прямого хода метода Гаусса эквивалентна разложению матрицы A системы (1.1) в произве- дение нижней треугольной матрицы с единичной главной диагональю L и верхней треугольной матрицы U. При этом элементы матрицы L вычисляются по формулам
(1.5), а элементы матрицы U суть u
kj
= a
(k?1)
kj
(1.15)

10
џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ
и вычисляются по формулам (1.3).
После разложения матрицы A в произведение двух треугольных для отыскания решения системы (1.1) нужно решить две системы с треугольными матрицами 
системы (1.12) и (1.13). Решение системы (1.12) заменяет преобразование вектора правой части системы (1.1) по формулам (1.4) прямого хода метода Гаусса. Решение же x системы (1.13) с учетом обозначений (1.9) и (1.15) определяется формулами (1.6)
обратного хода метода Гаусса.
Замечание 1.3. Соотношения (1.3) содержат формулы для u kj
(1.15) и промежуточ- ные значения, которые тоже нужно запоминать и хранить. Мы сейчас преобразуем эти формулы к такому виду, при котором хранение промежуточных значений не требуется.
Пусть матрица L имеет вид (1.11), т.е. ее элементы l
ik
= 0
при k > i,
(1.16)
а
U =
?
?
?
?
?
?
u
11
u
12
u
13
. . . u
1n u
22
u
23
. . . u
2n u
33
. . . u
3n u
nn
?
?
?
?
?
?
,
т.е.
u kj
= 0
при k > j.
(1.17)
Поскольку LU = A, то по правилу умножения матриц находим, что a
ij
=
n k=1
l ik u
kj
(1.18)
Преобразуем эту формулу двумя способами. В силу (1.11), (1.16)
n k=1
l ik u
kj
=
i?1
k=1
l ik u
kj
+ l
=1
ii u
ij
+
n k=i+1 0
l ik u
kj
=
=
i?1
k=1
l ik u
kj
+ u ij
,
а в силу (1.17)
n k=1
l ik u
kj
=
j?1
k=1
l ik u
kj
+ l ij u
jj
+
n k=j+1
l ik
0
u kj
=
=
j?1
k=1
l ik u
kj
+ l ij u
jj

1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ.
11
Отсюда и из (1.18) имеем u
ij
= a ij
?
i?1
k=1
l ik u
kj i = 1, . . . , n;
j = i, . . . , n;
l ij
=
1
u jj a
ij
?
j?1
k=1
l ik u
kj j = 1, . . . , n;
i = j + 1, . . . , n.
(1.19)
Очевидно, что реализация формул (1.19) возможна только тогда, когда все u jj
= a
(j?1)
jj в силу (1.15) отличны от нуля (ср. с (1.7)).
Замечание 1.4. Формулы (1.19) устроены так, что нельзя сначала вычислить все u
ij
, а затем все l ij или наоборот. Можно предложить следующий порядок вычислений по формулам (1.19):
u
1j
= a
1j
,
j = 1, 2, . . . , n;
l i1
= a i1
/u
11
,
i = 2, 3, . . . , n;
u
2j
= a
2j
? l
21
u
1j
,
j = 2, 3, . . . , n;
l i2
= (a i2
? l i1
u
12
)/u
22
,
i = 3, 4, . . . , n;
и т.д., т.е. чередовать вычисление строк матрицы U и столбцов матрицы L.
После построения матриц L и U решение систем (1.12) и (1.13) с треугольными матрицами находятся по формулам y
i
= b i
?
i?1
k=1
l ik y
k
,
i = 1, 2, . . . , n,
(1.20)
(вычисления ведутся сверху вниз).
x k
=
1
u kk y
k
?
n j=k+1
u kj x
j
,
k = n, n ? 1, . . . , 1
(1.21)
(вычисления ведутся снизу вверх).
Одной из важнейших характеристик любого численного метода является его тру- доемкость. Под трудоемкостью метода, предназначенного для решения системы (1.1),
обычно понимают число арифметических действий, необходимых для нахождения искомого решения. Часто в трудоемкость метода включают лишь действия умножения и деления, как наиболее трудоемкие операции с точки зрения работы компьютера. Так будем поступать и мы.

12
џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ
Легко видеть, что для вычислений по формулам (1.19) (для получения треуголь- ного разложения) требуется
Q =
n i=1
n j=i
(i ? 1) +
n j=1
n i=j+1
j =
n i=1
[(i ? 1)(n ? i + 1) + i(n ? i)] =
= 2
n i=1
[(n + 1)i ? i
2
] ? n(n + 1) = n(n + 1)
2
?
n(n + 1)(2n + 1)
3
? n(n + 1) =
=
n(n
2
? 1)
3
=
n
3 3
+ O(n) ?
n
3 3
(1.22)
действий умножения и деления.
Для вычислений по формулам (1.20) и (1.21) имеем соответственно
?
q =
n i=1
(i ? 1) =
n(n ? 1)
2
и q =
n k=1
(n ? k + 1) =
n(n + 1)
2
,
т.е. общее число действий для решения систем (1.12) и (1.13) по формулам (1.20),
(1.21) есть q =
?
q + q = n
2
(1.23)
Замечание 1.5. Из формул (1.22) и (1.23) следует, что при больших n основной объем работы, которую нужно выполнить для решения системы (1.1) описанным методом, падает на преобразование коэффициентов матрицы системы, т.е. на постро- ение треугольного разложения, в то время как преобразование вектора правой части
(решение системы (1.12)) и на отыскание самого решения трудозатраты сравниетельно невелики. В связи с этим при больших n решение нескольких систем с различными правыми частями и одной и той же матрицей оказывается по трудоемкости практи- чески таким же как и решение одной системы.
Выясним теперь условия, при которых вычисления по формулам (1.19)-(1.21) воз- можны, т.е. все u jj отличны от нуля.
Теорема 1.1. Пусть A  невырожденная матрица, L  нижняя треугольная мат- рица с единичной главной диагональю, а U  невырожденная верхняя треугольная матрица. Тогда, если A = LU, то это представление единственно.
Для доказательства теоремы 1.1 нам потребуется
Лемма 1.1. Произведение нижних (верхних)треугольных матриц есть нижняя
(верхняя) треугольная матрица. Обратная к невырожденной нижней (верхней) тре- угольной матрице есть нижняя (верхняя) треугольная матрица.
Упражнение 1.1. Доказать лемму 1.1.

1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ.
13
Доказательство теоремы 1.1. Пусть A = L
1
U
1
= L
2
U
2
. Тогда
L
2
= L
1
U
1
U
?1 2
и
L
?1 1
L
2
= U
1
U
?1 2
Слева стоит произведение нижних треугольных матриц, а справа  верхних. Поэтому произведение есть диагональная матрица D, т.е. L
?1 1
L
2
= D
. Отсюда находим, что
L
2
= L
1
D
. Поскольку главные диагонали L
1
и L
2
единичные, то главная диагональ
L
1
D
совпадает с главной диагональю D, следовательно, D = I. Отсюда L
1
= L
2
и
U
1
= U
2
. Теорема доказана.
Теорема 1.2. Пусть A  квадратная невырожденная матрица, L  нижняя тре- угольная матрица с единичной главной диагональю, а U  невырожденная верхняя треугольная матрица. Разложение A = LU существует тогда и только тогда,
когда все угловые миноры матрицы A отличны от нуля.
Напомним, что угловыми минорами матрицы A называются величины
?
1
= a
11
,
?
2
= det a
11
a
12
a
21
a
22
, . . . ,
?
n
= det[A].
Доказательство. 1
?
. (Необходимость)
Пусть разложение A = LU существует. Тогда по теореме 1.1 оно единственно.
Представим матрицы A, L и U в блочном виде
A =
A
m
A
12
A
21
A
22
,
L =
L
m
0
L
21
L
22
,
U =
U
m
U
12 0
U
22
,
где A
m
, L
m
, U
m и A
22
, L
22
, U
22
 квадратные матрицы размерностей m Ч m и (n ?
m) Ч (n ? m)
соответственно, а m < n  произвольное число. Разложение A = LU в блочном представлении имеет вид
A
m
A
12
A
21
A
22
=
L
m
0
L
21
L
22
U
m
U
12 0
U
22
=
L
m
U
m
L
m
U
12
L
21
U
m
L
21
U
12
+ L
22
U
22
(1.24)
Отсюда следует, что
A
m
= L
m
U
m
(1.25)
Поскольку матрица U треугольная и невырожденная, то все ее диагональные элемен- ты отличны от нуля. Поэтому невырождена и треугольная матрица U
m
. Тем самым
?
m
= det[A
m
] = det[L
m
] det[U
m
] = u
11
. . . u mm
= 0
при m = 1, . . . , n.
2
?
. (Достаточность) Пусть теперь ?
1
?
2
. . . ?
n
= 0
. Для доказательства суще- ствования треугольного разложения воспользуемся методом полной математической индукции по порядку системы n. При n = 1 матрица A = a
11
= ?
1
= 0
, матрица L = 1

14
џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ
и поэтому U = u
11
= a
11
= det U = 0
. Существование искомого разложения при n = 1
доказано.
Пусть A
k
 матрица порядка k и разложение A
k
= L
k
U
k существует с det U
k
= 0
при k = 1, . . . , m ? 1. Докажем, что существует и A
m
= L
m
U
m
, причем det U
m
= 0
Пусть a
·m
= [a
1m
. . . a m?1 m
]
T
 столбец, a m·
= [a m1
. . . a m m?1
]
 строка и разложение
A
m будем искать в виде
A
m
=
A
m?1
a
·m a
m·
a mm
=
L
m?1 0
l m·
1
U
m?1
u
·m
0
u mm
=
=
L
m?1
U
m?1
L
m?1
u
·m l
m·
U
m?1
l m·
u
·m
+ u mm
Отсюда следует, что неизвестный столбец u
·m
, неизвестная строка l m·
и u mm опреде- ляются следующими соотношениями:
L
m?1
u
·m
= a
·m l
m·
U
m?1
= a m·
?
U
T
m?1
l
T
m·
= a
T
m·
,
u mm
= a mm
? l m·
u
·m
(1.26)
Поскольку L
m?1
и U
m?1
невырожденные, из первого и второго соотношений (1.26)
можно найти u
·m и l m·
, соответственно, после чего третье соотношение дает u mm
Существование разложения (1.24) для m доказано. Осталось доказать, что det U
m
= 0
Но с учетом (1.25)
0 = ?
m
= det A
m
= det U
m
,
что и требовалось доказать. Теорема полностью доказана.

џ 2
Ленточные методы
2.1
Метод прогонки
Рассмотрим систему линейных алгебраических уравнений
Ax = b,
(2.1)
матрица которой является трехдиагональной. Запишем эту систему в развернутом виде. Пусть b
1
x
1
+ c
1
x
2
= d
1
,
a
2
x
1
+ b
2
x
2
+ c
2
x
3
= d
2
,
a i
x i?1
+ b i
x i
+ c i
x i+1
= d i
,
a n
x n?1
+ b n
x n
= d n
(2.2)
Алгоритм метода прогонки  метода решения системы (2.2)  состоит в следующем
(см. курс "Введение в численные методы", но, может быть, с другими обозначениями!)
а) Нахождение прогоночных коэффициентов (прямая прогонка) по формулам
?
i
= ?c i
/?
i
,
i = 1, 2, . . . , n ? 1,
?
i
= b i
+ a i
?
i?1
,
i = 2, . . . , n,
?
1
= b
1
,
?
i
= (d i
? a i
?
i?1
)/?
i
,
i = 2, . . . , n,
?
1
= d
1
/b
1
(2.3)
б) Нахождение самого решения (обратная прогонка)
x i
= ?
i x
i+1
+ ?
i
,
i = n ? 1, . . . , 1,
x n
= ?
n
(2.4)
15

16
џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ
Из (2.3),(2.4) следует, что общее число умножений и делений при вычислении коэф- фициентов ?
i и ?
i
Q = 2(n ? 1, )
(2.5)
а при вычислении коэффициентов ?
i и решения x i
q = 3(n ? 1).
(2.6)
Сравнение (2.5),(2.6) с (1.22), (1.23) свидетельствуют о том, что прогонка существенно менее трудоемка по сравнению с общим методом Гаусса. Связано это с тем, что мы явным образом воспользовались тем, что значительная часть элементов матрицы A
равна нулю.
2.2 Ленточные матрицы
Определение 2.1. Матрица A называется ленточной с полушириной ленты p, если ее элементы a ij
= 0
при |i ? j| > p, но существует по крайней мере один элемент a
ij
= 0
при |i ? j| = p
Пример 2.1. Для диагональной матрицы a ij
= 0
при |i ? j| > 0 и, следовательно,
ее полуширина равна нулю. Ее лента состоит из одной диагонали и ширина равна 1.
Условно диагональную матрицу можно изобразить как на рис. 1.
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
Рис. 1.
Пример 2.2. Полная матрица имеет 2n?1 диагоналей. Это и есть ширина ее ленты,
а полуширина будет p = n ? 1.
?
?
?
?
?
?
?
?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
?
?
?
?
?
?
?
?
Рис. 2.

2.3. ЛЕНТОЧНЫЙ ВАРИАНТ ТРЕУГОЛЬНОГО РАЗЛОЖЕНИЯ
17
Пример 2.3. На рис. 3 изображена матрица с полушириной p = 1.
?
?
?
?
?
?
?
?
? ?
? ? ?
? ? ?
? ? ?
? ? ?
? ?
?
?
?
?
?
?
?
?
Рис. 3.
Матрицы такой структуры называются трехдиагональными. Ширина ленты 3.
Пример 2.4. Матрица, изображенная на рис.4, имеет полуширину p = 1 и ширину ленты 2. Матрицы такой структуры называются трапециевидными. Изображенная матрица также называется правой ленточной.
?
?
?
?
?
?
?
?
? ?
? ?
? ?
? ?
? ?
?
?
?
?
?
?
?
?
?
Рис. 4.
Пример 2.5. У матрицы на рис. 5
?
?
?
?
?
?
?
?
? 0 ?
0 ? 0 ?
? 0 ? 0 ?
? 0 ? 0 ?
? 0 ? 0
? 0 ?
?
?
?
?
?
?
?
?
Рис. 5.
полуширина p = 2, ширина 5 и всего три ненулевых диагонали.
Определение 2.2. Рисунок, на котором (звездочками) отмечены позиции, где толь- ко и могут располагаться ненулевые элементы матрицы A, называется ее портретом.
2.3 Ленточный вариант треугольного разложения
Модифицируем алгоритм исключения Гаусса на случай ленточных матриц, т.е. зара- нее отбросим те вычисления, которые заведомо приводят к нулевым элементам. Это

18
џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ
позволит нам сэкономить в трудозатратах на решение системы. Обратимся сразу к варианту, основанному на треугольном разложении матрицы A.
Нам потребуется
Лемма 2.1. Если полуширина матрицы A равна p, то в треугольном разложении
A = LU
полуширина L (U) не больше p.
Доказательство. Пусть a
ij
= 0
при |i ? j| > p.
(2.7)
Докажем, что l
ij
= 0
при i ? j > p.
(2.8)
Для доказательства применим метод полной математической индукции по номерам столбцов матрицы L. При j = 1 из (1.19) находим, что l
i1
= a i1
/a
11
,
i = 2, . . . , n.
Отсюда с учетом (2.7) приходим к (2.8) с j = 1, т.е. l i1
= 0
при i ? 1 > p. Пусть теперь утверждение (2.8) верно для столбцов матрицы L с номерами k = 1, 2, . . . , j ? 1, т.е.
l ik
= 0
при i ? k > p, k = 1, 2, . . . , j ? 1.
(2.9)
Докажем справедливость (2.8) для j-ого столбца. Пусть i ? j > p.
(2.10)
Тогда в силу (1.19) и (2.7)
l ij
=
1
u jj
0
a ij
?
j?1
k=1
l ik u
kj
= ?
1
u jj j?1
k=1
l ik u
kj
(2.11)
Оценим разность индексов i ? k у первых сомножителей под знаком суммы в (2.11).
С учетом (2.10) и (2.11) будем иметь i ? k > p + j ? k p + 1.
Но тогда в силу (2.9) первые сомножители в сумме (2.11) обращаются в нуль и соотношение (2.8) установлено. Лемма доказана.
Преобразуем формулы (1.19)-(1.21) на случай ленточной матрицы A. Сначала вы- ясним, для каких значений индексов i и j нужно проводить вычисления по формулам
(1.19). Так как в силу леммы 2.1
u ij
= 0
при j ? i > p,
(2.12)

2.3. ЛЕНТОЧНЫЙ ВАРИАНТ ТРЕУГОЛЬНОГО РАЗЛОЖЕНИЯ
19
то ненулевые элементы u ij могут быть лишь при j?i p
, т.е. при j p+i
. Аналогично l
ij
= 0
при i ? j > p
(2.13)
и, следовательно, ненулевые элементы l ij могут быть только при i ? j p
, т.е. при i
p + j
. Отсюда u
ij
= 0
i = 1, . . . , n,
j = i, . . . , min[n, p + i],
l ij
= 0
j = 1, . . . , n,
i = j + 1, . . . , min[n, p + j].
Теперь преобразуем суммы в (1.19). В силу (2.12) ненулевые слагаемые в суммах
(1.19) могут быть только при j ? k p
, т.е. при k
j ? p.
а в силу (2.13)  только при i ? k p
, т.е. при k
i ? p.
Объединяя эти неравенства и принимая во внимание, что k  натуральное, будем иметь k
max[1, i ? p, j ? p].
Поскольку в формулах для u ij индексы i, j подчинены ограничению j i
, то в этих формулах k
max[1, j ? p].
В формулах же для l ij наоборот i > j и поэтому в них k
max[1, i ? p].
С учетом сказанного, для ленточной матрицы A с полушириной p формулы (1.19)
принимают вид u
ij
= a ij
?
i?1
k=max[1,j?p]
l ik u
kj i = 1, . . . , n,
j = i, . . . , min[n, p + i],
l ij
=
1
u jj a
ij
?
j?1
k=max[1,i?p]
l ik u
kj j = 1, . . . , n,
i = j + 1, . . . , min[n, p + j].
(2.14)
Преобразуем формулы (1.20) и (1.21). В формулах (1.20) в силу (2.13) i ? k p
, а в формулах (1.21) в силу (2.12) j ? k p
. Поэтому y
i
= b i
?
i?1
k=max[1,i?p]
l ik y
k
,
i = 1, . . . , n,
x k
=
1
u kk y
k
?
min[n,k+p]
j=k+1
u kj x
j
,
k = n, . . . , 1.
(2.15)

20
џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ
2.4 Оценка трудоемкости
Оценим трудоемкость LU - разложения ленточной матрицы A, имеющей полуширину p
, и трудоемкость решения системы (2.1) с такой матрицей. Для этого подсчитаем число умножений и делений, необходимых для реализации формул (2.14) и (2.15). На рис. 6 изображены портреты матриц L и U
n ? p
1
p p + 1
n
?
?
?
?
?
?
?
?
|1
|? 1
|? ? 1
? ? ? 1
? ? ? 1
? ? ? 1
?
?
?
?
?
?
?
?
L
1
p n?p p+1
n
?
?
?
?
?
?
?
?
? ? ?| ?
? ?| ? ?
?| ? ? ?
? ? ?
? ?
?
?
?
?
?
?
?
?
?
U
Рис. 6
Перепишем формулы (2.14), отделив формулы для элементов u ij и l ij
, обведенных треугольниками на рис. 6.
u ij
= a ij
?
i?1
k=1
l ik u
kj
,
i = 1, . . . , p,
j = i, . . . , p,
(2.16)
u ij
= a ij
?
i?1
k=j?p l
ik u
kj
,
j = p + 1, . . . , n,
i = j ? p, . . . , j,
(2.17)
l ij
=
1
u jj a
ij
?
j?1
k=1
l ik u
kj
, j = 1, . . . , p,
i = j + 1, . . . , p,
(2.18)
l ij
=
1
u jj a
ij
?
j?1
k=i?p l
ik u
kj
, i = p + 1, . . . , n,
j = i ? p, . . . , i ? 1.
(2.19)
Подсчитаем число умножений и делений, необходимых для реализации формул (2.16)-
(2.19). Из сравнения (2.16), (2.18) с (1.19) следует, что эти формулы при p = n совпадают. Поэтому, с учетом (1.22)
Q
16
(U
p
) + Q
18
(L
p
) =
p(p
2
+ 1)
3
(2.20)

2.4. ОЦЕНКА ТРУДОЕМКОСТИ
21
Здесь Q  трудоемкость, нижний индекс  номер формулы, а аргумент  объект вычислений.
Далее
Q
17
(u ij
) = i ? 1 ? j + p + 1 = i ? j + p,
Q
17
(U) =
n j=p+1
j i=j?p
(i ? j + p),
Q
19
(l ij
) = j ? 1 ? i + p + 1 + 1 = j ? i + p + 1,
Q
19
(L) =
n i=p+1
i?1
j=i?p
(j ? i + p + 1) =
n j=p+1
j?1
i=j?p
(i ? j + p + 1),
и, следовательно,
Q
17
(U) + Q
19
(L) =
n j=p+1
p +
j?1
i=j?p
(i ? j + p + i ? j + p + 1) =
= 2
n j=p+1
p +
p?1
k=1
k = (n ? p)p(p ? 1).
Отсюда и из (2.20) находим, что общее число умножений и делений при построении
L U
-разложения есть
Q =
p(p + 1)
3
(3n ? 2p ? 1) = n(p
2
+ p) ?
2 3
p
3
+ O(p
2
).
(2.21)
Замечание 2.1. Полуширина полной матрицы p = n ? 1. Подставляя это значе- ние p в найденное выражение, получим выражение для трудоемкости треугольного разложения, совпадающее с (1.22). Полагая же здесь p = 1, получим (2.5).
Обратимся к формулам (2.15).
Q
15
(y i
) =
i ? 1,
i = 1, . . . , p,
i ? 1 ? i + p + 1 = p,
i = p + 1, . . . , n,
Q
15
(y) =
n i=1
Q
15
(y i
) =
p(p ? 1)
2
+ p(n ? p) =
p(2n ? p ? 1)
2
,
Q
15
(x) =
p(2n ? p ? 1)
2
+ n и следовательно q = Q
15
(x) + Q
15
(y) = p(2n ? p ? 1) + n = (2p + 1)n ? p(p + 1)
(2.22)
(ср. с (1.23) при p = n ? 1 и (2.6) при p = 1).
Проанализируем формулы (2.21), (2.22). Рассмотрим три случая.

22
џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ
1
?
. p = O(n)
, например, p = ?n, ? < 1. Тогда
Q ? ?
2
n
3
?
2 3
?
3
n
3
= ?
2 1 ?
2?
3
n
3
Легко проверить, что при 0 < ? < 1 коэффициент при n
3
меньше 1/3.
2
?
. p = o(n)
, но p ? ? при n ? ?. В этом случае
Q ? p
2
n,
q ? 2pn.
В частности, при p =
?
n
Q ? n
2
,
q = 2n
3/2 3
?
. p =
const (p = 1, 2, . . . ).
Q ? p(p + 1)n,
q ? (2p + 1)n.
При p = 1 Q < q, при p
2 Q > q
2.5 LU - разложение для трехдиагональной матрицы.
Линейная система (2.1) с трехдиагональной матрицей A представляет особый интерес в силу того, что часто встречается в приложениях. Получим формулы LU - разложе- ния для этого случая, не апеллируя к (2.14), (2.15). Запишем систему (2.1) в виде (2.2).
Найдем такое LU - разложение матрицы A, в котором U (а не L!) имеет единичную главную диагональ (см. замечание 1.1). Матрица системы (2.2) имеет вид
A =
?
?
?
?
?
?
b
1
c
1 0 . . .
0
a
2
b
2
c
2 0
0 a
3
b
3 0
0 0
0 . . . b n
?
?
?
?
?
?
(2.23)
Пусть
L =
?
?
?
?
?
?
?
1 0
0 . . .
0
?
2
?
2 0 . . .
0 0
?
3
?
3 0
0 0
0 . . . ?
n
?
?
?
?
?
?
,
U =
?
?
?
?
?
?
1 ??
1 0
. . . 0 0
1
??
2
. . . 0 0
0 1
. . . 0 0
0 0
. . . 1
?
?
?
?
?
?
Тогда
LU =
?
?
?
?
?
?
?
1
??
1
?
1
?
2
??
2
?
1
+ ?
2
??
2
?
2
?
3
??
3
?
2
+ ?
3
??
3
?
3
?
4
??
4
?
3
+ ?
4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
?
?
?
?
?
?

2.5. LU - РАЗЛОЖЕНИЕ ДЛЯ ТРЕХДИАГОНАЛЬНОЙ МАТРИЦЫ.
23
Сравнивая эту формулу с матрицей (2.23), находим формулы для элементов матриц
U
и L:
?
i
= ?c i
/?
i
,
i = 1, . . . , n ? 1,
?
i
= b i
+ ?
i
?
i?1
,
i = 2, . . . , n,
?
1
= b
1
,
?
i
= a i
,
i = 2, . . . , n.
(2.24)
Обратимся к решению системы (2.2) с использованием LU - разложения. Пусть
Ux = ?
. Тогда L? = d. Записывая последнюю систему в развернутом виде, будем иметь
?
1
?
1
= d
1
,
?
2
?
1
+ ?
2
?
2
= d
2
,
?
i
?
i?1
+ ?
i
?
i
= d i
,
?
n
?
n?1
+ ?
n
?
n
= d n
Отсюда
?
i
= (d i
? ?
i
?
i?1
)/?
i
,
i = 2, . . . , n,
?
1
=
d
1
?
1
(2.25)
Аналогично находим, что решение системы Ux = ? определяется формулами x
i
= ?
i x
i+1
+ ?
i
,
i ? n ? 1, . . . , 1,
x n
= ?
n
(2.26)
Сравнение формул (2.24) - (2.26) с формулами (2.3),(2.4) показывает, что метод про- гонки является не чем иным как частным случаем ленточного варианта метода ис- ключения Гаусса.
Определение 2.3. Говорят, что матрица A ленточная с лентой нижней полуширины p
1
и верхней полуширины p
2
, если a ij
= 0
при i ? j > p
1
и j ? i > p
2
Упражнение 2.1. Построить модификацию формул (2.14), (2.15) для случая мат- рицы с разной полушириной.

24
џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ

џ 3
Методы Холецкого и блочного исключения. Вычисление обратной матрицы
3.1
Метод Холецкого (квадратных корней)
Вновь обратимся к системе
Ax = b.
(3.1)
На этот раз будем предполагать, что матрица A симметрична и положительно опре- делена, т.е.
A = A
T
и A > 0.
(3.2)
Последнее означает, что квадратичная форма x
T
Ax > 0
для любого ненулевого век- тора x. Напомним, что симметричная матрица имеет только действительные соб- ственные значениия, а положительно определенная  только положительные. В силу критерия Сильвестра необходимым и достаточным условием положительной опреде- ленности матрицы A является положительность всех ее угловых миноров ?
i
> 0,
i = . . . , n.
Построим алгоритм решения системы (3.1), который использует свойства (3.2)
матрицы A. Это будет метод Холецкого. Основой метода Холецкого является
Теорема 3.1. Если A = A
T
> 0
, то существует единственное разложение
A = LL
T
,
(3.3)
где L  нижняя треугольная матрица с положительными диагональными элемен- тами.
Определение 3.1. Разложение (3.3) называется разложением Холецкого, а матрица
L
 множителем Холецкого.
25

26
џ 3. МЕТОДЫ ХОЛЕЦКОГО И БЛОЧНОГО ИСКЛЮЧЕНИЯ
Доказательство теоремы. Сначала докажем единственность. Пусть существуют два разложения
A = L
1
L
T
1
= L
2
L
T
2
Обращая матрицы L
2
и L
T
1
, будем иметь
L
?1 2
L
1
= L
T
2
L
T
1
?1
Принимая во внимание, что (AB)
T
= B
T
A
T
, а для невырожденных матриц (AB)
?1
=
B
?1
A
?1
, найдем, что
L
?1 1
L
2
?1
= L
?1 2
L
1
= L
T
2
L
T
1
?1
= L
?1 1
L
2
T
(3.4)
В силу леммы 1.1 обратная к нижней треугольной матрице есть нижняя треуголь- ная матрица и произведение таких матриц есть снова нижняя треугольная матрица.
Из сказанного следует, что в левой части (3.4) стоит нижняя треугольная матрица, а справа  верхняя. Равенство (3.4) возможно только тогда, когда обе матрицы диаго- нальные. Но диагональная матрица совпадает со своей транспонированной. Поэтому из (3.4) следует, что
L
?1 1
L
2
?1
= L
?1 1
L
2
= D.
Это соотношение утверждает, что диагональная матрица D совпадает со своей об- ратной, что возможно только в том случае, если у этой матрицы диагональными элементами являются числа ±1. Поскольку L
2
= L
1
D
, а диагональные элементы L
1
и
L
2
положительны, то диагональные элементы D тоже должны быть положительны,
т.е. D ? I и, следовательно, L
2
= L
1
. Единственность доказана.
Построим теперь формулы для вычисления элементов L, откуда и будет следовать существование. Так как a ij
= a ji
, а l ij
= 0
при i < j, то будем считать, что i
j.
Тогда a
ij
=
n k=1
l ik l
T
kj
=
n k=1
l ik l
jk
=
j?1
k=1
l ik l
jk
+ l ij l
jj
+
n k=j+1
l ik
0
l jk
=
=
j?1
k=1
l ik l
jk
+ l ij l
jj
При i = j находим, что l
jj
=
a jj
?
j?1
k=1
l
2
jk
,
j = 1, . . . , n.
(3.5)

3.1. МЕТОД ХОЛЕЦКОГО (КВАДРАТНЫХ КОРНЕЙ)
27
Далее,
l ij
=
1
l jj a
ij
?
j?1
k=1
l ik l
jk
,
i = j + 1, . . . , n,
j = 1, . . . , n ? 1
(3.6)
Вычисления можно вести по столбцам j = 1, . . . , n для i = j + 1, . . . , n.
j = 1 :
l
11
=
?
a
11
,
l i1
=
a i1
l
11
,
i = 2, . . . , n j = 2 :
l
22
=
a
22
? l
2 21
,
l i2
=
1
l
22
[a i2
? l i1
l
21
] ,
i = 3, . . . , n и т.д.
Осталось доказать, что все l jj положительны, т.е. положительны подкоренные выражения. Докажем, что l
jj
=
?
j
/?
j?1
,
где ?
0
= 1.
Пусть, как раньше,
A =
A
j
A
12
A
21
A
22
,
L =
L
j
0
L
21
L
22
Тогда
A
j
= L
j
L
T
j
Отсюда
?
j
= det A
j
= (det L
j
)
2
=
j k=1
l kk
2
Аналогично
?
j?1
=
j?1
k=1
l kk
2
и, следовательно,
l
2
jj
= ?
j
/?
j?1
> 0,
j = 1, . . . , n.
Теорема доказана.
Упражнение 3.1. Показать, что для реализации формул (3.5), (3.6) при всех i и j требуется
Q =
n(n + 1)(n + 2)
6
?
n
3 6
(3.7)
операций умножения, деления и извлечения корня.
Замечание 3.1. Из (3.7) следует, что разложение Холецкого в два раза более эко- номично, чем треугольное разложение.

28
џ 3. МЕТОДЫ ХОЛЕЦКОГО И БЛОЧНОГО ИСКЛЮЧЕНИЯ
Обратимся теперь к решению системы (3.1). Поскольку Ax = LL
T
x = b
, то,
полагая L
T
x = y
, получим Ly = b. При этом y
i
=
1
l ii b
i
?
i?1
k=1
l ik y
k
,
i = 1, . . . , n,
x i
=
1
l ii y
i
?
n k=i+1
l ki x
k
,
i = n, . . . , 1.
(3.8)
Замечание 3.2. Очень часто используется модификация разложения Холецкого,
называемая LDL
T
-разложением. Суть ее в том, что вместо разложения (3.3) строится разложение матрицы A вида
A = LDL
T
,
где L  попрежнему нижняя треугольная матрица, но в отличие от (3.3) ее диагональ- ные элементы равны 1, а D  диагональная матрица. Достоинство LDL
T
-разложения состоит в том, что при его вычислении не требуется находить квадратные корни, а потому оно существует не только для положительно определенных матриц. Условием существования такого разложения для симметричной матрицы является отличие от нуля всех ее угловых миноров.
Упражнение 3.2. Построить формулы для вычисления элементов матриц L и D
в LDL
T
-разложении.
Ответ.
d j
= a jj
?
j?1
k=1
l
2
jk d
k
,
j = 1, 2, . . . , n.
l ij
= a ij
?
j?1
k=1
l ik d
k l
jk
/d j
,
i = j + 1, . . . , n, j = 1, . . . , n ? 1.
3.2 Ленточный вариант метода Холецкого
Как и в случае треугольного разложения можно построить разложение Холецкого в ленточном варианте. Пусть A  симметричная положительно определенная ленточ- ная матрица с полушириной p. Справедлива
Лемма 3.1. Если полуширина матрицы A = A
T
> 0
равна p, то и полуширина множителя Холецкого не больше p.
На доказательстве этой леммы мы не останавливаемся, ибо оно почти дословно повторяет доказательство аналогичной леммы 2.1.

3.3. МЕТОД БЛОЧНОГО ИСКЛЮЧЕНИЯ
29
В силу леммы 3.1 ненулевыми элементами l ik матрицы L могут быть только те, у которых индексы подчинены условиям j ? p k
j
(l jk
= 0).
(3.9)
Поэтому формула (3.5) принимает вид l
jj
=
a jj
?
j?1
k=max[1,j?p]
l
2
jk
,
j = 1, . . . , n.
(3.10)
В силу (3.9) в формулах (3.6) i ? p k
i и i ? p j
i
. Поэтому формулы (3.6)
принимают вид l
ij
=
1
l jj
?
?a ij
?
j?1
k=max[1,i?p]
l ik l
jk
?
? ,
i = j + 1, . . . , min[n, p + j],
j = 1, . . . , n.
(3.11)
Ленточный вариант разложения Холецкого построен.
Для отыскания решения системы (3.1) нужно еще преобразовать формулы (3.8).
Они принимают вид y
i
=
1
l ii
?
?b i
?
i?1
k=max[1,i?p]
l ik y
k
?
? , i = 1, . . . , n,
x i
=
1
l ii
?
?y i
?
min[n,i+p]
k=i+1
l ki x
k
?
? , i = n, . . . , 1.
(3.12)
Эти формулы полностью аналогичны формулам (??).
Упражнение 3.3. Подсчитать число действий умножения, деления и извлечения корня, необходимых для реализации формул (3.10)-(3.11).
Ответ.
Q =
(p + 1)(p + 2)(3n ? 2p)
6 3.3 Метод блочного исключения (метод частичного исключения неизвестных)
В методе исключения Гаусса из системы (3.1) последовательно исключаются неиз- вестные  компоненты вектора x
T
= [x
1
, x
2
, . . . , x n
]
. В ряде случаев бывает полезным процедуру исключения неизвестных произвести блочно. Пусть
A =
A
11
A
12
A
21
A
22
,
b =
b
1
b
2
,
x =
x
1
x
2
,
(3.13)

30
џ 3. МЕТОДЫ ХОЛЕЦКОГО И БЛОЧНОГО ИСКЛЮЧЕНИЯ
где A
11
 квадратная невырожденная матрица размеров mЧm, а b
1
и x
1
 m-мерные векторы. С учетом (3.13) система (3.1) принимает вид
A
11
A
12
A
21
A
22
x
1
x
2
=
b
1
b
2
или после блочного перемножения
A
11
x
1
+ A
12
x
2
= b
1
,
A
21
x
1
+ A
22
x
2
= b
2
(3.14)
Из первого уравнения (3.14) находим, что x
1
= A
?1 11
(b
1
? A
12
x
2
).
(3.15)
Подставляя это представление x
1
во второе уравнение (3.14), получим
A
21
A
?1 11
(b
1
? A
12
x
2
) + A
22
x
2
= b
2
или после преобразования
(A
22
? A
21
A
?1 11
A
12
)x
2
= b
2
? A
21
A
?1 11
b
1
(3.16)
В результате система (3.14) преобразовалась к системе
A
11
x
1
+ A
12
x
2
= b
1
,
(A
22
? A
21
A
?1 11
A
12
)x
2
= b
2
? A
21
A
?1 11
b
1
(3.17)
(Неизвестные x
1
исключены из второй группы уравнений).
Из (3.17) вроде бы следует, что для реализации блочного исключения нужно вы- числять A
?1 11
. На самом деле явно это делать вовсе не обязательно. Принимая во внимание (3.15), введем следующие обозначения:
A
?1 11
b
1
=
?
x
1
,
A
?1 11
A
12
= Z
12
(3.18)
Тогда вторая группа уравнений (3.17) примет вид
(A
22
? A
21
Z
12
)x
2
= (b
2
? A
21
?
x
1
).
(3.19)
Соотношения (3.18) можно переписать в виде системы уравнений
A
11
?
x
1
= b
1
,
A
11
Z
12
= A
12
,
(3.20)
а из (3.15) и (3.18) находим, что x
1
=
?
x
1
? Z
12
x
2
(3.21)
Итак, чтобы найти решение системы (3.14) нужно:

3.3. МЕТОД БЛОЧНОГО ИСКЛЮЧЕНИЯ
31 1
?
решить (m + 1) систему (3.20) с матрицей A
11
для отыскания вектора
?
x
1
и столбцов матрицы Z
12
,
2
?
по найденным
?
x
1
и Z
12
сформировать матрицу и правую часть системы (3.19) и решить полученную систему  найти вектор x
2
,
3
?
найти вектор x
1
по формулам (3.21).
Замечание 3.3. В трактовке (3.19),(3.20) метода блочного исключения фактически исключенными оказываются не неизвестные x
1
, а неизвестные x
2
. Из системы (3.14)
как бы исключается часть неизвестных (именно x
2
), затем она решается относительно оставшихся неизвестных (3.20) (но не полностью  нужен еще шаг (3.21)) и лишь потом находится x
2
из (3.19). Отсюда второе название метода  метод частичного исключения неизвестных (исключение x
2
).
Пример 3.1. Пусть матрица A имеет портрет, изображенный на рис. 1. Матрица
A
не является ленточной, хотя ее подматрица A
11
, расположенная в первых (n ? 1)
строках и (n ? 1) столбцах является ленточной с полушириной p = 2. Для решения системы (3.1) с такой матрицей не годится ленточный вариант исключения Гаусса,
а применение общего метода требует O(n
3
)
умножений и делений. Но если можно применить алгоритм блочного исключения,
?
?
?
?
?
?
?
?
?
?
? ? ?
?
? ? ? ?
?
? ? ? ? ?
?
? ? ? ?
?
?
? ? ? ?
? ? ? ? ? · · · ? ? ? ?
?
?
?
?
?
?
?
?
?
?
Рис. 1
то для решения двух систем (3.20) с пятидиагональными матрицами с использованием ленточного варианта исключения потребуется O(n) действий. Столько же действий потребуется для вычислений по формулам (3.19) и (3.21). В результате система будет решена за O(n) действий.
Пример 3.2. Матрица имеет портрет, изображенный на рис. 2.
?
?
?
?
?
?
?
?
?
?
?
?
? ? ? ? ? ? · · · ? ? ?
? ? ?
?
? ? ? ?
?
?
? ? ?
?
?
? ? ?
?
?
? ? ?
? ? ? ? ? ? · · · ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
Рис. 2

32
џ 3. МЕТОДЫ ХОЛЕЦКОГО И БЛОЧНОГО ИСКЛЮЧЕНИЯ
Переставляя первую строку на последнее место и то же делая с первым столбцом,
получим матрицу с портретом, изображенным на рис. 3.
?
?
?
?
?
?
?
?
?
?
?
?
? ?
? ?
? ? ?
? ?
? ? ?
? ?
? ? ?
? ?
? ? ? ?
? ? ? ? ? · · · ? ? ? ?
? ? ? ? ? · · · ? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
Рис. 3
Теперь в качестве A
11
следует выбрать трехдиагональную матрицу размеров (n?2)Ч
(n ? 2)
, стоящую в левом верхнем углу.
3.4 Обращение матрицы
Сразу же заметим, что обращение матрицы в действительности нужно не так часто,
как это может казаться. Предположим,например, что, последующее (за вычислением)
использование A
?1
предусматривает только формирование ее произведений с вектора- ми: u = A
?1
r
, v = A
?1
s и т.д. В таком случае можно было бы вовсе не вычислять A
?1
в явном виде. Вместо этого достаточно провести для A прямой ход метода Гаусса,
а затем находить векторы u, v, . . . как решения линейных систем с матрицей A и правыми частями r, s, . . . . Мы видели в лекции 1, что решение нашей системы требует q ? n
2
операций умножения, но такова же стоимость вычисления произведения A
?1
на вектор. Предварительная же работа резко различается по объему: Q ? n
3
/3
операций умножения для треугольного разложения A и (как мы сейчас покажем) ? n
3
операций для вычисления A
?1
, т.е. во втором случае втрое больше.
Может случиться, что нужны в явном некоторые элементы обратной матрицы,
причем они расположены в одном или нескольких (небольшом числе) столбцов A
?1
Если номера этих столбцов k, l и т.д., то указанные столбцы опять-таки проще всего вычислить как решения систем с матрицей A, в которых правыми частями служат единичные векторы e k
, e l
и т.д. (столбцы единичной матрицы).
И только если необходимы все или большая часть элементов A
?1
, применение процедуры численного обращения A оправдана.
Опишем одну из возможных процедур вычисления A
?1
, основанную на использо- вании LU-разложения. Пусть A = LU. Тогда A
?1
= U
?1
L
?1
или
UA
?1
= L
?1

3.4. ОБРАЩЕНИЕ МАТРИЦЫ
33
Воспользуемся этим соотношением для вычисления A
?1
. Допустим сначала, что L
?1
 известна. Обозначим A
?1
= X
, L
?1
= Y
. Пусть x
.j и y
.j
 j-е столбцы матриц X
и Y , соответственно, т.е. X = [x
.1
x
.2
. . . x
.n
]
, x
.j
= [x
1j x
2j
. . . x nj
]
T
. Тогда получим n систем вида
Ux.
j
= y
.j
,
j = 1, . . . , n
(3.22)
с треугольной матрицей, решения которых могут быть найдены по формулам (1.21)
и x
kj
=
1
u kk y
kj
?
n m=k+1
u km x
mj
,
k = n, n ? 1, . . . 1.
Найдем теперь L
?1
= Y
. Поскольку LY = I, то имеем n систем
Ly.
j
= e j
,
j = 1, . . . , n.
(3.23)
Заметим, что матрица L
?1
 нижняя треугольная и поэтому у столбца y.
j первые
(j ? 1)
элемента известны и равны нулю, т.е. y
T
.j
= [0 . . . 0 y jj y
j+1 j
. . . y nj
] = [0 y
T
.j
]
Отсюда следует, что для отыскания истинных неизвестных вектора y.
j нужно решить систему с треугольной матрицей размеров (n ? j + 1) Ч (n ? j + 1) относительно y
.j
Для этого можно воспользоваться формулами типа (1.21).
Оценим объем работы по вычислению A
?1
. В силу (1.22) факторизация A = LU
требует ? n
3
/3
умножений, решение одной системы (3.22) (см. (1.23) )  ? n
2
/2
, а всех ? n
3
/2
, решение всех систем (3.23)
n j=1
(n ? j + 1)
2 2
=
1 2
n j=1
i
2
=
n(n + 1)(2n + 1)
12
? n
3
/6.
Складывая, находим, что для вычисления матрицы A
?1
при помощи описанного алго- ритма требуется ? n
3
умножений. Это всего лишь в три раза больше, чем для решения системы (3.1).

34
џ 3. МЕТОДЫ ХОЛЕЦКОГО И БЛОЧНОГО ИСКЛЮЧЕНИЯ

џ 4
Устойчивость вычислительных алгоритмов линейной алгебры
4.1
Введение
Исследуем вопрос об устойчивости решения линейной системы по отношению к воз- мущению правой части. Пусть рассматривается система с квадратной невырожденной матрицей
Ax = b
(4.1)
и система с возмущенной правой частью
A?
x = ?b.
(4.2)
Обозначим ?b ? b = ?b, ?x ? x = ?x и оценим ?x через ?b. Вычитая (4.1) из (4.2), будем иметь
A?x = ?b
?
?x = A
?1
?b.
(4.3)
Пусть ·  некоторая норма вектора. В линейной алгебре наиболее часто использу- ются следующие нормы x
?
= max i
|x i
|,
x
1
=
n i=1
|x i
|,
x
2
=
n i=1
x
2
i
Как известно, норма матрицы, подчиненная векторной норме · , определяется со- отношением
A = sup x=0
Ax x
= sup x =1
Ax .
(4.4)
Указанным векторным нормам починены следующие матричные нормы:
A
?
= max i
n j=1
|a ij
|, A
1
= max j
n i=1
|a ij
|, A
2
=
?
max
(AA
T
).
35

36
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ
Очевидно, что для симметричных матриц
A = A
T
,
·
?
= ·
1
,
а A
2
= |?
max
|,
ибо Ax = ?x, A
2
x = ?Ax = ?
2
x
Упражнение 4.1. Доказать, что подчиненные нормы задаются именно этими соот- ношениями.
Из определения (4.4) матричной нормы, в частности, следует, что
A
Ax x
?
Ax
A
x .
(4.5)
Применяя это неравенство ко второму соотношению (4.3), будем иметь
?x
A
?1
?b .
(4.6)
Соотношение (4.6) дает оценку абсолютной погрешности решения через абсолютную погрешность правой части. При этом множителем (коэффициентом усиления) высту- пает норма обратной матрицы. Чем больше эта норма, тем на меньшую точность мы можем рассчитывать.
Получим теперь оценку относительной погрешности. Из (4.1) в силу (4.5)
A
x b .
(4.7)
Деля (4.6) на (4.7), получим
?x x
A
A
?1
?b b
(4.8)
Это и есть оценка относительной погрешности, которая, вообще говоря, несет больше информации, чем оценка (4.6). Здесь коэффициентом усиления выступает число
A
A
?1
=:
cond A,
называемое числом обусловленности матрицы A.
Если число обусловленности матрицы A большое, то про матрицу A говорят, что она плохо обусловлена. В противном случае говорят о хорошо обусловленной матрице.
Поскольку AA
?1
= I
, то A A
?1 1
, т.е. число обусловленности не может быть меньше единицы. Имея систему с хорошо обусловленной матрицей, мы вправе рас- считывать на то, что при небольших возмущениях правой части возмущение решения не будет слишком велико.

4.2. ПРИМЕРЫ СИСТЕМ С ПЛОХО ОБУСЛОВЛЕННОЙ МАТРИЦЕЙ
37 4.2
Примеры систем с плохо обусловленной матри- цей
Пример 4.1.
x
1
= 1,
x
1
+ 0.01x
2
= 1.
(4.9)
Очевидно, что эта система невырождена и ее единственным решением является вектор
[1, 0]
T
. Возмутим правую часть системы (4.9) и найдем решение возмущенной задачи
?
x
1
= 1,
?
x
1
+ 0.01?
x
2
= 1.01,
?b
2
= 0.01.
(4.10)
Очевидно, что решением этой системы является вектор ?x = [1, 1]
T
, который мало похож на невозмущенный вектор x, ибо
?x
2
= 1,
?x
1
= 0,
?x
1
= 1.
Это значение абсолютной погрешности решения полностью согласуется с оценками
(4.6), (4.8), ибо
A
?1
=
1 0
?100 100
,
A
1
= max j
2
i=1
|a ij
| = 2,
A
?1 1
= 101,
x
1
= 1,
b
1
= 2,
?b
1
= 0.01
и, следовательно, в силу (4.6)
?x
1 101 · 0.01 = 1.01,
а в силу (4.8)
?x
1
x
1 2 · 101 ·
0.01 2
= 1.01.
При заданном (4.10) уровне погрешности правой части обусловленность матрицы этой системы (cond A = 202) следует признать плохой.
Пример 4.2.
A =
?
?
?
?
?
?
?
?
1 ?1 ?1 . . . ?1 ?1 0
1 ?1 . . . ?1 ?1 0
0 1 . . . ?1 ?1 0
0 0 . . .
1 ?1 0
0 0 . . .
0 1
?
?
?
?
?
?
?
?
,
b =
?
?
?
?
?
?
?
?
?1
?1
?1
·
?1 1
?
?
?
?
?
?
?
?

38
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ
В развернутом виде система запишется так x
1
?x
2
?x
3
? . . . ?x n
= ?1,
x
2
?x
3
? . . . ?x n
= ?1,
x n?1
?x n
= ?1,
x n
=
1.
(4.11)
Очевидно, что решением системы (4.11) является вектор x = [0, 0, . . . , 0, 1]
T
Легко видеть, что det A = 1.
Возмутим последнюю компоненту вектора b
?b = [?1, ?1, . . . , ?1, 1 + ?]
T
и оценим погрешность решения.
Вычитая из возмущенной системы систему (4.11), для погрешности решения по- лучим
?x
1
??x
2
? . . . ??x n
= 0,
?x
2
? . . . ??x n
= 0,
?x n?1
??x n
= 0,
?x n
= ?.
Отсюда находим, что
?x n
= ?,
?x n?1
= ?,
?x n?2
= ?x n
+ ?x n?1
= 2?,
?x n?3
= ?x n
+ ?x n?1
+ ?x n?2
= 4? = 2 2
?.
Погрешность в каждой из последующих компонент, начиная с ?x n?2
, удваивается, так что
?x n?k
= ?x n
+ ?x n?1
+ · · · + ?x n?(k?1)
= 2
k?1
?,
а
?x
1
= 2
n?2
?.
Таким образом,
?x
?
= 2
n?2
|?|,
x
?
= 1,
?b
?
= |?|,
b
?
= 1,
A
?
= n.
Поскольку в силу (4.8)
?x
?
x
?
cond A
?b
?
b
?
,
а в рассматриваемом случае
?x
?
x
?
= 2
n?2
|?|,

4.3. ПРИМЕР ХОРОШО ОБУСЛОВЛЕННОЙ СИСТЕМЫ
39
то cond A = A
?1
?
A
?
2
n?2
и, следовательно,
A
?1
?
n
?1 2
n?2
При n = 102, A
?
= 102
, cond A
2 100
> 10 30
, A
?1
> 10 27
. Если ? = 10
?15
, то
?x
?
> 10 15
. Матрица рассматриваемой системы очень плохо обусловлена.
Понятие числа обусловленности введено нами только для невырожденных матриц.
Условие det A = 0 означает вырожденность матрицы A, и может сложиться впечатле- ние, что, если det A ? 0, то матрица плохо обусловлена. Однако, прямой связи между величиной определителя матрицы A и ее обусловленностью нет. Так, определитель матрицы из примера (4.2) равен единице, а cond A
2
n?2
С другой стороны, хорошо обусловленная матрица может иметь очень маленький определитель. Например, у матрицы
A =
?
?
?
?
?
10
?1 10
?1 0
0 10
?1
?
?
?
?
?
cond A = 1, хотя det A = 10
?n
4.3 Пример хорошо обусловленной системы
Рассмотрим следующую систему
2x
1
?
9x
2
+
5x
3
=
?4,
1.2x
1
? 5.4x
2
+
6x
3
=
0.6,
x
1
?
x
2
? 7.5x
3
= ?8.5.
(4.12)
Легко проверить, что решением этой системы является вектор x = [0 1 1]
T
Возмутим матрицу и вектор правой части системы (4.12) так, чтобы
?A =
?
?
0 0 0 0 ? 0 0 0 0
?
? , ?b =
?
?
0
?
0
?
? .
(4.13)

40
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ
Очевидно, что решение возмущенной системы совпадает с точным решением, т.е.
?
x = x = [0 1 1]
T
Будем теперь решать систему (4.12) с возмущением (4.13) при ? = 10
?4
на шестираз- рядном десятичном калькуляторе методом последовательных исключений, используя формулы (1.3)  (1.6).
Прямой ход. 1-й шаг. Вычисляем по формуле (1.5) множители l
21
=
a
21
a
11
= 0.6,
l
31
=
a
31
a
11
= 0.5,
являющиеся элементами первого столбца левой треугольной матрицы L, и вычитая из второго и третьего уравнений возмущенной системы первое уравнение, умноженное на l
21
и l
31
, соответственно, т.е. преобразовывая систему при помощи формул (1.3),
(1.4), получим
2 x
1
?
9 x
2
+
5 x
3
=
?4,
0.0001 x
2
+
3 x
3
= 3.0001,
3.5 x
2
? 10 x
3
=
?6.5.
(4.14)
Все вычисления на этом шаге выполнены точно, без округлений, ибо все числа в промежуточных вычислениях, равно как и окончательные числа, имели мантиссы с меньшим, чем шесть, числом разрядов.
2-й шаг. После вычисления по формуле (1.5) множителя l
32
= a
(1)
32
/a
(1)
22
= 3.5/0.0001 = 35000
последнее уравнение системы (4.14) согласно (1.2) должно быть преобразовано к виду a
(2)
33
x
3
= b
(2)
3
,
(4.15)
где согласно (1.3)
a
(2)
33
= a
(1)
33
? l
32
a
(1)
23
= ?10 ? l
32
· 3 = ?10 ? 105000 = ?105010,
а согласно (1.4)
b
(2)
3
= b
(1)
3
? l
32
b
(1)
2
= ?6.5 ? l
32
· 3.0001 =
= ?6.5 ? 105003.5 = ?105010,
т.е.
?105010 x
3
= ?105010.
Однако на используемом калькуляторе будет получено уравнение
?105010 x
3
= ?105011.
(4.16)
В самом деле, коэффициент a
(2)
33
вычисляется точно, так как при его вычислении не возникает чисел, мантиссы которых имеют больше шести разрядов. В то же время

4.3. ПРИМЕР ХОРОШО ОБУСЛОВЛЕННОЙ СИСТЕМЫ
41
при вычислении b
(2)
3
умножение 3.0001 на l
32
дает семизначное число 105003.5, после округления которого до шести разрядов получим 105004. Вычисление b
(2)
3
завершается выполнением операции вычитания b
(2)
3
? ?6.5 ? 105004 = ?105010, 5 ? 105011,
которая также проводится с округлением, что и приводит к (4.16).
Обратный ход. Из (4.16)
?
x
3
= 1.000009522 · · · ? 1.00001.
Сравнение с истинным значением x
3
показывает хорошую точность. Далее, согласно
(1.6), (4.14)
x
2
= (3.0001 ? 3 x
3
)/0.0001 = (3.0001 ? 3.00003)/0.0001 = 0.7.
Здесь все вычисления выполнены точно. Наконец,
?
x
1
= (?4 + 9 x
2
? 5 x
3
)/2 = (?4 + 6.3 ? 5.00005)/2 =
= ?1.350025 ? ?1.35003.
Итак, приближенное решение [?1.35003 0.7 1.00001] мало похоже на точное реше- ние.
В чем причина появления столь значительной погрешности? Говорить о накопле- нии ошибок округления не приходится, так как всего было выполнено 28 арифметиче- ских операций и лишь в четырех случаях потребовалось округление. Предположение о плохой обусловленности системы также не подтверждается, ибо, как показывают вычисления,
A =
?
?
2
?9 5
1.2 ?5.4 6
1
?1
?1.5
?
? , det A = ?21,
A
?1
=
1 21
?
?
?46.5 72.5 27
?15 20 6
?4.2 7
0
?
?
и, следовательно,
A
1
= 18.5,
A
?1 1
=
99.5 21
,
cond A
10 2
Действительная причина состоит в том, что метод исключения Гаусса в том виде,
в каком он был описан, является неустойчивым методом. Чтобы определить, в чем именно его слабость, рассмотрим более внимательно процедуру вычисления b
(2)
3
, где и появились первые округления. При вычислении произведения l
32
· b
(1)
2
= 105003.5
из-за того, что его мантисса содержит более шести знаков, пришлось прибегнуть к округлению. А так как это произведение к тому же оказалось очень большим,
погрешность округления также оказалась больше 0.5.

42
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ
4.4 Метод Гаусса с выбором ведущего элемента
Из-за отмеченной неустойчивости метод Гаусса в вычислительной практике обычно применяется в сочетании с некоторой схемой выбора ведущего элемента. Например,
схема выбора ведущего элемента по столбцу состоит в следующем. Перед началом первого шага среди коэффициентов a
11
, a
21
, . . . , a n1
, образующих первый столбец мат- рицы A, выбирается коэффициент с наибольшим модулем; пусть это будет a k,1
. Если k > 1
, то в системе (4.1) переставляются 1-е и k-е уравнения, при k = 1 перестановка не нужна. После этой предварительной работы обычным образом проводится 1-й шаг прямого хода. До начала 2-го шага среди коэффициентов a
(1)
22
, a
(1)
31
, . . . , a
(1)
n2
(т.е.
во втором столбце текущей матрицы) выбирается коэффициент a
(1)
l2
с наибольшим модулем. В случае l > 2 переставляются 2-е и l-е уравнения, затем выполняется 2-й шаг. И т.д.
Что достигается выбором ведущего элемента по столбцу? Мы можем теперь га- рантировать, что множители l ij всех шагов по абсолютной величине ограничены еди- ницей. Формулы (1.3) показывают, что, во-первых, добавки l ik a
(k?1)
kj к текущим зна- чениям коэффициентов имеют тот же порядок величины, что и сами коэффициенты,
во-вторых, за один шаг уровень коэффициентов матрицы может вырасти не более,
чем в два раза. Действительно, согласно (1.3)
a
(k)
ij
= a
(k?1)
ij
? l ik a
(k?1)
kj a
(k?1)
ij
+ a
(k?1)
kj
2 max ij a
(k?1)
ij
(4.17)
Упражнение 4.2. Решить систему (4.12) методом Гаусса с выбором ведущего эле- мента по столбцу на 6-разрядном десятичном калькуляторе.
Иногда используется и другая схема выбора ведущего элемента, а именно, схема выбора по строке. Здесь до начала 1-го шага определяется наибольший по модулю среди коэффициентов a
11
, a
12
, . . . , a
1n
. Пусть им будет коэффициент a
1k
. Если k > 1, то производится перенумерация неизвестных: 1-е и k-е неизвестные меняются номерами.
Это соответствует перестановке столбцов матрицы системы. При k = 1 перестановка не нужна. Теперь обычным образом проводится 1-й шаг прямого хода. И т.д.
Переходя к выбору ведущего элемента по столбцу, мы получили для системы (4.12)
приближенное решение хорошего качества. Но это не значит, что описанные схемы с выбором ведущего элемента придают методу Гаусса гарантированную устойчивость.
Хотя обычно схемы с выбором ведущего элемента по столбцу или по строке действи- тельно обеспечивают устойчивое вычисление.
В каких же случаях утрачивается устойчивость? Чтобы понять это, заметим, что во многих численных методах ошибки промежуточных вычислений в совокупности равносильны тому, как если бы тем же методом точно решали исходную задачу,
предварительно изменив ее входные данные. Это относится и к методу Гаусса. Можно показать, что решение линейной системы (4.1), вычисленное методом Гаусса (с той или иной схемой выбора ведущего элемента или вообще без выбора) при наличии ошибок

4.4. МЕТОД ГАУССА С ВЫБОРОМ ВЕДУЩЕГО ЭЛЕМЕНТА
43
округления точно удовлетворяет измененному уравнению
(A + ?A)?
x = b.
(4.18)
В пояснение сказанного рассмотрим умножение двух треугольных матриц с учетом ошибок округления
AB =
a
11
a
12 0
a
22
b
11
b
12 0
b
22
=
=
a
11
b
11
(1 + ?
1
) (a
11
b
12
(1 + ?
2
) + a
12
b
22
(1 + ?
3
))(1 + ?
4
)
0
a
22
b
22
(1 + ?
5
)
Если положить
A =
a
11
a
12
(1 + ?
3
)(1 + ?
4
)
0
a
22
(1 + ?
5
)
,
B =
b
11
(1 + ?
1
) b
12
(1 + ?
2
)(1 + ?
4
)
0
b
22
,
то легко проверить, что
AB = AB.
Для нормы матрицы возмущения в (4.18) справедлива оценка
?A
?
f (n)g(A) A
?
p
?t
(4.19)
Здесь f(n)  некоторая медленно растущая функция от порядка n системы (типа степенной с небольшим показателем); p  основание машинной арифметики; t  число разрядов, отведенных для представления мантиссы. Чтобы определить оставшийся сомножитель g(A), обозначим a = a
0
= max ij
|a ij
|,
a k
= max ij>k
|a
(k)
ij
|,
k = 1, 2, . . . , n ? 1.
Тогда g(A) = max
0 6k6n?1
a k
/a.
Таким образом, величина g(A) измеряет насколько могли вырасти элементы проме- жуточных матриц метода по сравнению с уровнем элементов в исходной матрице A.
Итак, при прочих равных условиях ошибки ?A, вносимые методом в исходную информацию, тем больше, чем больший рост элементов матриц допускается в прямом ходе. Если в версии метода Гаусса из первой лекции рост элементов может быть сколь угодно велик, то в схемах с выбором ведущего элемента он ограничен. Действитель- но, за шаг максимальный модуль элемента матрицы может вырасти самое большее в два раза (см. (4.17)). Так как маловероятно, чтобы возмущение происходило на каждом шаге, то обычно коэффициент роста g(A) невелик. В этом случае благодаря множителю p
?t в правой части (4.19) матрицу-возмущение ?A можно считать малой по сравнению с A.
На вопрос о том, как сказывается возмущение (4.18) на решении системы (4.1),
отвечает

44
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ
Теорема 4.1. Пусть матрица A системы (4.1) имеет обратную и для ее возму- щения ?A справедлива оценка
?A < A
?1 ?1
(4.20)
Тогда для относительной погрешности решения справедлива оценка
?x x
cond A
1 ?
cond A
?A
A
?A
A
+
?b b
,
(4.21)
где
(A + ?A)(x + ?x) = b + ?b
(4.22)
 возмущенная система.
Доказательство. Из (4.22)
Ax + ?Ax + A?x + ?A?x = b + ?b.
Вычитая отсюда (4.1), получим
A?x = ?b ? ?Ax ? ?A?x,
или
?x = A
?1
[?b ? ?Ax ? ?A?x].
Отсюда
?x
A
?1
( ?b + ?A
x + ?A
?x ).
Разрешим это неравенство относительно ?x
(1 ? A
?1
?A ) ?x
A
?1
( ?b + ?A
x ).
С учетом (4.20)
?x
A
?1
( ?b + ?A
x )
1 ? A
?1
?A
(4.23)
Но
1 ? A
?1
?A = 1 ?
cond A
?A
A
Учитывая это и деля (4.23) на x , будем иметь
?x x
A
?1
A
?b
A
x
+
?A
A
1 ?
cond A
?A
A
(4.24)
Поскольку A x b
, то из (4.24) вытекает (4.21). Теорема доказана.

4.4. МЕТОД ГАУССА С ВЫБОРОМ ВЕДУЩЕГО ЭЛЕМЕНТА
45
Приведем примеры матриц, преобразование которых при помощи метода Гаусса с выбором ведущего элемента приводит к максимально возможному увеличению коэф- фициентов промежуточных матриц прямого хода.
Пример 4.3. Выбор по столбцу
A =
?
?
?
?
1 0
0 1
?1 1
0 1
?1 ?1 1
1
?1 ?1 ?1 1
?
?
?
? ,
A
1
=
?
?
?
?
1 0
0 1
0 1
0 2
0 ?1 1
2 0 ?1 ?1 2
?
?
?
? ,
A
2
=
?
?
?
?
1 0 0
1 0 1 0
2 0 0 1
4 0 0 ?1 4
?
?
?
? ,
A
3
= U =
?
?
?
?
1 0 0 1 0 1 0 2 0 0 1 4 0 0 0 8
?
?
?
? ,
L =
?
?
?
?
1 0
0 0
?1 1
0 0
?1 ?1 1
0
?1 ?1 ?1 1
?
?
?
? .
Упражнение 4.3. Показать, что выбор ведущего элемента по строке для матрицы
A =
?
?
?
?
1 ?1 ?1 ?1 0
1
?1 ?1 0
0 1
?1 1
1 1
1
?
?
?
?
приводит к максимальному росту элементов промежуточных матриц.
Замечание 4.1. Возможна еще одна схема выбора ведущего элемента: выбор по всей матрице. В этом случае гарантируется полная устойчивость метода Гаусса, однако сама процедура выбора такого ведущего элемента очень трудоемка  для ее реа- лизации требуется O(n
3
)
действий, что сравнимо с трудоемкостью самого метода и,
следовательно, существенно удорожает решение. Отметим, что при выборе ведущего элемента по столбцу или по строке требуется лишь O(n
2
)
дополнительных операций.

46
џ 4. УСТОЙЧИВОСТЬ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ

џ 5
Методы вращений и отражений
5.1
Метод вращений
Вновь обратимся к решению системы
Ax = b
(5.1)
с невырожденной матрицей. Рассмотрим метод вращений решения системы (5.1),
который позволяет получить представление матрицы A в виде произведения орто- гональной матрицы Q и верхней треугольной матрицы R.
Как и в методе Гаусса, в методе вращений на первом шаге неизвестное x
1
исклю- чается из всех уравнений кроме первого. Для того, чтобы исключить x
1
из второго уравнения, умножим первое уравнение на некоторое число c
12
, а второе  на s
12
и заменим первое уравнение суммой вновь полученных уравнений. Затем умножим первое уравнение на s
12
, второе  на c
12
и вычтем первое из второго:
(c
12
a
11
+ s
12
a
21
)x
1
+ (c
12
a
12
+ s
12
a
22
)x
2
+ . . . = c
12
b
1
+ s
12
b
2
,
(?s
12
a
11
+ c
12
a
21
)x
1
+ (?s
12
a
12
+ c
12
a
22
)x
2
+ . . . = ?s
12
b
1
+ c
12
b
2
Числа c
12
и s
12
выберем из условий c
2 12
+ s
2 12
= 1,
?s
12
a
11
+ c
12
a
21
= 0.
(5.2)
Решение нелинейной системы (5.2) при a
2 11
+ a
2 21
= 0
дают, например, формулы c
12
=
a
11
a
2 11
+ a
2 21
,
s
12
=
a
21
a
2 11
+ a
2 21
(5.3)
Если же a
2 11
+a
2 21
= 0
, то решение (5.2) дают числа c
12
= 1
, s
12
= 0
. Второе из уравнений
(5.2) как раз и означает, что в новом втором уравнении неизвестного x
1
не будет.
47

48
џ 5. МЕТОДЫ ВРАЩЕНИЙ И ОТРАЖЕНИЙ
В результате система (5.1) преобразуется к виду a
(1)
11
x
1
+ a
(1)
12
x
2
+ a
(1)
13
x
2
+ . . . + a
(1)
1n x
n
= b
(1)
1
,
a
(1)
22
x
2
+ a
(1)
23
x
3
+ . . . + a
(1)
2n x
n
= b
(1),
2
a
31
x
1
+ a
32
x
2
+ a
33
x
3
+ . . . + a
3n x
n
= b
3
,
a n1
x
1
+ a n2
x
2
+ a n3
x
3
+ . . . + a nn x
n
= b n
,
(5.4)
в котором a
(1)
1j
= c
12
a
1j
+ s
12
a
2j
,
a
(1)
2j
= ?s
12
a
1j
+ c
12
a
2j
,
j = 1, 2, . . . , n,
b
(1)
1
= c
12
b
1
+ s
12
b
2
,
b
(1)
2
= ?s
12
b
1
+ c
12
b
2
(5.5)
Как легко видеть, преобразование исходной системы (5.1) к виду (5.4) эквивалент- но умножению матрицы системы A и вектора правой части b слева на матрицу T
12
,
имеющую вид
T
12
=
?
?
?
?
?
?
?
c
12
s
12
?s
12
c
12 0
0 1
1
?
?
?
?
?
?
?
Для исключения неизвестного x
1
из третьего уравнения системы (5.4) новое первое уравнение и третье уравнение, которое пока не было затронуто преобразованиями,
заменяются их линейными комбинациями с коэффициентами c
13
, s
13
и ?s
13
, c
13
, соот- ветственно. При этом коэффициенты c
13
и s
13
находятся из системы c
2 13
+ s
2 13
= 1,
?s
13
a
(1)
11
+ c
13
a
31
= 0
(5.6)
и суть c
13
=
a
(1)
11
a
(1)
11 2
+ a
2 31
,
s
13
=
a
31
a
(1)
11 2
+ a
2 31
,
а коэффициенты преобразованных первого и третьего уравнений находятся по фор- мулам a
(2)
1j
= c
13
a
(1)
1j
+ s
13
a
3j
,
a
(1)
3j
= ?s
13
a
(1)
1j
+ c
13
a
3j
,
j = 2, 3, . . . , n,
b
(2)
1
= c
13
b
(1)
1
+ s
13
b
3
,
b
(1)
3
= ?s
13
b
(1)
1
+ c
13
b
3

5.1. МЕТОД ВРАЩЕНИЙ
49
Это преобразование эквивалентно умножению слева матрицы системы (5.4) и ее век- тора правой части на матрицу
T
13
=
?
?
?
?
?
?
?
?
?
c
13 0 s
13 0
1 0
?s
13 0 c
13 0
0 1
1
?
?
?
?
?
?
?
?
?
и приводит к тому, что коэффициент при x
1
в преобразованном третьем уравнении обращается в нуль.
Продолжая подобным образом, мы исключим x
1
из всех остальных уравнений.
Первое уравнение изменяется на каждом таком "малом"шаге, которых будет n ? 1.
Поэтому, по завершении первого шага система (5.1) примет вид a
(n?1)
11
x
1
+ a
(n?1)
12
x
2
+ a
(n?1)
13
x
3
+ . . . + a
(n?1)
1n x
n
= b
(n?1)
1
,
a
(1)
22
x
2
+ a
(1)
23
x
3
+ . . . + a
(1)
2n x
n
= b
(1)
2
,
a
(1)
n2
x
2
+ a
(1)
n3
x
3
+ . . . + a
(1)
nn x
n
= b
(1)
n
Упражнение 5.1. Вывести соотношения для a
(l?1)
1j
, a
(1)
lj
, b
(l?1)
1
, b
(1)
l и преобразующих коэффициентов c
1l и s
1l
, l = 2, . . . , n; j = 1, . . . , n.
Ответ.
a
(l?1)
1j
= c
1l a
(l?2)
1j
+ s
1l a
lj
,
a
(1)
lj
= ?s
1l a
(l?2)
1j
+ c
1l a
lj
,
j = 1, 2, . . . , n,
b
(l?1)
1
= c
1l b
(l?2)
1
+ s
1l b
l
,
b
(1)
l
= ?s
1l b
(l?2)
1
+ c
1l b
l
,
l = 2, 3, . . . , n,
c
1l
=
a
(l?2)
11
a
(l?2)
11 2
+ a
2
l1
,
s
1l
=
a l1
a
(l?2)
11 2
+ a
2
l1
,
l = 2, 3, . . . , n,
a
(0)
11
= a
11
Замечание 5.1. Поскольку в силу невырожденности матрицы A по крайней мере один из коэффициентов a i1
= 0
, то a
(n?1)
11
= 0
В матричной записи эта система имеет вид
A
(1)
x = b
(1)
,

50
џ 5. МЕТОДЫ ВРАЩЕНИЙ И ОТРАЖЕНИЙ
где
A
(1)
= T
1
A,
b
(1)
= T
1
b,
T
1
= T
1n
T
1 n?1
. . . T
13
T
12
На втором шаге метода вращений из третьего, четвертого и т.д. уравнений полу- ченной системы исключается неизвестное x
2
. Шаг состоит из (n ? 2) "малых"шагов,
и в каждом из них второе уравнение комбинируется с одним из нижележащих. После выполнения второго шага система преобразуется к виду a
(n?1)
11
x
1
+ a
(n?1)
12
x
2
+ a
(n?1)
13
x
3
+ . . . + a
(n?1)
1n x
n
= b
(n?1)
1
,
a
(n?1)
22
x
2
+ a
(n?2)
23
x
3
+ . . . + a
(n?2)
2n x
n
= b
(n?1)
2
,
a
(2)
33
x
3
+ . . . + a
(2)
3n x
n
= b
(2)
3
,
a
(2)
n3
x
3
+ . . . + a
(2)
nn x
n
= b
(2)
n
Упражнение 5.2. Вывести соотношения для a
(l?1)
2j
, a
(2)
lj
, b
(l?1)
2
, b
(2)
l
, c
2l и s
2l
Ответ.
a
(l?1)
2j
= c
2l a
(l?2)
2j
+ s
2l a
(1)
lj
,
a
(2)
lj
= ?s
2l a
(l?2)
2j
+ c
2l a
(1)
2j
,
j = 2, 3, . . . , n,
b
(l?1)
2
= c
2l b
(l?2)
2
+ s
2l b
(1)
l
,
b
(2)
l
= ?s
2l b
(l?2)
2
+ c
2l b
(1)
l
,
l = 3, 4, . . . , n,
c
2l
=
a
(l?2)
22
a
(l?2)
22 2
+ a
(1)
l2 2
,
s
2l
=
a
(1)
l2
a
(l?2)
22 2
+ a
(1)
l2 2
В матричной форме эта система имеет вид
A
(2)
x = b
(2)
,
где
A
(2)
= T
2
A
(1)
,
b
(2)
= T
2
b
(1)
,
T
2
= T
2n
T
2 (n?1)
. . . T
24
T
23
После (n ? 1) шагов получим систему a
(n?1)
11
x
1
+ a
(n?1)
12
x
2
+ a
(n?1)
13
x
3
+ . . . + a
(n?1)
1n x
n
= b
(n?1)
1
,
a
(n?1)
22
x
2
+ a
(n?2)
23
x
3
+ . . . + a
(n?2)
2n x
n
= b
(n?1)
2
,
a
(n?1)
nn x
n
= b
(n?1)
n
,
(5.7)

5.1. МЕТОД ВРАЩЕНИЙ
51
где a
(l?1)
kj
= c kl a
(l?2)
kj
+ s kl a
(k?1)
lj
,
a
(k)
lj
= ?s kl a
(l?2)
kj
+ c kl a
(k?1)
lj
,
j = k, k + 1, . . . , n,
b
(l?1)
k
= c kl b
(l?2)
k
+ s kl b
(k?1)
l
,
b
(k)
l
= ?s kl b
(l?2)
k
+ c kl b
(k?1)
l
,
k = 1, . . . , n,
l = k + 1, . . . , n,
(5.8)
а c
kl
=
a
(l?2)
kk a
(l?2)
kk
2
+ a
(k?1)
lk
2
,
s kl
=
a
(k?1)
lk a
(l?2)
kk
2
+ a
(k?1)
lk
2
(5.9)
В матричной записи полученная система имеет вид
A
(n?1)
x = b
(n?1)
,
где
A
(n?1)
= T
n?1
A
(n?2)
,
b
(n?1)
= T
n?1
b
(n?2)
,
T
n?1
= T
n?1n
Обозначим через R полученную верхнюю треугольную матрицу A
(n?1)
. Она свя- зана с исходной матрицей A равенством
R = T A,
(5.10)
где T = T
n?1
T
n?2
. . . T
1
Замечание 5.2. Действие матрицы T
kl на вектор x эквивалентно его повороту против хода часовой стрелки в координатной плоскости Ox k
x l
на угол ?
kl такой, что cos ?
kl
= c kl
,
sin ?
kl
= s kl
Существование такого угла гарантируется соотношениями (5.9). Эта геометрическая интерпретация и дала название методу.
С учетом (5.2), (5.6) и т.д. легко видеть, что
T
kl
T
T
kl
= I,
т.е. матрицы T
kl ортогональные. Произведение ортогональных матриц есть матрица ортогональная, и поэтому T есть ортогональная матрица, равно как и T
T
= T
?1
= Q
Отсюда и из (5.10) находим, что
A = QR.
(5.11)
Упражнение 5.3. Показать, что для построения разложения (5.11) с использова- нием формул (5.8), (5.9), требуется ? 4n
3
/3
действий умножения.

52
џ 5. МЕТОДЫ ВРАЩЕНИЙ И ОТРАЖЕНИЙ
Замечание 5.3. Принимая во внимание результаты вычислений из упражнения 5.3,
заключаем, что метод вращений примерно в четыре раза более трудоемок, чем метод
Гаусса.
Выясним, какими же достоинствами обладает этот метод по сравнению с мето- дом Гаусса, благодаря которым он заслуживает право на существование несмотря на большую трудоемкость.
С учетом (5.2) из (5.5) находим, что
[a
(1)
1j
]
2
+ [a
(1)
2j
]
2
= c
2 12
a
2 1j
+ s
2 12
a
2 2j
+ 2c
12
s
12
a
1j a
2j
+
+ s
2 12
a
2 1j
+ c
2 12
a
2 2j
? 2s
12
c
12
a
1j a
2j
= a
2 1j
+ a
2 2j
,
j = 1, 2, . . . , n,
т.е. сумма квадратов первых элементов j-го столбца не изменилась. Если учесть,
что на первом малом шаге коэффициенты остальных уравнений не изменились, то полученное равенство означает: длина любого столбца матрицы не изменилась. Точно так же из формул (5.6), (5.7) выводим, что
[a
(2)
1j
]
2
+ [a
(1)
3j
]
2
= [a
(1)
1j
]
2
+ a
2 3j
,
j = 1, . . . , n,
т.е. длины столбцов неизменны и на 2-ом малом шаге. Это верно и по отношению к каждому из последующих шагов, и по отношению к прямому ходу в целом. Таким образом, в отличие от метода Гаусса метод вращений застрахован от роста элементов промежуточных матриц: на протяжении всего процесса исключения абсолютная вели- чина коэффициента a
(k)
ij не может превосходить длину j-го столбца исходной матрицы
A
. Как следствие, матрица ?A  матрица эквивалентных возмущений, учитывающая ошибки промежуточных возмущений, для метода вращений будет мала, т.е. метод всегда устойчив.
5.2 Метод отражений
Рассмотрим еще один метод, который дает разложение матрицы A в виде произведе- ния ортогональной и верхней треугольной матриц. Это будет метод отражений.
Пусть w  некоторый вектор (столбец) единичной длины w
2 2
= (w, w) = w
T
w = 1.
(5.12)
Введем в рассмотрение матрицу
U = I ? 2ww
T
,
(5.13)
которую назовем матрицей отражения и изучим ее свойства.

5.2. МЕТОД ОТРАЖЕНИЙ
53 1
?
Матрица U симметрична, т.е.
U = U
T
(5.14)
В самом деле, так как
(ww
T
)
T
= (w
T
)
T
w
T
= ww
T
,
то ww
T
есть симметричная матрица, а в силу (5.13) вместе с ней симметричной является и матрица U.
2
?
. Матрица U есть ортогональная матрица, т.е.
U
?1
= U
T
(5.15)
Принимая во внимание (5.14), (5.13) и (5.12), имеем
UU
T
= UU = (I ? 2ww
T
)(I ? 2ww
T
) =
= I ? 4ww
T
+ 4w w
T
w
1
w
T
= I,
что и означает справедливость (5.15).
3
?
. Число ? = ?1 является однократным собственным значением матрицы U, кото- рому отвечает собственный вектор w из (5.13). Число ? = 1 является (n ? 1)-кратным собственным значением матрицы U, которому отвечает (n ? 1)-мерное собственное подпространство, состоящее из всех векторов v, ортогональных w.
В силу (5.14), (5.15) имеем U
2
= UU
T
= I
. Так как все собственные значения матрицы I равны 1, то собственные значения матрицы U удовлетворяют соотношению
?
2
U
= 1
, и, следовательно, равны либо +1, либо ?1.
Далее, принимая во внимание (5.13), (5.12), находим, что
Uw = (I ? 2ww
T
)w = w ? 2w w
T
w
1
= ?w.
(5.16)
Наконец, пусть
(v, w) = w
T
v = 0.
Тогда
Uv = (I ? 2ww
T
)v = v ? 2ww
T
v = v,
(5.17)
что и требовалось доказать.
4
?
Вектор Uy есть зеркальное отражение вектора y относительно плоскости, ор- тогональной вектору w.
Пусть y  произвольный вектор. Представим его в виде y = z + v,
(5.18)

54
џ 5. МЕТОДЫ ВРАЩЕНИЙ И ОТРАЖЕНИЙ
где z  проекция y на w, т.е. z = (w, y)w, а вектор v ортогонален w: (w, v) = 0).
Принимая во внимание (5.16), (5.17), находим, что
Uy = U(z + v) = ?z + v
(5.19)
(см. рис.1)
7
-
6 6
?
w w
z y
v
?z
Uy
Рис. 1 5
?
Векторы y ? Uy и y + Uy ортогональны друг другу. При этом, если y = z + v,
как в (5.18), то y ? Uy = 2z,
z w,
(5.20)
y + Uy = 2v,
v ? w.
(5.21)
Утверждения следуют из (5.18), (5.19).
6
?
Пусть x и y  произвольные векторы единичной длины. Тогда, если w =
y +
sign (x, y) x y +
sign (x, y) x
2
,
(5.22)
то матрица отражения U, построенная по вектору w, переводит вектор y в вектор,
коллинеарный вектору x.
В самом деле, пусть матрица отражения U обладает искомым свойством, т.е. Uy =
??x
. Найдем вектор w, образующий эту матрицу. Согласно свойству 5
?
(см. (5.20))
искомый вектор w коллинеарен вектору y ? Uy, а поскольку w
2
= 1
, то w =
y + ?x y + ?x
2
(5.23)
В силу 2
?
матрица U ортогональна и поэтому
Uy
2
= y
2
= 1 = |?| x
2
= |?|,

5.2. МЕТОД ОТРАЖЕНИЙ
55
т.е. ? = ±1. Если y = ±x, то знаменатель (5.23) будет равен нулю при ? = ?sign (x, y) =
1
. Этого не произойдет, если ? = sign (x, y). Если вектор y близок к x, то мы будем избавлены от неприятностей, связанных с делением на малое число, выбрав и здесь
? =
sign (x, y). Из сказанного следует, что в общем случае целесообразно выбирать ?
согласно (5.22), а если (y, x) = 0, то, например, ? = 1.
Воспользуемся матрицей отражения для приведения квадратной матрицы к тре- угольному виду. На первом шаге приведения рассмотрим в качестве вектора y из свойства 6
?
нормированный первый столбец матрицы A
y
1
= [a
11
a
21
. . . a n1
]
T
/
n i=1
a
2
i1
,
(5.24)
а в качестве x  вектор e
1
= [1 0 . . . 0]
T
. Если a
21
= a
31
= · · · = a n1
= 0
, то переходим к следующему шагу, положив A
(1)
= A
, U
1
= I
и введя обозначения a
(1)
ij
= a ij
. В
противном случае умножим матрицу A слева на матрицу отражения
U
1
= I ? 2w
1
w
T
1
= I
n
? 2w
1
w
T
1
,
(5.25)
где вектор w
1
вычисляется согласно формуле (5.22)
w
1
=
y
1
+
sign(e
1
, y
1
)e
1
y
1
+
sign(e
1
, y
1
)e
1 2
(5.26)
В результате получим матрицу
A
(1)
= U
1
A,
в первом столбце которой стоят нули во всех позициях, кроме первой. Этим заканчи- вается первый этап.
Пусть мы уже осуществили l?1 > 0 шагов и пришли к матрице A
(l?1)
с элементами a
(l?1)
ij такими, что a
(l?1)
ij
= 0
при i > j, j = 1, . . . , l ? 1. В пространстве R
n?l+1
векторов размерности n ? l + 1 рассмотрим вектор y
l
= [a
(l?1)
ll a
(l?1)
l+1,l
. . . a
(l?1)
nl
]
T
/ (a
(l?1)
ll
)
2
+ · · · + (a
(l?1)
nl
)
2
Если a
(l?1)
l+1,l
= a
(l?1)
l+2,l
= · · · = a
(l?1)
nl
= 0
, то переходим к следующему шагу, положив
A
(l)
= A
(l?1)
,
U
l
= I.
В противном случае строим матрицу отражения
V
l
= I
n?l+1
? 2w l
w
T
l
(размеры матрицы V
l и вектора w l
равны (n?l +1)), переводящую вектор y l
в вектор,
коллинеарный e l
= [1 0 . . . 0]
T
? R
n?l+1
, и переходим к матрице
A
(l)
= U
l
A
(l?1)
,
(5.27)

56
џ 5. МЕТОДЫ ВРАЩЕНИЙ И ОТРАЖЕНИЙ
где
U
l
=
I
l?1 0
0
V
l
После (n ? 1) шагов мы приходим к матрице
A
(n?1)
= U
n?1
U
n?2
. . . U
1
A,
имеющей правую треугольную форму. Обозначим
U
n?1
. . . U
1
= U.
Тогда
A
(n?1)
= UA,
A = U
T
A
(n?1)
Если нужно решить систему (5.1), то после описанных преобразований приходим к эквивалентной системе
A
(n?1)
x = Ub
(5.28)
с треугольной матрицей.
Упражнение 5.4. Построить алгоритм метода отражений, при котором для разло- жение матрицы в произведение ортогональной и верхней треугольной матриц доста- точно ?
2 3
n
3
действий умножения.

џ 6
Разностные уравнения
6.1 Разностные уравнения
Пусть y(n) = y n
 функция целочисленного аргумента n ? Z. Будем ее называть сеточной функцией. Обозначим через
(набла) оператор левой конечной разности
(разности назад), т.е.
y n
= y n
? y n?1
(6.1)
Степень оператора определим рекуррентным образом k
=
(
k?1
).
(6.2)
Пусть a(n), b(n), c(n), d(n) и f(n)  заданные сеточные функции. Рассмотрим урав- нение a(n)
3
y n
+ b(n)
2
y n
+ c(n) y n
+ d(n)y n
= f (n)
(6.3)
относительно сеточной функции y n
. Уравнение (6.3) называется разностным урав- нением. Разностные уравнения являются аналогами дифференциальных уравнений и в значительной степени повторяют свойства последних. Как и в случае диффе- ренциальных уравнений, важным является понятие порядка разностного уравнения.
Если a(n) = 0, то казалось бы естественным объявить порядком уравнения (6.3)
число три. Однако при таком определении порядка разностного уравнения нас ждут неприятности. Чтобы убедиться в этом, положим в (6.3) a(n) = 1, b(n) = 0, c(n) = ?3,
d(n) = 2
. В результате получим уравнение
3
y n
? 3 y n
+ 2y n
= f (n).
(6.4)
Принимая во внимание (6.1) и (6.2), находим, что y
n
= y n
? y n?1
,
3
y n
= y n
? 3y n?1
+ 3y n?2
? y n?3
,
а, подставляя эти выражения в (6.4), будем иметь y
n
? 3y n?1
+ 3y n?2
? y n?3
? 3y n
+ 3y n?1
+ 2y n
= f (n)
57

58
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
или
3y n?2
? y n?3
= f (n).
(6.5)
Вводя новый индекс m = n ? 2, уравнение (6.5) преобразуем к виду
3y m
? y m?1
= f (m + 2).
(6.6)
Это уравнение эквивалентно уравнению (6.4), и назвать его разностным уравнени- ем третьего порядка просто не поворачивается язык. И дело, конечно, не просто в названии. От удачно введенного определения зависит простота последующих утвер- ждений, использующих это определение. Поскольку запись (6.3) не содержит явным образом информации о числе, которым следовало бы определить порядок разностного уравнения, то будем разностное уравнение записывать в виде
?(n, y n
, y n?1
, . . . , y n?k
) = 0.
(6.7)
Определение 6.1. Уравнение (6.7) называется разностным уравнением.
Определение 6.2. Разностное уравнение (6.7), если оно явно зависит от y n
и от y
n?k
, назывется уравнением k-го порядка.
Определение 6.3. Разностное уравнение k-го порядка называется линейным, если оно линейно зависит от y n
, y n?1
, . . . , y n?k
Мы будем изучать только линейные разностные уравнения, которые будем запи- сывать в виде k
j=0
?
j
(n)y n?j
= f (n),
n ? Z.
(6.8)
Пока мы предполагаем, что уравнение (6.8) задано при всех n ? Z. Уравнение (6.8)
будет уравнением k-го порядка, если коэффициенты ?
0
(n)
и ?
k
(n)
не обращаются в нуль ни при одном n ? Z.
Определение 6.4. Сеточная функция y n
, n ? Z
называется решением уравнения
(6.8), если при подстановке ее в (6.8) последнее превращается в тождество.
Определение 6.5. Сеточная функция y n
, n ? Z
называется общим решением разностного уравнения (6.8), если в ней содержится любое решение (6.8).
Для того, чтобы определить какое-либо решение уравнения (6.8) (частное решение)
достаточно указать его значения в любых k последовательных точках, например,
n
0
, n
0
+ 1, . . . , n
0
+ k ? 1

6.2. УРАВНЕНИЯ ПЕРВОГО ПОРЯДКА
59 6.2 Линейные разностные уравнения первого поряд- ка
Эти уравнения имеют вид
?
0
(n)y n
+ ?
1
(n)y n?1
= f (n).
(6.9)
Поскольку ?
0
(n) = 0
, то на этот коэффициент уравнение можно поделить. Пусть
?
1
(n)/?
0
(n) = ?q n
, а f(n)/?
0
(n)
снова обозначим через f(n) = f n
. Тогда разностное уравнение первого порядка (6.9) можно переписать так y
n
= q n
y n?1
+ f n
(6.10)
Разрешить разностное уравнение  значит выразить y n
через известные величины.
Чтобы можно было решить (6.10), нужно задать начальное условие y
0
= a.
(6.11)
Используя теперь рекуррентные соотношения (6.10), можно последовательно опреде- лить y n
при всех последующих значениях n:
y
1
= q
1
y
0
+ f
1
= q
1
a + f
1
,
y
2
= q
2
y
1
+ f
2
= q
2
(q
1
a + f
1
) + f
2
= q
1
q
2
a + q
2
f
1
+ f
2
и т.д.
Часто бывает полезно иметь не рекуррентное соотношение для последовательно- го вычисления решения, а некоторую формулу, представляющую решение. Найдем представление решения уравнения (6.10). Для этого рассмотрим сначала отвечающее ему однородное уравнение y
n
= q n
y n?1
(6.12)
и найдем его решение. Имеем y
1
= q
1
y
0
,
y
2
= q
2
y
1
,
y n
= q n
y n?1
Перемножая последовательно полученные равенства и сокращая левую и правую части на y
1
y
2
. . . y n?1
, получим y
n
= q
1
. . . q n
y
0
= y
0
n j=1
q j
(6.13)
Величина y
0
есть начальное значение y n
и является произвольной постоянной. Реше- ние однородного уравнения (6.12) найдено.

60
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
Замечание 6.1. Напомним, что если линейное однородное дифференциальное урав- нение первого порядка записать в виде y = P (x)y, то его общее решение примет вид y(x) = c exp
?
?
?
x
0
P (?)d?
?
?
?
Обратимся теперь к неоднородному уравнению (6.10). Его решение будем искать,
используя решение однородного уравнения (6.12), методом вариации постоянной. Пусть
?
y n
=
n j=1
q j
(6.14)
Это  решение уравнения (6.12), а c
?
y n
 его общее решение. Заставим коэффициент c
зависеть от n и в таком виде будем искать решение уравнения (6.10)
y n
= c n
?
y n
(6.15)
Подставляя (6.15) в (6.10), получим c
n
?
y n
= q n
c n?1
?
y n?1
+ f n
Из (6.12)
?
y n
= q n
?
y n?1
и поэтому c
n
?
y n
= c n?1
?
y n
+ f n
,
т.е.
c n
= c n?1
+ f n
/
?
y n
Отсюда c
1
= c
0
+ f
1
?
y
1
,
c
2
= c
1
+ f
2
?
y
2
,
c n
= c n?1
+ f n
?
y n
,
Складывая эти соотношения, находим, что c
n
=
n k=1
f k
?
y k
+ c
0
,
а принимая во внимание (6.14), будем иметь c
n
=
n k=1
f k
k j=1
q
?1
j
+ c
0

6.3. УРАВНЕНИЯ K-ГО ПОРЯДКА С ПОСТОЯННЫМИ КОЭФФИЦИЕНТАМИ61
Подставляя это выражение в (6.15), получим общее решение неоднородного уравнения
(6.10)
y n
=
n j=1
q j
c +
n k=1
f k
k j=1
q
?1
j
(6.16)
Замечание 6.2. Напомним, что если линейное неоднородное дифференциальное урав- нение первого порядка записать в виде y = P (x)y+f(x), то его общее решение примет вид y(x) = exp
?
?
?
x
0
P (?)d?
?
?
?
?
?c +
x
0
f (?) exp
?
?
?
?
?
0
P (?)d?
?
?
?
d?
?
? .
Если коэффициент q n
=
const = q, то из (6.16) находим, что y
n
= q n
c +
n k=1
f k
q
?k
,
(6.17)
а если и f n
=
const = f, то при q = 1
y n
= q n
c + f n
k=1
q
?k
= q n
c + f q
?1
? q
?n+1 1 ? q
?1
= cq n
+ f
1 ? q n
1 ? q
(6.18)
6.3 Разностные уравнения k-го порядка с постоянны- ми коэффициентами
Если коэффициенты ?
j
(n)
из (6.8) не зависят от n, то мы имеем разностное уравнение с постоянными коэффициентами k
j=0
?
j y
n?j
= f (n),
n ? Z.
(6.19)
Решение отвечающего (6.19) однородного уравнения k
j=0
?
j y
n?j
= 0
(6.20)
можно искать в виде y
n
= q n
,
q = 0
ср. с y(x) = e
?x
,
(6.21)
где q = const = 0. Подставляя (6.21) в (6.20), получим q
n?k k
j=0
?
j q
k?j
= 0.

62
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
На q n?k можно сократить, в результате чего для отыскания q получим алгебраическое уравнение степени k
?
0
q k
+ ?
1
q k?1
+ · · · + ?
k?1
q + ?
k
= 0,
(6.22)
называемое характеристическим уравнением, отвечающим разностному уравнению
(6.20).
Характеристическое уравнение (6.22) имеет ровно k корней, включая кратные и комплексные. Обозначим их через q
1
, q
2
, . . . , q k
(6.23)
Очевидно, что сеточные функции q
n l
,
l = 1, . . . , k
(6.24)
являются решениями разностного уравнения (6.20).
Имеет место
Теорема 6.1. Если корни (6.23) характеристического уравнения (6.22) простые, то решения (6.24) разностного уравнения (6.20) линейно независимы, а общее решение этого уравнения имеет вид y
n
=
k l=1
c l
q n
l
Доказательство. Проведем доказательство линейной независимости (6.20) при k = 2
. Допустим противное, т.е. пусть c
1
q n
1
+ c
2
q n
2
? 0,
|c
1
| + |c
2
| = 0.
Но тогда и c
1
q n?1 1
+ c
2
q n?1 2
= 0.
Рассмотрим эти два тождества как систему уравнений для определения c
1
и c
2
. На- ходим, что определитель этой системы
? =
q n
1
q n
2
q n?1 1
q n?1 2
= (q
1
q
2
)
n?1
(q
1
? q
2
) = 0
и, следовательно, система имеет лишь тривиальное решение c
1
= c
2
= 0
. Это проти- воречит предположению, что и доказывает теорему.
Замечание 6.3. Если комплексное число q = |q|e i?
, ? = m?, m ? Z является корнем характеристичечкого уравнения (6.22), коэффициенты которого действительны, то

6.3. УРАВНЕНИЯ K-ГО ПОРЯДКА С ПОСТОЯННЫМИ КОЭФФИЦИЕНТАМИ63
число q = |q|e
?i?
, комплексно сопряженное к q, также является корнем характе- ристического уравнения (6.22), а наряду с комплексными решениями разностного уравнения (6.20)
q n
и q n
(6.25)
решениями указанного разностного уравнения будут и действительная и мнимая ча- сти решений (6.25), т.е.
|q|
n cos n?,
|q|
n sin n?.
(6.26)
Решения (6.26), как и (6.25), линейно независимы.
Пример 6.1. Найдем общее решение разностного уравнения y
n
? 2
ch ?y n?1
+ y n?2
= 0,
? = 0.
Характеристическое уравнение этого разностного уравнения имеет вид q
2
? 2
ch ?q + 1 = 0,
а его корни суть q
1,2
=
ch ? ±
?
ch
2
? ? 1 = e
±?
Эти корни различные, и поэтому y
n
= c
1
e
?n
+ c
2
e
??n
Теорема 6.2. Если q есть корень характеристического уравнения (6.22) кратности s
1
, то сеточная функция
P
s?1
(n)q n
,
где P
s?1
(n)
 произвольный многочлен, степень которого не выше s ? 1, является решением разностного уравнения (6.20). При этом решения n
l q
n
,
l = 0, . . . , s ? 1
линейно независимы.
Доказательство. Доказательство проведем для случая s = k = 2. Покажем сначала, что nq n
есть решение уравнения (6.20) при k = 2. Имеем
?
0
nq n
+ ?
1
(n ? 1)q n?1
+ ?
2
(n ? 2)q n?2
=
= q n?2
(?
0
q
2
+ ?
1
q + ?
2
)(n ? 2) + q(2?
0
q + ?
1
) =
= q n?1
(2?
0
q + ?
1
) = 0,
ибо 2?
0
q +?
1
= (?
0
q
2
+?
1
q +?
2
)
и обязано обращаться в нуль на корне кратности два.
Линейная независимость решений q n
и nq n
доказывается так же, как и в предыдущей теореме.

64
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
Теорема 6.3. Если правая часть f(n) неоднородного разностного уравнения (6.19)
имеет вид P
m
(n)
?
q n
, где P
m
(n)
 многочлен степени m, а
?
q является s-кратным,
s
0
, корнем характеристического уравнения (6.22), то уравнение (6.19) имеет решение вида y(n) = n s
Q
m
(n)
?
q n
(6.27)
Доказательство смотри, например, в [].
Упражнение 6.1. Найти общее решение неоднородного разностного уравнения y
n
? 2y n?1
+ y n?2
= n(1 + 2
n
).
Ответ.
y n
= c
1
+ c
2
n +
1 6
(n + 3)n
2
+ (n ? 2)2
n+2 6.4 Системы разностных уравнений
Рассмотрим систему двух линейных однородных разностных уравнений первого по- рядка с двумя неизвестными функциями u
n
= a
11
u n?1
+ a
12
v n?1
,
v n
= a
12
u n?1
+ a
22
v n?1
,
n ? Z,
(6.28)
где a ij
, i, j = 1, 2  постоянные. Введем в рассмотрение вектор-функцию y(n) = [u n
v n
]
T
и матрицу
A =
a
11
a
12
a
21
a
22
,
которую будем предполагать невырожденной, det A = 0. Используя введенные обо- значения, систему (6.28) можно переписать в векторном виде y(n) = Ay(n ? 1), det A = 0,
или A
?1
y(n) = y(n ? 1).
(6.29)
В записи (6.29) можно забыть, что y(n) был двумерный вектор, а A  матрица второго порядка. Будем мыслить систему (6.29) как систему m-го порядка. Решение этой системы будем искать в виде y(n) = ?q n
(6.30)
где q = const = 0, а ?  ненулевой m-мерный вектор. Подставляя (6.30) в (6.29),
получим
?q n
= A?q n?1
,

6.5. ЗАДАЧА НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ
65
а сокращая на q n?1
= 0
, приходим к системе
?q = A?.
Эта система однородна, и, чтобы у нее были нетривиальные решения, определитель ее матрицы должен быть равен нулю
A ? qI = 0.
(6.31)
Уравнение (6.31)  характеристическое уравнение системы (6.29)  является алгеб- раическим уравнением m-ой степени. Если все его корни q k
 собственные значения матрицы A  различны, то соответствующие собственные векторы ?
k линейно неза- висимы, и общее решение системы (6.29) принимает вид y(n) =
m k=1
c k
?
k q
n k
(6.32)
Матрица A может иметь полный набор линейно независимых собственных векторов и при наличии кратных корней характеристического уравнения (6.31). И в этом случае общее решение системы (6.29) имеет вид (6.32). Если же у канонической формы матрицы A имеются жордановы клетки, то для отыскания общего решения системы
(6.29) нужно поступать так же, как и в случае систем дифференциальных уравнений.
На этом мы останавливаться не будем.
Упражнение 6.2. Найти общее решение системы (6.29) с матрицей
A =
?1 ?1 1
?1
Ответ.
y(n) =
?
?
?
?
?
c
1
?
?
?
? sin
3?n
4
cos
3?n
4
?
?
? + c
2
?
?
?
cos
3?n
4
sin
3?n
4
?
?
?
?
?
?
?
?
q n/2 6.5 Разностная задача на собственные значения
До сих пор мы обсуждали вопросы отыскания общего решения разностного уравне- ния, которое зависит от k произвольных постоянных, если уравнение имеет порядок k
. Чтобы выделить единственное решение разностного уравнения, как и в случае дифференциального уравнения k-го порядка, нужно задать k линейно независимых начальных или граничных условий. Как и для дифференциального уравнения, для разностного уравнения можно поставить задачу на собственные значения.

66
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
Займемся этой задачей, решение которой понадобится нам при дальнейших иссле- дованиях. Пусть
?y n+1
+ 2y n
? y n?1
= ?y n
,
n = 1, . . . , N ? 1,
y
0
= y
N
= 0.
(6.33)
Требуется найти такие значения параметра ? (собственные значения), при которых однородная задача (6.33) имеет нетривиальные решения. Если исключить из первого уравнения (6.33) неизвестное y
0
, а из последнего уравнения  неизвестное y
N
, то полу- чим обычную алгебраическую задачу на собственные значения для трехдиагональной матрицы
A =
?
?
?
?
?
?
?
?
2 ?1
?1 2 ?1
?1 2 ?1
?1 2 ?1
?1 2
?
?
?
?
?
?
?
?
,
порядок которой равен N ? 1. Матрица A симметрична, и потому ее собственные зна- чения действительны, а собственные векторы, отвечающие различным собственным значениям, ортогональны в смысле скалярного произведения
(v, w) :=
N ?1
i=1
v i
w i
Найдем собственные значения и отвечающие им собственные функции задачи (6.33).
Перепишем уравнение (6.33) в виде
?y n+1
+ 2(1 ? ?/2)y n
? y n?1
= 0,
n = 1, . . . , N ? 1
и предположим, что
|1 ? ?/2|
1,
т.е. 0
?
4.
(6.34)
Тогда для некоторого ? можно положить
1 ? ?/2 = cos ?/N
(6.35)
и переписать уравнение так y
n+1
? 2 cos
?
N
y n
+ y n?1
= 0.
(6.36)
Характеристическое уравнение, отвечающее уравнению (6.36), есть q
2
? 2 cos
?
N
q + 1 = 0,
а его корни суть q
1,2
= cos
?
N
±
cos
2
?
N
? 1 = cos
?
N
± i sin
?
N
= e
±?i/N

6.5. ЗАДАЧА НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ
67
Тем самым,
y n
= c
1
sin
?n
N
+ c
2
cos
?n
N
(6.37)
есть общее решение уравнения (6.36). Потребуем, чтобы это решение удовлетворяло граничным условиям (6.34). Будем иметь y
0
= c
2
= 0,
y
N
= c
1
sin ? = 0.
(6.38)
Поскольку c
1
не может быть нулевым, то отсюда находим, что sin ? = 0 и, следова- тельно,
? = k?,
k ? Z.
(6.39)
Из (6.35) находим
? = ?
k
= 2 1 ? cos k?
N
= 4 sin
2
k?
2N
,
k ? Z,
(6.40)
а из (6.38) 
y n
= y
(k)
n
= c
1
sin k?n
N
(6.41)
При k = 0 решение y
(0)
n
? 0
и, следовательно, число ?
0
= 0
не является собственным значением. При k = N решение y
(N )
n
= c sin N? ? 0
, и ?
N
тоже не может быть собственным значением. Собственные значения
?
k
= 4 sin
2
k?
2N
,
k = 1, . . . , N ? 1
(6.42)
различны, ибо функция sin t при 0 < t < ?/2 является монотонной. Поскольку изучаемая задача эквивалентна алгебраической задаче на собственные значения для матрицы (N ? 1) порядка, то соотношения (6.42) задают все собственные значения задачи (6.33). Собственные функции y
(k)
n
= c
1
sin k?n
N
,
k = 1, . . . , N ? 1
(6.43)
ортогональны. Подсчитаем их нормы y
(k)
n
2
= y
(k)
n
, y
(k)
n
= c
2 1
N ?1
n=1
sin
2
k?n
N
=
= c
2 1
N ?1
n=1 1 ? cos
2k?n
N
2
= c
2 1
N ? 1 2
?
1 2
N ?1
n=1
cos
2k?n
N
Далее,
N ?1
n=1
cos
2k?n
N
=
Re
N ?1
n=1
e
2ik?
N
n
=
Re e
2ik?
? e
2ik?/N
e
2ik?/N
? 1
= ?1.

68
џ 6. РАЗНОСТНЫЕ УРАВНЕНИЯ
Тем самым,
y
(k)
n
2
= c
2 1
N/2 = 1
при c
1
=
2/N,
(6.44)
а ортонормированные собственные функции суть y
(k)
n
=
2/N sin k?n
N
,
k = 1, . . . , N ? 1.
(6.45)
Из (6.42) следует, что для всех собственных значений предположение (6.34) выполне- но. Поэтому в рассмотрении противоположного предположения смысла нет.
Упражнение 6.3. Решить следующую задачу на собственные значения
? y n+1
+ 2y n
? y n?1
,
n = 1, . . . , N ? 1,
? y
1
+ y
0
+
?
2
y
0
= 0,
y
N
? y
N ?1
+
?
2
y
N
= 0.
Ответ.
?
k
= 4 sin
2
k?
2N
,
k = 0, . . . , N,
y
(k)
n
=
2/N cos k?n
N

џ 7
Ортогональные многочлены
7.1 Общие ортогональные многочлены
Функцию ?(x) ? 0 будем называть весовой функцией на интервале (?1, 1), если на этом интервале она неотрицательна и интегрируема.
Пусть на (?1, 1) задана последовательность многочленов
P
0
(x), P
1
(x), . . . , P
n
(x), . . . ,
(7.1)
в которой каждый многочлен P
n
(x)
имеет степень n. Если для любых двух многочле- нов из этой последовательности выполняется условие
( P
m
, P
n
) :=
1
?1
?(x)P
m
(x)P
n
(x)dx = 0,
m = n,
то многочлены (7.1) называются ортогональными на (?1, 1) с весом ?(x).
Лемма 7.1. Если в системе из (n + 1) ортогональных многочленов
P
0
(x), P
1
(x), . . . , P
n
(x)
каждый многочлен P
k
(x)
имеет степень k, то всякий многочлен Q
n
(x)
степени n можно единственным образом представить в виде
Q
n
(x) = a
0
P
0
(x) + a
1
P
1
(x) + · · · + a n
P
n
(x).
(7.2)
Доказательство. Пусть ортогональные многочлены P
k
(x)
имеют вид
P
k
(x) = c
(k)
0
+ c
(k)
1
x + · · · + c
(k)
k x
k
,
c
(k)
k
= 0,
а многочлен
Q
n
(x) = c
0
+ c
1
x + · · · + c n
x n
69

70
џ 7. ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ
Подставляя эти представления в (7.2) и приравнивая коэффициенты при одинаковых степенях x k
, получим следующую систему линейных алгебраических уравнений для определения неизвестных коэффициентов a k
c
(0)
0
a
0
+ c
(1)
0
a
1
+ c
(2)
0
a
2
+ . . . + c
(n)
0
a n
= c
0
,
c
(1)
1
a
1
+ c
(2)
1
a
2
+ . . . + c
(n)
1
a n
= c
1
,
c
(n)
n a
n
= c n
По условию теоремы коэффициенты c
(k)
k рассматриваемой системы многочленов от- личны от нуля, и, следовательно, эта алгебраическая система имеет единственное решение. Лемма доказана.
Замечание 7.1. Умножая соотношение (7.2) на ?(x)P
k
(x)
и интегрируя результат по интервалу (?1, 1), легко находим, что a
k
=
( P
k
, Q
n
)
P
k
2
,
P
k
2
= ( P
k
, P
k
).
Лемма 7.2. Для всякой весовой функции ?(x) существует единственная последова- тельность ортонормированных многочленов {P
n
(x)}
, имеющих положительный коэффициент при старшей степени.
Доказательство. Обозначим коэффициент при старшей степени x многочлена
P
n
(x)
через µ
n
. Доказательство теоремы проведем методом полной математической индукции. Имеем, P
0
(x) = µ
0
> 0
и, следовательно,
( P
0
, P
0
) = µ
2 0
( 1 , 1 ) = 1.
Поэтому
µ
0
= 1/ ( 1 , 1 )
и многочлен P
0
(x)
определен.
Пусть определены ортонормированные многочлены
P
0
(x), P
1
(x), . . . , P
n?1
(x).
Определим многочлен P
n
(x)
. Будем его искать в виде P
n
(x) = µ
n x
n
+ Q
n?1
(x)
. В силу леммы 7.1 находим, что
P
n
(x) = µ
n x
n
+
n?1
k=0
a k
P
k
(x),
где числа µ
n и a k
подлежат определению. Умножая это соотношение скалярно на
P
m
(x)
,
m = 0, . . . , n ? 1
, находим, что
0 = µ
n
( x n
, P
m
(x) ) + a m
,
m = 0, . . . , n ? 1,

7.1. ОБЩИЕ ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ
71
т.е.
a m
= ?µ
n
( x n
, P
m
(x) )
и, следовательно,
P
n
(x) = µ
n x
n
?
n?1
k=0
( x n
, P
k
) x k
есть произвольный ортогональный многочлен степени n. Умножая его скалярно на самого себя и требуя нормированности, находим
1 = µ
2
n
?
? x n
?
n?1
k=1
( x n
, P
k
)x k
2
, 1
?
?
?
0
Отсюда определяем µ
n
> 0
. Лемма доказана.
Лемма 7.3. Если P
n
(x)
принадлежит совокупности ортогональных с весом ?(x)
многочленов, то для всякого многочлена Q
m
(x)
степени m < n
( Q
m
, P
n
(x) ) = 0,
m < n.
(7.3)
Доказательство. В силу леммы 7.1
Q
m
(x) = a
0
P
0
(x) + a
1
P
1
(x) + · · · + a m
P
m
(x).
Подставляя это представление в (7.3), получаем утверждение леммы.
Лемма 7.4. Если весовая функция ?(x) четная, то каждый ортогональный много- член P
n
(x)
содержит только те степени x, которые имеют одинаковую с номером n
четность, т.е.
P
n
(?x) ? (?1)
n
P
n
(x).
(7.4)
Доказательство.Пусть ?(x) = ?(?x) и
1
?1
?(x)P
n
(x)P
m
(x)dx = 0,
m = 0, . . . , n ? 1.
Заменой переменной интегрирования x = ?t эти условия приводятся к виду
1
?1
?(t)P
n
(?t)P
m
(?t)dt = 0,
m = 0, . . . , n ? 1,

72
џ 7. ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ
т.е. P
m
(?t)
тоже ортогональные многочлены. Но в силу леммы 7.2 любой ортогональ- ный многочлен определен с точностью до множителя, и поэтому
P
n
(?x) = c n
P
n
(x).
Отсюда в частности следует, что
(?1)
n a
n x
n
= c n
a n
x n
,
т.е. c n
= (?1)
n и поэтому
P
n
(?x) = (?1)
n
P
n
(x).
Лемма доказана.
Теорема 7.1. Для любых трех последовательных ортогональных многочленов спра- ведлива рекуррентная формула
?
n
P
n+1
(x) = (x ? ?
n
)P
n
(x) ? ?
n
P
n?1
(x).
(7.5)
Доказательство. Перепишем (7.5) в виде xP
n
(x) = ?
n
P
n+1
(x) + ?
n
P
n
(x) + ?
n
P
n?1
(x).
В левой части этого равенства стоит многочлен степени n + 1. В силу леммы 7.1 он может быть разложен по многочленам P
0
, . . . , P
n+1
xP
n
(x) =
n+1
k=0
c
(n+1)
k
P
k
(x),
(7.6)
где в силу замечания 7.1
c
(n+1)
k
=
( xP
n
, P
k
)
P
k
2
(7.7)
Но тогда с учетом (7.6)
c
(n+1)
k
=
1
P
k
2
( P
n
, xP
k
) =
1
P
k
2
P
n
,
k+1
j=0
c
(k+1)
j
P
j
=
=
1
P
k
2
k+1
j=0
c
(k+1)
j
( P
n
, P
j
) .
Отсюда следует, что если k + 1 < n, то c
(n+1)
k
= 0,
k + 1 < n и поэтому xP
n
(x) = c
(n+1)
n+1
P
n+1
(x) + c
(n+1)
n
P
n
(x) + c
(n+1)
n?1
P
n?1
(x),
(7.8)
т.е.
?
n
= c
(n+1)
n+1
,
?
n
= c
(n+1)
n
,
?
n
= c
(n+1)
n?1
(7.9)
Теорема доказана.

7.2. МНОГОЧЛЕНЫ ЧЕБЫШЕВА ПЕРВОГО РОДА
73
Теорема 7.2. Все нули ортогонального многочлена P
n
(x)
действительны, различны и расположены на интервале (?1, 1).
Доказательство. Достаточно показать, что многочлен P
n
(x)
на (?1, 1) меняет знак n раз. Допустим противное, т.е. что многочлен P
n
(x)
меняет знак только в точках
?
1
, ?
2
, . . . , ?
m
, где m < n. Тогда многочлен
Q
m
(x) = (x ? ?
1
)(x ? ?
2
) . . . (x ? ?
m
)
также меняет знак точно в этих же самых точках. Следовательно, произведение
P
n
(x)Q
m
(x) ? 0
сохраняет знак на (?1, 1) и следовательно
1
?1
?(x)P
n
(x)Q
m
(x)dx = 0,
что противоречит лемме 7.3. Противоречие снимается, если n = m. Теорема доказана.
Замечание 7.2. В силу доказанной теоремы для нулей x
(n)
k ортогонального много- члена P
n
(x)
имеют место неравенства
?1 < x
(n)
1
< x
(n)
2
< · · · < x
(n)
k
< · · · < x
(n)
n
< 1.
(7.10)
7.2 Многочлены Чебышева первого рода
Рассмотрим следующее однородное разностное уравнение второго порядка с постоян- ными коэффициентами y
n
? 2xy n?1
+ y n?2
= 0.
(7.11)
Здесь x  параметр. Поставим для (7.11) начальные условия y
0
= 1,
y
1
= x.
(7.12)
Тогда y
2
= 2x · x ? 1 = 2x
2
? 1,
y
3
= 2x(2x
2
? 1) ? x = 4x
3
? 3x,
y
4
= 8x
4
? 8x
2
+ 1, . . . .
(7.13)
Очевидно, что значение решения задачи (7.11), (7.12) в узле n есть многочлен от x степени n.
Найдем решение задачи (7.11), (7.12) в явном виде. Характеристическое уравнение разностного уравнения (7.11) имеет вид q
2
? 2xq + 1 = 0,

74
џ 7. ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ
а его корни суть q
1
= q = x +
?
x
2
? 1
и q
2
= 1/q.
(7.14)
Поэтому общее решение уравнения (7.11) есть y
n
= c
1
q n
+ c
2
q
?n
Полагая здесь n = 0 и n = 1 и принимая во внимание начальные условия (7.12),
находим, что y
0
= c
1
+ c
2
= 1,
y
1
= c
1
q + c
2
q
?1
= x.
(7.15)
Преобразем второе из уравнений (7.15) с учетом (7.14) и первого уравнения,
y
1
= c
1
(x +
?
x
2
? 1) + c
2
(x ?
?
x
2
? 1) =
= (c
1
+ c
2
)x + (c
1
? c
2
)
?
x
2
? 1 = x + (c
1
? c
2
)
?
x
2
? 1 = x.
Отсюда следует, что c
1
= c
2
, а с учетом (7.15) находим, что c
1
= c
2
= 1/2,
и поэтому y
n
=
q n
+ q
?n
2
(7.16)
есть решение задачи (7.11), (7.12).
Как было замечено раньше, это есть многочлен от x степени n. Пусть |x| < 1.
Тогда в силу (7.14)
q = x + i
?
1 ? x
2
и, следовательно, |q| = 1. Пусть q = e i?
. Тогда x = cos ?,
? = arccos x,
(7.17)
y n
=
e in?
+ e
?in?
2
= cos n? = cos[n arccos x].
Определение 7.1. Алгебраические многочлены
T
n
(x) = cos[n arccos x],
|x| < 1,
n = 0, 1, . . .
(7.18)
называются многочленами Чебышева первого рода.
Они принадлежат к семейству ортогональных многочленов. Определим весовую функцию ?(x), при которой многочлены T
n
(x)
будут ортогональными. Из (7.17) сле- дует, что d? = ?
dx
?
1 ? x
2
,
x = 1
при ? = 0 и x = ?1 при ? = ?.

7.3. СВОЙСТВА МНОГОЧЛЕНОВ ЧЕБЫШЕВА
75
Принимая теперь во внимание (7.18), находим, что при m = n
0 =
?
0
cos m? cos n? d? =
1
?1 1
?
1 ? x
2
T
m
(x)T
n
(x)dx.
Тем самым
?(x) = (1 ? x
2
)
?1/2
(7.19)
7.3 Свойства многочленов Чебышева
1
?
. При четном n многочлен T
n
(x)
является четной функцией x, а при нечетном n 
нечетной.
Доказательство следует из (7.19) и леммы 7.4.
2
?
. Коэффициент при старшей степени многочлена T
n
(x)
для n
1
равен 2
n?1
, т.е.
µ
n
= 2
n?1
, а
T
n
(x) = 2
n?1
x n
+ . . . .
Доказательство. См. рекуррентную формулу (7.11).
3
?
. Нули многочлена T
n
(x)
расположены в точках x
k
= ? cos
(2k ? 1)?
2n
,
k = 1, 2, . . . , n.
(7.20)
Доказательство. Из (7.18) находим, что n arccos x k
= ?
?
2
+ k? =
(2k ? 1)?
2
или arccos x k
=
(2k ? 1)?
2n
,
т.е.
x k
= cos
(2k ? 1)?
2n
,
k = 1, n.
Так как функции T
n
(x)
являются либо четными, либо нечетными, то нули T
n
(x)
расположены симметрично относительно начала координат x
n+1?k
= ?x k
= ? cos
(2k ? 1)?
2n
Перенумеровывая нули в обратном порядке, приходим к (7.20).
4
?
. max
[?1,1]
|T
n
(x)| = 1
, причем
T
n
(x m
) = (?1)
m
,

76
џ 7. ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ
где x
m
= ? cos m?
n
,
m = 0, . . . , n.
(7.21)
Доказательство очевидно.
5
?
. Среди всех многочленов степени n с единичным коэффициентом при старшей степени многочлен
T
n
(x) =
1 2
n?1
T
n
(x),
n
1
на [?1, 1] имеет наименьшее значение максимума модуля.
Доказательство. Допустим противное, т.е. допустим существование такого мно- гочлена P
n
(x) = x n
+ . . .
, что max
[?1,1]
|P
n
(x)| < max
[?1,1]
|T
n
(x)|.
(7.22)
Тогда T
n
(x) ? P
n
(x) ? 0
и это есть многочлен степени не выше (n ? 1). Более того, в (n + 1) точке (7.21) этот многочлен принимает отличные от нуля значения с чередующимися знаками.
?1
P
n
T
x m?1
Ч
Ч
Ч
x m
x m+1
x m
Рис. 1
Но это означает, что алгебраический многочлен T
n
(x) ? P
n
(x)
степени меньшей n обращается в нуль по крайней мере в n точках, что невозможно.
Замечание 7.3. Можно доказать, что если P
n
(x) = x n
+ . . . ,
n
1
, и max
[?1,1]
|P
n
(x)| = 2
?n+1
,
то P
n
(x) ? T
n
(x) = 2
?n?1
T
n
(x)

7.4. МНОГОЧЛЕНЫ ЛЕЖАНДРА
77
Благодаря свойству 5
?
многочлены Чебышева T
n
(x)
называются многочленами,
наименее уклоняющимися от нуля.
6
?
. Если x > 1, то
T
n
(x) =
ch n Arch x,
где
Arch x = ln(x +
?
x
2
? 1).
Доказательство. В силу (7.16)
T
n
(x) =
q n
+ q
?n
2
=
e n ln q
+ e
?n ln q
2
=
ch n ln q = ch n ln(x +
?
x
2
? 1).
Замечание 7.4. Arch x  обратная функция к ch x.
Упражнение 7.1. Доказать, что ch Arch x = x, Arch ch x = x.
7.4 Многочлены Лежандра
Многочленами Лежандра называются многочлены, которые ортогональны друг другу на [?1, 1] с весом ? ? 1. Обозначаются они через P
n
(x)
1
?1
P
m
(x)P
n
(x)dx = 0,
m = n.
Если P
0
(x) ? 1
, то P
1
(x) ? x
Трехточечное рекуррентное соотношение для многочленов Лежандра имеет вид
(n + 1)P
n+1
(x) ? (2n + 1)xP
n
(x) + nP
n?1
(x) = 0
и следовательно
P
2
(x) =
1 2
(3x
2
? 1),
P
3
(x) =
1 2
(5x
3
? 3x),
P
4
(x) =
1 8
(35x
4
? 30x
2
+ 3),
P
5
(x) =
1 8
(63x
5
? 70x
3
+ 15x)
и т.д.

78
џ 7. ОРТОГОНАЛЬНЫЕ МНОГОЧЛЕНЫ

џ 8
Итерационные методы
В предыдущих лекциях для системы линейных алгебраических уравнений
Ax = b
(8.1)
с квадратной невырожденной матрицей были рассмотрены четыре прямых метода отыскания решения:
а) метод Гаусса (LU-разложение, треугольное разложение) и его модификация с выбором ведущего элемента,
б) метод Холецкого, применяемый в случае симметричной положительно опреде- ленной матрицы,
в) метод вращений,
г) метод отражений.
Все эти методы позволяют в принципе (при отсутствии ошибок округления) найти точное решение за конечное число действий. Это число действий было оценено нами величиной O(n
3
)
, где n  порядок системы. Если матрица A системы имеет ленточную структуру с полушириной ленты p много меньшей n, то ленточные варианты первых двух методов позволяют найти точное решение с меньшей, чем O(n
3
)
, затратой дей- ствий.
В этой лекции мы рассмотрим другой класс методов решения системы (8.1) 
итерационных. Эти методы, как правило, если и позволяют найти точное решение системы (8.1), то только как предел при стремлении числа итераций (а, следовательно,
и действий) к бесконечности. Однако для широкого класса задач, встречающихся в приложениях, те или иные итерационные методы могут оказаться предпочтительнее с точки зрения используемых трудозатрат, чем описанные прямые.
79

80
џ 8. ИТЕРАЦИОННЫЕ МЕТОДЫ
8.1 Одношаговые итерационные методы
Из курса "Введение в численные методы" известно, что многие одношаговые итера- ционные методы могут быть записаны в так называемой канонической форме
B
x k+1
? x k
?
+ Ax k
= b,
k = 0, 1, . . . ,
(8.2)
где B  некоторая матрица, определяющая итерационный метод, а ?  итерационный параметр. В частности, в виде (8.2) могут быть записаны метод Якоби, метод Гаусса-
Зейделя, метод последовательной верхней релаксации, метод простых итераций. Пред- положим, что
A = A
T
> 0,
B = B
T
> 0.
(8.3)
При этих предположениях в курсе "Введение в численные методы"доказано, что если
B >
?
2
A,
(8.4)
т.е. если для любого ненулевого вектора x справедливо
( Bx , x ) >
?
2
( Ax , x )
, то итерационный метод (8.2) является сходящимся.
Для метода простых итераций, который имеет вид (8.2) с
B = I
, кроме того, дана и оценка скорости сходимости. Именно, если
? = 2/(?
1
+ ?
n
),
(8.5)
где ?
1
и ?
n
, соответственно, наименьшее и наибольшее собственные значения матрицы
A
, то для погрешности итераций справедлива оценка x ? x k
q x ? x k?1
,
q =
?
n
? ?
1
?
n
+ ?
1
=
1 ? ?
1
/?
n
1 + ?
1
/?
n
< 1,
(8.6)
где ·  евклидова длина.
Из (8.6) следует, что если ?
n
?
1
, то число q очень близко к единице, а скорость сходимости итераций очень низкая. Но при указанном выборе нормы вектора
A = max i
?
i
(A) = ?
n
,
A
?1
= max i
?
i
(A
?1
) = ?
?1 1
и
?
n
/?
1
= A
A
?1
=
cond A := ?(A) := ?.
Таким образом, если матрица A плохо обусловлена, а это типичная ситуация, то метод простых итераций будет сходиться очень медленно.