Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление

• Название: Конспект лекций. Санкт петербург 2011 2 Оглавление
• Анкор: Вычислительная математика лекции.pdf
• Дата: 01.04.2017
• Размер: 3.86 Mb.
• Имя файла: Вычислительная математика лекции.pdf
• Тип: Конспект; #4404
• страница: 16 из 19

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

169
Назовем итерационным циклом последовательность из
(
1)
2
n n

шагов (количество циклов k=1,…) , состоящих в переборе элементов a
lk
в нижней поддиагональной части матрицы, то есть
1,...,
1,...,
k
n
l
k
n

 
На i – ом шаге цикла осуществляется подобное преобразование
A
i
→A
i+1
=T
kl
A
i
T
kl
T
(T
kl
T
=T
kl
-1
). T
kl
подбирается так, чтобы обнулить элемент a
lk
матрицы T
kl
A
i
. Заметим, что при подобном преобразовании симметричность матрицы сохраняется. Поэтому условия обнуления элемента a
lk
совпадают с таковыми для элемента a
kl
. Последнюю задачу мы решали ранее, когда рассматривали QR разложение методом вращения.
Умножение матрицы T
kl
A
i на T
kl
T
справа может сделать элемент a
lk
снова ненулевым. Однако доказывается, что A
i
→D при k→∞.
Произведение матриц T
kl
справа от A
i формируют матрицу, столбцы которой есть собственные векторы. Перед очередным вращением проверяют выбранный элемент a
lk
на малость. Если выполнено условие
100| a
lk
|/|a kk
|
ε
маш.
, полагают элемент малым и вращение не производят.
Критерий окончания нет вращений в цикле, то есть все элементы в поддиагональной части малы.
14.4. Степенной метод.
Возьмем произвольный начальный вектор x
0
и построим последовательность векторов, используя формулы
( )
(
1)
(
1)
( )
1
(
1)
(
1)
(
,
)
(
,
)
k
k
k
k
k
k
k
x
Ax
x x
x
x







Справедлива следующая теорема.
Теорема. Пусть А – матрица простой структуры, для которой
выполнено условие │λ
1
│>│λ
2
│
≥
│λ
3
│
≥
…
≥
│λ
n
│. Предположим, что в
разложении x
0
= c
1
e
1
+c
2
e
2
+…+c
n
e
n
по базису из собственных векторов c
1

170
не равен нулю. Тогда λ
1
(k)
→
λ
1
при k
→∞
и верна оценка относительной
погрешности
( )
( )
1 1
2 1
1 1
|
|
(
)
|
|
k
k
k
C



 





Доказательство. Вектор x
(k)
= A
k x
0
. Так как A
k e
i
= λ
i k
e i
, то x
(k)
= A
k x
0
= λ
1
k
c
1
e
1
+λ
2
k
c
2
e
2
+…+λ
n k
c
n
e
n
.
Положим
( )
( )
1
k
k
k
x
z


и заметим, что z
(k)
→ c
1
e
1
при k
→ ∞, причем
( )
( )
1 1 2
1
k
n
k
k
i
i i
i
w
z
c e
c e



 



 
 

Так как
2 1
1
i





для всех i
≥ 2, то
( )
2 1
||
||
k
k
w
c



. Действительно
2 2
1 2
1 1
1 1
2 2
||
||
* |
| * ||
||
|
| * ||
||
|
| * ||
|| .
k
k
k
k
n
n
i
i
k
i
i
i
i
i
i
k
n
i
i
i
i
w
c
e
c
e
c
c
c
e




























Теперь нетрудно установить, что
(
1)
( )
(
1)
( )
1 1 1 1 1
1 1
1
(
1)
(
1)
(
1)
(
1)
1 1 1 1
(
,
)
(
,
)
(
,
)
при k
(
,
)
(
,
)
(
,
)
k
k
k
k
k
k
k
k
k
x x
z
z
c e c e
x
x
z
z
c e c e














 
Кроме того
( )
( )
(
1)
(
1)
( )
(
1)
(
1)
1 1
(
1)
(
1)
(
1)
(
1)
1
(
,
)
(
,
)
=
(
,
)
(
,
)
k
k
k
k
k
k
k
k
k
k
k
z
z
z
w
w
z
z
z
z
z















