Практикум j практическое примщенше численных методов

Ai = min

А„ = max

(Ах, х)


тих / v ^
Х?0 (х, X)


(Ах,х)


.ИСЬЛ j v •
х^о (Х,Х)

ажны также экстремальные свойства Релея
Для локализации собственных значений произвольной матрицы А привл' ются круги Гершгорина: любое собственное значение матрицы А лежи] крайней мере в одном из кругов
п
| А (1ц | ^ ^ ^ \o>ij |) i ⁼ 1,2 ... ,71.
*Фз=1
Приведем теперь некоторые о возмущении собственных значений при во: щении элементов матрицы. Помимо (6.1) рассмотрим задачу
Ар> = А (р. (
Ограничимся случаем, когда все собственные значения простые. С то¹ сгью до членов второго порядка возмущение собственных значений за < возмущения матрицы дается оценкой
|А - А| < Ci\\A - А\\, (
где ||х|| = (.т,ж)^1/2. Мерой усггойчивости собственного значения А* слу: величина
которая называется коэффициентом перекоса матрицы А, соответствующим данному собственному значению. Здесь ^ — собственное значение матрицы А*. Для нормированных собственных векторов задач (6.1) и (6.5) соответ­ствующая оценка устойчивости имеет вид
ij£j=l I * I
В частности, для симметричной матрицы все коэффициенты перекоса рав­ны единице и оценки устойчивости вычисления собственных значений опти­мальны.
Итерационные методы решения частичной проблемы
собственных значений
Для нахождения минимального по модулю (максимального) собственного значения используется прямые и обратные итерации. Пусть матрица А яв­ляется симметричной, все ее собственные значения простые и упорядочены следующим образом
N < |А₂| < ••• < |А„|.
Определим последовательность векторов
x
(6.7)

^k+l = Ах^к, /г = 0,1,...
п
(х^к+1,х^к)
(х^к₎х^к)


Ап + О



(6.8)

ри некотором заданном х° (прямые итерации). Рассматривая последователь­ности скалярных произведений (х^к,х^к), (х^к+1,х^к), при ограничении (■Х°,(р^п)Ф О, получим
Тем самым при использовании итерационного процесса (6.7) находится мак­симальное по модулю собственное значение матрицы А.
Принимая во внимание, что собственные значения матрицы А~¹ есть 1/А*, i = 1,2,... ,п, при использовании последовательности (обратные итерации)
2/^fc+1=A-y, (2/VyO, к = 0,1,... (6.9)
имеем
Шп ^{уШ’^ук) = 1.
(_?/У) АГ
Тем самым при обратных итерациях находится минимальное по модулю соб­ственное значение матрицы.
Заметим, что в прямых и обратных итерациях нет необходимости в специ­альном вычислении соответствующих собственных векторов, так как
l
(6.10)

iin х^к = (р^п, lirn у^к = (р¹.
Вычислительная реализация обратных итераций (6.9) может быть оспог на однократном LU разложении матрицы А. После этого каждая обрат итерация по вычислительным затратам эквивалентна прямой итерации.
Отметим процедуру ускорения сходимости обратных итераций при извест хорошем приближении к собственному значению Ах к собственному значе] A_x. В этом случае рассматриваются обратные итерации со сдвигом, когд;
z^k+1 = (А - \_xE)~¹z^k, (z°, у?¹) Ф 0, & = 0,1, —
Скорость сходимости обратных итераций со сдвигом и без сдвига опред( ется отношениями

Ах — Ах		Ai
ю 1	5	а2

соответственно. В более общей ситуации обратные итерации со сдвигом пользуются для нахождения ближайшего к заданному числу собствен* значения соответствующего собственного вектора.
Решение полной проблемы собственных значений
Прямые и обратные итерации хорошо приспособлены для определения дельных собственных значений и собственных векторов. Для решения а тральной задачи в целом используется QR алгоритм. Он основан на предел лении матрицы А в произведение А = QR, где Q — ортогональная Q*Q = a R — верхняя треугольная матрицы.
Строится последовательность ортогональных матриц Qk и верхних треугс ных матриц Rk по рекуррентным формулам
А = Q\Ri₎ А\ — RiQu А_х = Q2R2, А₂ = R2Q2,
Ак-1 = QkRk) А_к = RkQk,
Процесс построения по матрице А матриц Qk, Rk, Ак называется QR-aj ритмом.
Пусть для невырожденной матрицы А собственные значения различны модулю и
|Ai| > |Ах| > • • • > |А_П|
и существует представление
А = Т\Т~\ T~^l = LU, A = diag{A₁,A₂,...,A_n}.
Тогда последовательность матриц А_к <ЗЯ-алгоритма сходится к верхней тре­угольной матрице, а диагональные элементы — к собственным значениям матрицы А.
При решении полной спектральной задачи на основе QR для минимизации вычислительной работы проводится предварительное преобразовании исход­ной матрицы к верхней почти треугольной матрице, в которой а^ = 0, г > j + 1. При рекуррентных преобразованиях Q 72-алгоритма матрицы А_к оста­ются почти верхнетреугольными.
Решение спектральной проблемы для симметричной вещественной матрицы может осуществляться методом вращений (методом Якоби). Вещественная матрица, отличающаяся от единичной матрицы четырьмя элементами, распо­ложенными на пересечении строк и столбцов с номерами z,j, и имеющая вид
T(ij) =
-s • • • С
где с² + s² = 1 называется матрицей вращения. Заметим, что матрица враще­ния является ортогональной и при умножении вектора на матрицу вращения T(ij) меняются только г и j координаты вектора.
Для любой матрицы А и любой пары индексов г, j, г ф j всегда можно найти такую матрицу вращения T(ij), что элемент Ъц матрицы В = T*(ij)AT(ij) равен нулю. Определим последовательность матриц
Aq = Л, А\, А_2}...,
каждая из которых получается из предыдущей с помощью преобразования подобия, определяемой матрицей вращения. На каждом шаге этого процесса обнуляется отдельный внедиагональный элемент. При таком преобразовании сумма квадратов внедиагональных элементов убывает. Последовательность матриц А_к сходится к диагональной матрице и диагональные элементы мат­рицы А_к в рассматриваемом методе вращений (метод Якоби) являются соот­ветствующими приближениями для собственных значений матрицы А.
Оптимизация метода вращений достигается за счет выбора элемента для уц_и^ чтожения на каждом шаге преобразований. Это может быть максимальны^ по модулю виедиагональный элемент всей матрицы А_к или на выбранное столбце.

3. 
   Упражнения
   

Упражнение 6.1 Напишите программу для пахооюдения минимального по модулю собственного значения и соответствующего собственного вектора симметричной матрицы при использовании обратных итераций. С ее по- мощью найдите минимальное по модулю собственное значение матрицы Гильберта А, когда
dij = т—т—г = 1,2,..., гг, j = l,2,...,n, i + j- l
при значениях п от 2 до 10.
В модуле invlter функция invIterO обеспечивает нахождение минимального по модулю собственного значения и соответствующего собственного вектора с помощью обратных итераций (см. (6.9)). Для решения системы линейных уравнений используется функция solveLU из модуля lu.
import numpy as np
import math as mt
from lu import decLU, solveLU
def invIter(A, tol = 1.0e-6):
и и и
The inverse power method.
Returns the smallest eigenvalue and the corresponding eigenvector.
и и и
n = len(A)
iterMax = 100
x = np.random.rand(n)
xNorn = mt.sqrt(np.dot(x, x))
x = x / xNorn
for i in range(0, iterMax): xl = x.copyO x = solveLU(A, x) xNorn = mt.sqrt(np.dot(x, x)) x = x / xNorn
if mt.sqrt(np.dot(xl-x, xl-x)) < tol: return 1 / xNorn, x
print ’Method did not converge (100 iterations).’
Решение частичной проблемы собственных значений для матрицы Гильберта дается следующей программой.
import numpy as np from invlter import invlter print ’n eigenvalue’ for n in range(2, 11):
A = np.zeros((n, n), ’float’) for i in range(0, n):
for j in range(0, n):
A[i,j] = 1./(i+j+1) lam, x = invIter(A) print n, lam
Wllige n v a 1 u 0 /.; SI
1||Й7б2ШШ&Ш
Н^^ИШНШЙВ