Методическое пособие для практических и лабораторных работ по курсу «Численные методы»

128
Рис.2
Рис. 3.

Рис. 4.
Рис. 5.
Обозначения: d - точные входные данные, d% - входные данные с малой погрешностью, т.е.
1
d
d
ε
− ≤
%
= , x - точное решение с точными входными данными,
в
x - вычисленное решение с точными входными данными, x
% - точное решение, отвечающее входным данным d% ,
в
x
% - вычисленное решение с входными данными d% . На рисунках: слева – области входных данных, а справа – области решений.
Рисунок 2 соответствует случаю, когда и задача и численный метод хорошо обусловлены. В этом случае все четыре результата лежат близко друг к другу.
Рисунок 3 соответствует случаю, когда задача хорошо обусловлена, а вычисления плохо обусловлены. Здесь вычисленные
в
x и
в
x
% далеко
«убежали» от точных x и
x
%
Рисунок 4 соответствует случаю, когда плохо обусловлена задача, а численный метод «хороший». Из-за малых погрешностей во входных данных точные решения x и
x
%
далеко «убежали» друг от друга.
На рисунке 5 представлен случай, когда и задача и вычисления плохо обусловлены. Здесь уже все четыре переменных x , x
% ,
в
x и
в
x
% далеко
«разбежались» друг от друга.
Замечание. Иногда причиной появления плохо обусловленных задач становится отсутствие у пользователей ЭВМ элементарного представления об их существовании. В последнее время, к сожалению, получила развитие опасная тенденция пренебрежительного отношения к математическим знаниям вообще и к знанию вычислительной математики в частности.

4.3.
Вычисление определителей

Пусть дана матрица

129