Используя неравенство Коши – Буняковского |(x,y)|
||x|| ||y||, имеем
( )
(
1)
(
1)
( )
(
1)
(
1)
| (
,
) | ||
||
|| .
k
k
k
k
k
k
w
w
z
w
w
z







Откуда следует
( )
( )
(
1)
( )
(
1)
1 1
2 0
(
1)
(
1)
1 1
|
|
||
||
||
|| ||
||
|
|
||
||
||
||
k
k
k
k
k
k
k
k
w
w
w
w
C
z
z















Действительно,
2 1
||
||
k
k
w
c



Тогда

171 1
*
2 2
2 2
2 1
1 1
1 1
1 1
*
(
1)
*
1 2
0
*
2
||
|| ||
||
1
,
||
||
;
k
k
k
k
k
k
k
w
w
c
c
c
c
z
c C
c



































З а м е ч а н и е 1. Теорема справедлива для произвольных матриц при указанном соотношении модулей собственных чисел.
З а м е ч а н и е 2. Если |λ
1
| > 1, то ||x
(k)
||
→
∞
при k
→
∞
и при
вычислениях на ЭВМ возможно переполнение. Если же |λ
1
| < 1
возможно исчезновение порядка. Для предупреждения этих ситуаций
обычно x
(k)
нормируют так, чтобы ||x
(k)
|| = 1 для всех k
≥
1. Возможна
следующая модификация: y
(k)
=Ax
(k-1)
, λ
(k)
= (y
(k)
, x
(k-1)
), x
(k)
= y
(k)
/||y
(k)
||.
Предполагается, что ||x
(0)
|| = 1. Последовательность y
(k)
сходится к
собственному вектору по направлению.
Важными достоинствами метода являются простота,
возможность эффективного использования разреженности матрицы
и отсутствия необходимости ее преобразования.
14.5. Метод обратных итераций для вычисления собственных
векторов.
Большинство методов вычисления собственных значений плохо приспособлены для вычисления собственных векторов. Собственные векторы являются ненулевыми решениями вырожденной однородной системы (A – λ
j
E)x = 0. Если вычислено приближенное значение λ
j
, матрица этой системы оказывается невырожденной и единственное ее решение равно тривиальному. В методе обратных итераций приближения к собственному вектору определяются последовательным решением систем уравнений (A – λ
j
*E)y
(k+1)
= x
(k)
с последующей нормировкой решения x
(k+1)
= y
(k+1)
/||y
(k+1)
||. В качестве начального приближения x
(0)
берут произвольный, часто единичный вектор.
Поясним механизм действия метода на примере матрицы простой структуры и в предположении, что λ
j простое собственное число. Векторы x
(0)
и y
(1)
представим в виде линейной комбинации собственных векторов

172
(0)
(1)
1 1
, y
n
n
i i
i i
i
i
x
c e
d e






Так как
*
(1)
*
1
(
)
(
)
n
j
i
i
j
i
A
E y
d

 





, то система уравнений метода принимает вид
*
1 1
(
)
n
n
i
i
j
i
i i
i
i
d
e
c e
 






Приравнивая коэффициенты при e
i
, получим
*
i
i
i
j
c
d
 


. Следовательно
*
(1)
*
*
*
1 1
(
).
n
n
j
j
i
i
j
j
i i
i
i j
i
j
i
j
i
j
c
y
e
c e
c e
 
 
 
 











Если |λ
j
– λ
j
*|
≪
|λ
i
– λ
j
*| для всех i
≠
j, то второе слагаемое в
правой части последней формулы мало по сравнению с первым.
Поэтому
(1)
*
1
j
j
i
j
y
c e
 


оказывается близким по направлению к вектору e
j
.
Можно показать, что вычисляемый на k – ой итерации вектор
*
( )
( )
*
(
)
n
j
j
k
k
j
j
i i
i j
i
j
x
c e
c e
 

 






,
где |β
(k)
|
→ |c j
-1
| при k
→ ∞. Вектор x
(k)
сходится к e
j
по направлению со скоростью геометрической прогрессии, знаменатель которой
*
*
|
|
max
|
|
j
j
i j
i
j
q
 
 




Если |λ
j
– λ
j
*|
≪
|λ
i
– λ
j
*| для всех i
≠
j, то метод обратных
итераций сходится очень быстро.
З а м е ч а н и е 1. Так как λ
j
* почти совпадает с λ
j
, то матрица
A – λ
j
*E плохо обусловлена и возникают опасения, что большие
погрешности при решении систем уравнений повлияют на точность
итерационных приближений. К счастью это не так. Возникающая
ошибка оказывается почти пропорциональной вектору e
j
. В данном
случае плохая обусловленность не ухудшает, а улучшает ситуацию.

173
З а м е ч а н и е 2. Записав формулу метода в виде
y
(k+1)
=(A – λ
j
*E)
-1
x
(k)
, замечаем, что метод обратных итераций это
степенной метод, примененный к матрице (A – λ
j
*E)
-1
14.6. Упражнения.
Объектом исследования является метод Якоби для действительных симметрических матриц. Такие матрицы ортогонально подобны диагональной матрице и имеют вещественные собственные значения.
Задача заключается в детализации метода до уровня алгоритма и его программная реализация на языке Matlab.
Как известно, итерационный цикл в методе Якоби состоит в переборе n(n-1)/2 элементов в поддиагональной части матрицы. После выделения очередного элемента a ik
(i > k) формируют матрицу вращения
T
ki
, обнуляющую этот элемент при умножении T
ki
А. Используя матрицу
T
ki проводят подобное преобразование T
ki
АT
ki
Т
. Очередное преобразование не производят, если элемент a ik мал. Критерием малости является выполнение условия 100|a ik
|/|a ii
|
ε
маш
. Критерием окончания процесса является отсутствие преобразований в цикле.
Алгоритм.
1. Ввод исходных данных: матрица А, заготовка матрицы собственных векторов q1=E, заготовка диагональной матрицы собственных значений D=A, вспомогательная переменная, хранящая число циклов w=0, вспомогательная переменная, подсчитывающая количество отсутствующих преобразований в цикле, s2=0. Критерием окончания процесса является s2
≥n(n-1)/2.
2. Пока одновременно выполнены условия s2
, положить s2=0, w=w+1,проводить итерационный процесс. Процесс заканчивается , если условия перестают выполняться.
В последующих пунктах описываются операции итерационного процесса.

174 3. Осуществить перебор элементов ниже главной диагонали. Для этого:
3.1 Выбрать очередной столбец k=1,..,n-1.
3.2 В очередном столбце выбрать очередную строку i=k+1, …,n и соответствующий элемент a ik
3.2.1. Проверить a ik на малость 100|a ik
|/|a ii
|
ε
маш
. При выполнении условия малости положить s2=s2+1 и вернутся к п. 3.2.
В противном случае формировать матрицу T
ki
, осуществлять подобное преобразование, формировать матрицу собственных векторов, положить s2=0.
2 2
2 2
, s
, t=E, T
,
,
, 1 1.
kk
ik
ki
kk
ik
kk
ik
T
kk
ii
ki
ik
a
a
c
t
a
a
a
a
t
t
c t
s t
s D
t D t q
t q





 

 


Вернуться к п. 3.2.
Ниже приведена соответствующая программа на языке Matlab. function
EiJacob3
a=[1 2 3 4 5;2 8 -7 -2 3 ;3 -7 2 1 5;4 -2 1 7 2;5 3 5 2 0];
%a=[6.25 -1 .5;-1 5 2.12;.5 2.12 3.6];
%a=[1 2 3 4;2 8 -7 -2;3 -7 2 1;4 -2 1 3];
n=length(a);q1=eye(n);s2=0;
w=0;D=a;
while
(s2 w=w+1;
s2=0;
for k=1:n-1
for i=k+1:n if
100*(abs(D(i,k)))/abs(D(i,i))else
%Расчет элементов матрицы Tki d=sqrt(D(k,k)^2+D(i,k)^2);
c=D(k,k)/d;s=D(i,k)/d;
%Формирование матрицы Tki ( t )
t=eye(n);
t(k,k)=c;t(k,i)=s;
t(i,i)=c;t(i,k)=-s;
%Расчет собственных чисел (диаг. матрица D)
% и собственных векторов (матрица q1')
D=t*D*t';q1=t*q1;
s2=0; end end

175 end end
D,q1',s2,w, [e,L]=eig(a)
Методы минимизации.
15.
Введение в минимизацию функций.
Требуется определить минимум вещественной функции n переменных f(x
(1)
,…,x
(n)
)=f(x), где x n-мерный вектор. Максимум f(x) совпадает с минимумом (-f(x)). Функция f(x
(1)
,…,x
(n)
) определена на множестве S n-мерного пространства. Если S совпадает со всем n-мерным пространством, задачу минимизации называют безусловной. В противном случае задача имеет ограничения, заданные системой равенств и неравенств. Функцию f(x) часто называют целевой функцией. Если f(x) и все ограничения линейны, говорят о задаче линейного программирования.
Если хотя бы одна из этих функций не линейна, имеем дело с задачей нелинейного программирования. Задача с ограничениями существенно сложнее безусловной задачи. Ограничимся рассмотрением задачи безусловной минимизации. Обозначим x* точку минимума функции f(x).
Если f(x*) R
n
, то x* - глобальный минимум. Будем рассматривать методы поиска локальных минимумов.
15.1.
Минимизация функции одной переменной.
15.1.
1. Постановка задачи. Определения.
Предположим, что для f(x), определенной на [a,b] имеется единственное значение x=x*, отвечающее минимуму f на[a,b]. При этом для x>x* функция строго возрастает и для x

176 можно выделить промежуток [a,b], внутри которого функция окажется унимодальной.
Если точка минимума x* находится внутри отрезка [a,b], необходимое условие минимума имеет вид f'(x*)=0. Существование минимума гарантировано, если одновременно выполнено условие f''(x*)>0.
Обусловленность задачи минимизации.
Пусть значения функции f вычисляются с погрешностью
|f(x)-f*(x)|
Δ. Вблизи точки минимума справедливо приближенное равенство f(x)
≈f(x*)+(1/2)f''(x*)(x-x*)
2
Оценим минимальное расстояние между точками x и x*, начиная с которого заведомо будет выполнено неравенство f*(x)>f*(x*). Это расстояние назовем интервалом неопределенности точки x* локального минимума и будем обозначать его ρ.
Имеем f*(x)-f*(x*)=f(x)-f(x*)+[(f*(x)-f(x))-(f*(x*)-f(x*))]
≈
≈(1/2)f''(x*)(x-x*)
2
+[(f*(x)-f(x))-(f*(x*)-f(x*))]
≥
≥ (1/2)f''(x*)(x-x*)
2
-2Δ.
Следовательно, неравенство f*(x)>f*(x*) выполнено, если
(x-x*)
2
≳4Δ/f''(x Таким образом, (x-x*)≈2 Δ/ x
Если задача хорошо масштабирована, т е x*

1, |f(x*)|1, f''(x)1, то последнее соотношение можно записать в терминах относительных погрешностей δ(x*) δ Отсюда следует, что если δ(f*)10
-m
, то δ(x*)10
(-m/2)
потеря половины верных значащих цифр в результате
Полученный результат основан на анализе свойств функции f(x). Если для поиска минимума решается нелинейное уравнение f'(x)=0, то радиус неопределенности вычисляется по ранее найденной формуле ρ Δ/|f
(2)
(x*)|.
15.1.
2. Метод деления отрезка пополам.

177
В этом методе используется принцип последовательного сокращения отрезка локализации, основанный на следующем соображении.
Пусть f унимодальная на отрезке [a,b] функция и a
α<γ<β b. Тогда
1
0
)если f(α)
f(β), то x* [a,β],
2
0
) если f(α)
≥f(β), то x* [α,b],
3
0
) если f(α)
≥f(γ) и f(γ) f(β), то x* [α,β].
1
0
. Предположим противное x*>β. Тогда вследствие
унимодальности f получим f(α)
>f(β), что противоречит условию.
2
0
.Предположим противное x*<α. Тогда вследствие унимодальности
f получим f(α)
что противоречит условию.
3
0
. В силу п.1
0
имеем x*
[a,β], а всилу п.2 0
x*
[α,b]. Следовательно,
x*
[α,β].
Дополнительно заметим, если функция унимодальна на отрезке [a,b], то она унимодальна на любом отрезке [c,d]
⊂[a,b].
Обозначим отрезок [a,b] через [a
0
,b
0
]. Выберем на отрезке [a
0
,b
0
] две симметрично расположенные точки α
0
=[(a
0
+b
0
)/2]-δ, β
0
= [(a
0
+b
0
)/2]+δ.
Здесь δ параметр метода, δ>0. Далее вычисляют f(α
0
) и f(β
0
). Сравнение этих значений позволяет сократить отрезок локализации следующим образом если f(α
0
) f(β
0
), то x
[a
1
,b
1
]=[a
0
,β
0
];
0 0
0 1
(
);
2 2
b
a




 
 

если f(α
0
)
>f(β
0
), то x
[a
1
,b
1
]=[α
0
,b
0
]. Как и в варианте1 0
1
;
2


 

a0
b0 2
δ
Вариант
1; ((b0-a0)/2+
δ)
Вариант
2;
((b0-a0)/2+
δ)
α0
β0

178
За очередное приближение к точке минимума в первом случае принимают x
1
*
=α
0
, во втором x
1
*
=β
0
. Описанная процедура представляет один шаг итерационной последовательности.
Обозначим через Δ
n
=b n
-a n
длину отрезка [a n
,b n
]. Очевидно
Δ
n+1
= [Δ
n
/2]+δ=[(Δ
n
-2δ)/2]+2δ. Значения Δ
n образуют цепочку
Δ
1
=[(Δ
0
-2δ)/2]+2δ, Δ
2
=[(Δ
1
-2δ)/2]+2δ=[(Δ
0
-2δ)/2 2
]+2δ,…, Δ
n
=[(Δ
0
-2δ)/2
n
]+2δ.
Метод половинного деления сходится о скоростью геометрической прогрессии при q=1/2.
Δ
n
→2δ при n→∞. Если задана погрешность ε определения минимума, следует положить |x n
-x*|-a n
=Δ
n
< ε .Причем должно быть ε>2δ.
Рекомендуется задавать δ≤0.2ε.
15.1.
3. Метод золотого сечения.
Метод золотого сечения позволяет обойтись на каждом шаге итераций вычислением лишь одного значения функции вместо двух.
Золотым сечением отрезка называют такое его разбиение на две неравные части, что отношение длины всего отрезка к длине его большей части равно отношению длины большей части к длине меньшей части.
Золотое сечение отрезка [a,b] осуществляется каждой из двух симметрично расположенных относительно центра отрезка точек
2
(
)
0.382(
) ,
3 5
2
=
(
)
0.618(
) .
1 5
a
b
a
b
a
a
b
a
b
a




 




















179 a
b
α
β
0.618 0.382
Вариант
1
Вариант
2
Действительно
1 5
1 5
,
;
2 2
b
a
b
b
a
a
b
a
a
b


 
















Замечательно то, что точка α осуществляет золотое сечение не только отрезка [a,b], но и отрезка [a,β]
1 5
2
a
a
a



 







Точно также точка β осуществляет золотое сечение не только отрезка
[a,b], но и отрезка [α,b]. Этот факт далее используется.
Очередная (k+1)-я итерация осуществляется аналогично (k+1)-й итерации метода половинного деления. Новый отрезок локализации определяется по тому же правилу, как и в методе половинного деления.
Если f(α
k
)
f(β
k
), x
[a k+1
,b k+1
] =[a k
,β
k
] .
Если f(α
k
)
>f(β
k
), x
[a k+1
,b k+1
] =[α
k
,b k
].
Важно то, что какой бы из отрезков [a k
,β
k
] или [α
k
,b k
] не был выбран за очередной отрезок локализации, точка очередного приближения x k+1
совпадает с одной из точек α
k или β
k
(в первом случае x k+1
=α
k
, а во втором случае x k+1
=β
k
). Причем, указанные точки уже осуществляют золотое сечение новых отрезков. Требует определения лишь одна дополнительная симметричная точка золотого сечения для нового отрезка. Поэтому

180 достаточно вычислять значение функции лишь в этой одной недостающей точке.
Точка x k+1
отстоит от конца отрезка [a k+1
,b k+1
] на величину, не превышающую
2 2
(
)
1 5
1 5
k
k
k
b
a





. Таким образом,
|x k+1
-x*|
1 2
1 5
k
k

  

Каждая итерация сокращает длину отрезка локализации в
5 1 2

раз Поэтому b k
-a k
=Δ
k
=
0 2
5 1
k

 





Таким образом, метод сходится со скоростью геометрической прогрессии, знаменатель которой
2 0.62.
5 1
q



Критерий прекращения итерационного процесса Δ
k
<ε, где
ε – заданная погрешность