=
nn
n
n
n
n
a
a
a
a
a
a
a
a
a
A
...
...
...
2 1
2 22 21 1
12 11
и требуется вычислить ее определитель. Для вычисления определителей матриц можно применять алгоритмы точных методов решения систем линейных алгебраических уравнений
b
x
A
r r
=
. Например, выполняемые преобразования прямого хода в методе исключения Гаусса, приводящие матрицу А системы к треугольному виду, таковы, что они не изменяют определителя матрицы А. Учитывая, что определитель треугольной матрицы равен произведению диагональных элементов, имеем
)
1
(
)
1
(
22 11
)
1
(
)
1
(
2
)
1
(
2
)
1
(
22 11
)
1
(
)
1
(
2
)
1
(
2
)
1
(
22 1
12 11 2
1 2
22 21 1
12 11
...
...
...
0
...
...
0
...
...
...
...
det
−
⋅
⋅
=
=
=
=
n
nn
nn
n
n
nn
n
n
n
nn
n
n
n
n
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
A

Таким образом, определитель матрицы равен произведению всех ведущих элементов при ее преобразовании методом Гаусса.
Если известно
LU
– разложение матрицы:
LU
A
=
, где
L – нижняя, а
U
– верхняя треугольные матрицы, то, очевидно,
∏
=
=
⋅
=
n
i
ii
ii
u
l
U
L
A
1
det det det

4.4.
Обращение матриц

Определение 6. Обратной к матрице A называется такая матрица
1
−
A , для которой
E
A
A
A
A
=
⋅
=
⋅
−
−
1 1
, где
E – единичная матрица:












=
1 0
0
...
...
0 1
0 0
0 1
E

130
Определение 7. Квадратная матрица A называется неособенной или невырожденной, если ее определитель
A
det отличен от нуля.
Всякая неособенная матрица имеет обратную матрицу. Пусть дана неособенная матрица












=
nn
n
n
n
n
a
a
a
a
a
a
a
a
a
A
...
...
...
2 1
2 22 21 1
12 11

Требуется найти
1
−
A . Обозначим элементы обратной матрицы
)
(
1
ij
x
A
=
−
,
n
j
i
,...,
1
,
=
. По определению
E
AA
=
−
1
. Перепишем данное равенство в следующем виде:












=












⋅












1 0
0
...
...
0 1
0 0
0 1
...
...
...
...
...
...
...
2 1
2 22 21 1
12 11 2
1 2
22 21 1
12 11
nn
n
n
n
n
nn
n
n
n
n
x
x
x
x
x
x
x
x
x
a
a
a
a
a
a
a
a
a
(4.1)
Обозначим столбцы обратной матрицы
1
−
A как векторы
j
xr , а столбцы единичной матрицы E через
j
er . Из равенства (4.1) получаем
j
j
e
x
A
r r
=
,
n
j
,...,
1
=
.
(4.2)

Следовательно, чтобы найти элементы обратной матрицы нужно, найти решения n систем линейных алгебраических уравнений (4.2). Так как эти системы имеют одну и ту же матрицу коэффициентов A, то их можно решать одновременно (делая соответствующие преобразования сразу над n правыми частями).
4.5.
Применение метода итераций для уточнения элементов обратной
матрицы

Пусть для неособенной матрицы A найдены приближенные значения элементов обратной матрицы
1
−
A . Обозначим матрицу с такими элементами

131
через
1 0
−
≈
A
D
. Считаем, что полученная точность приближения элементов матрицы
1
−
A нас не удовлетворяет. Для уточнения элементов обратной матрицы строим следующий итерационный процесс:
1 1
−
−
−
=
k
k
AD
E
F
,
,...
2
,
1
=
k
,
(4.3)
)
(
1 1
−
−
+
=
k
k
k
F
E
D
D
,
,...
2
,
1
=
k
(4.4)

Доказано, что итерации сходятся, если начальная матрица
0
D достаточно близка к искомой матрице
1
−
A .
Матрица
1
−
k
F на каждом шаге характеризует в некотором смысле степень близости матрицы
1
−
k
D к матрице
1
−
A .
Обычно итерации продолжают до тех пор, пока элементы матрицы
k
F по модулю не станут меньше заданного числа
ε
, и тогда приближенно полагают
k
D
A
≈
−
1
Указанный итерационный процесс оказывается весьма полезным, так как точные методы вычисления обратной матрицы часто приводят к заметным погрешностям, вызванным неизбежными ошибками округления и большим количеством арифметических операций при расчете.

132
Лекция 5. Ортогональные преобразования. Матрицы вращения и отражения.
QR- и
HR-разложения матриц.
Метод ортогонализации. Метод отражений.
5.1. Ортогональные преобразования. Матрицы вращения и
отражения. QR- и HR-разложения матриц.

Большинство методов исключения решения системы линейных алгебраических уравнений
b
x
A
r r
=
неособыми преобразованиями приводят исходную матрицу к верхней треугольной матрице. Решение задачи (3.1) сводится к решению некоторой эквивалентной задачи
β
r r
=
x
U
(5.1) с верхней треугольной матрицей
U
. Решение последней задачи не представляет затруднений.
Переход от задачи (3.1) к задаче (5.1) осуществляется путем умножения матрицы A слева на последовательность некоторых матриц
n
B
B ,...,
1
, т.е. на матрицу
1 1
... B
B
B
B
n
n
⋅
⋅
⋅
=
−
. Итак
BA
U
=
По свойствам числа обусловленности матрицы
)
(
)
(
)
(
)
(
A
cond
B
cond
BA
cond
U
cond
⋅
≤
=
, т.е. число обусловленности произведения матриц, не меньше чисел обусловленности сомножителей. Тем самым при преобразованиях, осуществляемых умножением на последовательность матриц, обусловленность системы линейных алгебраических уравнений ухудшается, что может привести к существенному росту погрешности решения.
Такого нежелательного явления можно избежать, если осуществлять приведение матрицы A к верхнему треугольному виду с помощью ортогональных преобразований.
Определение 1. Матрица Q с вещественными коэффициентами
ij
q называется ортогональной, если
E
Q
Q
T
=
,
1
−
=
Q
Q
T
(5.2)
Отметим следующие свойства ортогональных матриц

133 2
2
x
x
Q
r r
=
,
2 2
A
QA
=
(5.3)
Действительно
2 1
2 2
)
,
(
)
,
(
)
,
(
x
x
x
x
x
Q
Q
x
x
Q
x
Q
x
Q
n
i
i
T
r r
r r
r r
r r
=
=
=
=
=
∑
=
и
2 1
1 1
2
)
(
max
)
(
max
)
)
((
max
A
A
A
QA
Q
A
QA
QA
QA
T
j
n
j
T
T
j
n
j
T
j
n
j
=
=
=
=
≤
≤
≤
≤
≤
≤
λ
λ
λ

Отсюда
)
(
)
(
)
(
2 1
2 2
1 2
A
cond
A
A
QA
QA
QA
cond
=
⋅
=
=
−
−
(5.4)
Итак, если задача (3.1) сведена к задаче (5.1) с помощью ортогонального преобразования, то
)
(
)
(
A
cond
U
cond
=
, т.е. обусловленность исходной задачи не изменилась.
Рассмотрим матрицы вращения
)
,
( l
k
Q
,
k
l
>
. Ненулевые элементы
ij
q матрицы
)
,
( l
k
Q
определены следующим образом






=
=
=
=
−
=
=
=
=
≠
≠
=
=
,
,
sin
,
,
,
sin
,
,
cos
,
,
,
1
k
j
l
i
если
l
j
k
i
если
l
j
i
или
k
j
i
если
l
i
k
i
и
j
i
если
q
ij
ϕ
ϕ
ϕ

Все остальные элементы матрицы
)
,
( l
k
Q
равны нулю. Нетрудно проверить, что это ортогональные матрицы.
Умножим слева (3.1) на
)
2
,
1
(
Q
:
b
Q
x
A
Q
r r
)
2
,
1
(
)
2
,
1
(
=

134
Обозначим
)
2
,
1
(
)
2
,
1
(
A
A
Q
=
и элементы матрицы
)
2
,
1
(
A
через
ij
a
)
2
,
1
(
Матрица
)
2
,
1
(
A
отличается от матрицы A только первыми двумя строками.
Пусть
j
j
j
a
a
a
2 12 1
12 1
sin cos
)
2
,
1
(
ϕ
ϕ
−
=
,
j
j
j
a
a
a
2 12 1
12 2
cos sin
)
2
,
1
(
ϕ
ϕ
+
=
,
ij
j
a
a
=
2
)
2
,
1
(
при
2
>
i

Пусть хотя бы один из коэффициентов
11
a ,
21
a отличен от нуля.
Положив
2 21 2
11 21 12
sin
a
a
a
+
−
=
ϕ
,
2 21 2
11 11 12
cos
a
a
a
+
=
ϕ
, получаем
0
)
2
,
1
(
21
=
a
,
2 21 2
11 11
)
2
,
1
(
a
a
a
+
=

На этом шаге мы занулили коэффициент при
1
x во втором уравнении системы
b
x
A
r r
=
. Теперь умножаем слева преобразованную систему на
)
3
,
1
(
Q
, получаем
b
Q
Q
x
A
Q
Q
r r
)
2
,
1
(
)
3
,
1
(
)
2
,
1
(
)
3
,
1
(
=

Полагая
2 31 2
11 31 13
)
2
,
1
(
)
2
,
1
(
)
2
,
1
(
sin
a
a
a
+
−
=
ϕ
,
2 31 2
11 11 13
)
2
,
1
(
)
2
,
1
(
)
2
,
1
(
cos
a
a
a
+
=
ϕ
, получаем
0
)
3
,
1
(
31
=
a
,
0
)
2
,
1
(
)
2
,
1
(
)
3
,
1
(
2 31 2
11 11
>
+
=
a
a
a
, т.е мы занулили коэффициент при
1
x в третьем уравнении. Продолжая таким образом, мы в итоге придем к системе с верхней треугольной матрицей.
Сравнивания проведенные преобразования с рассмотренным ранее методом вращений, видим, что мы сейчас просто по иному изложили этот метод. Итак, в методе вращений исходную матрицу приводят к верхней

135
треугольной последовательным домножением на ортогональные матрицы вращений, не ухудшая обусловленность системы. Итог таких преобразований можно сформулировать в виде теоремы.
Теорема 1. Пусть
0
det
≠
A
. Тогда существует ортогональная
матрица Q такая, что

U
QA
=

(5.5)

При этом все диагональные элементы верхней треугольной матрицы
U , за исключением, возможно,
nn
u , положительны.
Из (5.5) получаем, что
U
Q
A
1
−
=
. Так как матрица
Q – ортогональна, то предыдущую теорему можно переформулировать так.
Теорема 2. Пусть
0
det
≠
A
. Матрица A разлагается в произведение

QR
A
=
,
(5.6)

где Q - ортогональная, а R – верхняя треугольная матрица.
Пусть матрица A системы (3.1) – ортогональная, тогда
b
A
x
T
r r
=
.
(5.7)
Отсюда понятна идея методов ортогонализации: нужно построить такое преобразование, которое ортогонализирует строки исходной матрицы.
Пусть
i
ar – строки матрицы A . Если
0
det
≠
A
, то система векторов
i
ar является линейно независимой и, следовательно, образует базис в
n
R .
Следовательно отыскание нужного преобразования сводится к задаче об ортогонализации базиса. Искомый ортогональный базис можно строить с помощью алгоритма Грамма-Шмидта. Изложим метод ортогонализации решения СЛАУ.

5.2. Метод ортогонализации
Систему уравнений
b
x
A
r r
=
запишем в виде
0
)
,
(
1
=
y
a r r
,
0
)
,
(
2
=
y
a r r
,…,
0
)
,
(
=
y
a
n
r r
,
(5.8) где
)
,
,...,
,
(
1
,
2 1
+
=
n
k
kn
k
k
k
a
a
a
a
ar
,
n
k
,...,
1
=
;
)
(
1
,
k
n
k
b
a
−
=
+
,
)
1
,
,...,
,
(
)
(
)
2
(
)
1
(
n
x
x
x
y
=
r

136
Из (5.8) видно, что для того, чтобы найти решение исходной системы, нужно найти вектор yr , который ортогонален линейно независимым векторам
n
a
a
r r
,...,
1
и имеет последнюю координату, равную 1.
Ортогональность векторов
yr и
n
a
a
r r
,...,
1
влечет за собой ортогональность вектора yr ко всему подпространству
n
P , натянутому на вектора
n
a
a
r r
,...,
1
, и, следовательно, всякому базису этого подпространства.
Верно и обратное. Поэтому для нахождения решения исходной системы достаточно построить какой-либо ортогональный базис подпространства
n
P и найти вектор zr , ортогональный этому базису.
Если
)
1
(
+
n
z
– последняя координата вектора
zr , то вектор
)
1
(
/
+
=
n
z
z
y
r r
с отброшенной последней координатой и есть решение исходной системы, т.е.
)
1
(
)
(
)
(
/
+
=
n
i
i
z
z
x
,
n
i
,...,
1
=
Итак, алгоритм построения решения следующий.
К системе линейно независимых векторов
n
a
a
r r
,...,
1
добавляем еще один линейно независимый вектор
1
+
n
ar
:
{ )
1
,
0
,...
0
(
1
n
n
a
=
+
r
Строим систему ортонормированных векторов
1 1
,...,
+
n
b
b
r r
таких, что для любого k (
1 1
+
≤
≤
n
k
) последовательность векторов
k
b
b
b
r r
r
,...,
,
2 1
является ортонормированным базисом подпространства
k
P , натянутого на
k
a
a
r r
,...,
1
. В этом случае вектор
1
+
n
b
r ортогонален к подпространству
n
P , натянутому на вектора
n
a
a
r r
,...,
1
, и искомое решение будет иметь вид
)
1
(
1
)
(
1
)
(
/
+
+
+
=
n
n
i
n
i
b
b
x
,
n
i
,...,
1
=
,
(5.9)
где
)
(
1
k
n
b
+
(
1
,...,
2
,
1
+
=
n
k
) – компоненты вектора
1
+
n
b
r
Сначала строим ортогональный базис
{ }
i
ur , а через него ортонормированный базис
{ }
i
b
r
, по следующим формулам:
1 1
a
u
r r
=
,
1 1
1
u
u
b
r r
r
=
,
)
,
(
1 1
2 1
u
u
u
r r
r
=
,
∑
=
+
+
+
−
=
k
i
i
i
k
k
k
b
b
a
a
u
1 1
1 1
)
,
(
r r
r r
r
,
1 1
1
+
+
+
=
k
k
k
u
u
b
r r
r
,
( )
n
k
,...,
1
=
.
(5.10)

Построив вектор
1
+
n
b
r
, решение вычисляем по формулам (5.9).

137
Отметим, что этот метод требует, как и все точные методы,
)
(
3
n
O
операций. Метод простой, однако неустойчивый, если в матрице A строки близки к линейно зависимым.

5.3. Метод отражений

Данный метод основан на переходе от заданной системы
b
x
A
r r
=
к новой системе
b
U
x
UA
r r
=
такой, что система
d
Rx
=
, где
UA
R
=
и
Ub
d
=
, решается проще, чем исходная. При выборе матрицы
U
нужно учитывать, по крайней мере, следующие два фактора. Во-первых, ее вычисление не должно быть чересчур сложным и трудоемким. Во-вторых, умножение на матрицу
U
не должно «портить» матрицу A (мера обусловленности матрицы A не должна значительно увеличиваться).
Такое приведение матрицы A к верхней треугольной матрице R можно осуществить с помощью последовательных ортогональных преобразований отражения.
В основе преобразования матрицы A к виду R (определяемому условием
0
=
ij
r
при
i
j
<
) лежит преобразование Хаусхолдера, или,
преобразование отражения,
T
E
U
ω
ω
r r
2
−
=
,
(5.11)
где
ω
r
– некоторый вектор-столбец единичной длины,
1
)
,
(
=
ω
ω
r r
. Под
T
ω
ω
r r
здесь понимается матрица, являющаяся произведением вектора-столбца
ω
r на вектор-строку
T
ω
r
, т.е.
)
(
ij
T
ω
ω
ω
=
r r
, где
j
i
ij
ω
ω
ω =
Из определения следует, что
T
ω
ω
r r
– симметричная матрица.
Воспользовавшись тем, что
1
)
,
(
=
=
ω
ω
ω
ω
r r
r r
T
,
(5.12)
убедимся, что
T
U
U
=
:
E
E
E
E
UU
T
T
T
T
T
T
T
T
=
+
−
−
=
−
−
=
ω
ω
ω
ω
ω
ω
ω
ω
ω
ω
ω
ω
r r
r r
r r
r r
r r
r r
4 2
2
)
2
)(
2
(

Таким образом, матрица
U
– ортогональная и симметричная. Значит, матрицы A и B , связанные соотношением
UAU
B
=
(
)
1
−
=
=
UAU
UAU
T
, являются подобными, этот факт используется для решения задач на собственные значения.
Теперь нетрудно догадаться, как нужно распорядиться свободой задания элементов векторов
ω
r при построении матриц отражения, чтобы за

138
конечное число шагов преобразований Хаусхолдера произвольно заданную матрицу A привести к верхней треугольной матрице R .
А именно, можно показать, что начатым с
A
R
=
1
процессом
m
m
m
R
U
R
=
+
1
,
1
,
2
,...,
2
,
1
−
−
=
n
n
m
,
(5.13)
где
T
m
m
m
E
U
ω
ω
r r
2
−
=
, данная
n
n
×
матрица
A за
1
−
n
шаг будет приведена к треугольной матрице R , если задающие матрицы Хаусхолдера
m
U векторы
m
ω
r по данной матрице A строить следующим образом.
При
1
=
m
вектор
1
ω
r определяется равенством
)
,...,
,
(
1 21 1
11 1
1
n
T
a
a
s
a
−
= µ
ω
r
,
(5.14)
где
∑
=
⋅
−
=
n
i
i
a
a
sign
s
1 2
1 11 1
)
(
,
)
(
2 1
11 1
1 1
a
s
s
−
=
µ
(5.15)

Такое задание
1
ω
r обеспечивает ортогональность симметричной матрицы
T
E
U
1 1
1 2
ω
ω
r r
−
=
и одновременное получение с ее помощью нужных
1
−
n
нулей в первом столбце матрицы
1 1
2
R
U
R
=
(
)
A
U
1
=
. Элементы матрицы
2
R обозначим через
)
2
(
ij
r . Далее при
2
=
m
строим вектор
2
ω
r следующим образом:
)
,...,
,
,
0
(
)
2
(
2
)
2
(
32 2
)
2
(
22 2
2
n
T
r
r
s
r
−
= µ
ω
r
,
(5.16) где
∑
=
⋅
−
=
n
i
i
r
r
sign
s
2 2
)
2
(
2
)
2
(
22 2
)
(
)
(
,
)
(
2 1
)
2
(
22 2
2 2
r
s
s
−
=
µ
(5.17)

Такое задание
T
2
ω
r обеспечивает ортогональность симметричной матрицы
T
E
U
2 2
2 2
ω
ω
r r
−
=
и одновременное получение с ее помощью нужных
2
−
n
нулей во втором столбце матрицы
2 2
3
R
U
R
=
Вектор
3
ω
r по матрице
3
R строится совершенно аналогично, только фиксируется нулевым не одна, а две первые его координаты, и определяющую роль играет теперь не второй, а третий столбец матрицы
3
R и

139
его третий элемент. При этом у матрицы
3 3
4
R
U
R
=
окажется
3
−
n
нулевых элемента в третьем столбце и сохранятся полученные на предыдущем шаге нули в первом и втором столбцах.
Этот процесс очевидным образом может быть продолжен до исчерпания и без особого труда может быть описан общими формулами типа (5.13) – (5.17), для чего нужно лишь ввести обозначения для элементов последовательности матриц
m
R .
Матрицу отражений часто обозначают буквой Н. Результат проделанных преобразований можно сформулировать в виде теоремы.

Скачать 1.29 Mb.

Поделитесь с Вашими друзьями: