Численные методы. Соснин (2005) (1160462), страница 7

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 7)

ДВУХШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ ВАРИАЦИОННОГО ТИПА35Метод сопряженных градиентовВ данном методе берется D = A, поэтому матрица A должна быть симметричной и положительноопределенной. Как уже показывалось ранее, в этом случае для τk+1 можно получить такую формулу: k kω ,rτk+1 =.hAω k , ω k iВ свою очередь, для αk+1 будет справедливо такое выражение: k k !−1ω ,rτk+1αk+1 = 1 −· k−1 k−1τk αk hω,ri— в состав этих формул входит уже не погрешность, а невязка и поправка, которые мы вычислятьумеем.Метод сопряженных невязокВ данном случае берется D = AT A, и несложно проверить, что D = DT > 0. Получающиеся формулыдля τk+1 и αk+1 мы опускаем (вывод предоставляется читателю).Кроме того, в данном методе и во всех последующих вводится дополнительное требование:DB −1 A = (DB −1 A)T > 0.(1.42)В данном случае оно эквивалентно тому, что B T A > 0.

Это небольшое ограничение на матрицу B.Метод сопряженных поправокЗдесь мы берем D = AT B −1 A. Легко проверить, что D = DT > 0, если B = B T > 0. Ограничение(1.42) приводит к требованию положительной определенности матрицы A. Формулы для τk+1 и αk+1также становятся корректными.Метод сопряженных погрешностейЗдесь все тоже просто: D = B0 , где B0 = B0T > 0, причем матрица B0 должна быть легкообратимой.Матрица B задается строго:B = (AT )−1 B0 .Ограничение (1.42) приводит к требованию невырожденности матрицы A.Общие замечания(1) Можно показать, что любой из выше перечисленных методов сходится не медленнее, чем соответствующий ему одношаговый итерационный метод с чебышевским набором параметров.(2) Кроме того, если количество шагов итерации превысит размерность системы, то последовательность итерационных приближений выйдет на точное решение — то есть вышеперечисленные методыфактически являются прямыми.

Другое дело, что нужная точность будет достигнута намного раньше.(3) Учитывая, что α1 = 1, можно не угадывать x1 (первое приближение), а просто подставить в(1.41) k = 0 :x1 − x0B+ Ax0 = fα1 τ1— это простая одношаговая схема для поиска x1 .Глава 2Задачи на собственные значенияСначала о терминологии. Пусть A = (aij ) — матрица размера n × n, а x = (x1 , x2 , . . . , xn )T — векторнеизвестных. Тогда поиск таких констант λ и векторов x, чтоAx = λx,называется задачей на собственные значения.

Нетрудно заметить, что эта задача эквивалентнапоиску таких x и λ, для которых(A − λE)x = 0.При этом λ называются собственными значениями, а соответствующие им вектора x — собственными векторами.Из курса линейной алгебры известно, что если det(A − λE) = 0, то решение существует. Для егонахождения надо найти такие λ, чтобы определитель |A − λE| обращался в ноль.

Этот определительявляется полиномом от λ с коэффициентами из A — поэтому, если n мало, то корни этого полиномалегко найти. Это также просто сделать, если матрица A является диагональной или верхнетреугольной: в этом случае определитель будет равен произведению диагональных элементов.Далее мы будем рассматривать различные методы нахождения собственных значений.Определение.

Метод решает полную проблему собственных значений, если он находит ихвсе. В противном случае (например, надо найти только границы спектра матрицы A) говорят, чторешается частичная проблема.2.1Поиск собственных значений методом вращенийДанный метод, предложенный К. Якоби, позволяет найти все собственные значения (решает полнуюпроблему поиска). Для этого требуется симметричность матрицы A : A = AT . Известно, что в этомслучае для A справедливо такое представление (оно единственно):A = QT Λ Q,(2.1)где Q — некая ортогональная матрица (QT = Q−1 ), а Λ — диагональная матрица, причем на диагонали у нее стоят именно собственные значения матрицы A : Λ = diag (λ1 , .

. . , λn ) .Если домножить равенство (2.1) на (QT )−1 = Q слева, а на Q−1 = QT — справа, то получим такуюформулу:Λ = QAQT(2.2)— то есть для нахождения собственных значений матрицы A нам необходимо найти соответствующуюматрицу Q и провести два матричных умножения.362.1. ПОИСК СОБСТВЕННЫХ ЗНАЧЕНИЙ МЕТОДОМ ВРАЩЕНИЙ37Матрицу Q будем находить с помощью ортогональных преобразований матрицы A, постепенноуменьшая абсолютные значения ее внедиагональных элементов:A1 = Vij1 A (Vij1 )T ;A2 = Vij2 A1 (Vij2 )T ;···Ak = Vijk Ak−1 (Vijk )T ,и так далее. Здесь Vijk — некие ортогональные матрицы, а индексы в них говорят о номере преобразования — (k ) и об индексе уменьшаемого элемента из Ak−1 — ( ij ). Произведение ортогональныхматриц дает ортогональную матрицу, поэтому, если мы на некотором шаге придем к диагональнойматрице, то это будет означать получение преобразования (2.2).Матрицы Vijk будут задаваться так (ϕ — пока свободный параметр):10. ..1cos ϕ− sin ϕ10k.Vij = ..000 1sin ϕcos ϕ1..

.01— здесь на диагонали стоят единицы везде, кроме vii и vjj , где стоят косинусы, а в позициях vji и vijстоят синусы. Все остальные элементы — нулевые. Легко проверить, что все эти матрицы ортогональны.Индексы i и j задаются на каждом шаге заново. Они обозначают индекс максимального по модулювнедиагонального элемента, то есть|aij | = max |aklm |,l,ml6=mгде aklm — элементы матрицы Ak (k-е приближение к Λ ). Если i и j можно выбрать несколькими способами, то берется произвольная пара. Если же все внедиагональные элементы — нулевые, то процесспрекращается.Назовем количественной характеристикой диагональности матрицы Ak такое число:t(Ak ) =n XnX(akij )2i=1 j=1j6=ik→∞Понятно, что если t(Ak ) −→ 0, то Ak сходится к диагональному виду.

Таким образом, через t(Ak )можно наглядно оценивать скорость сходимости данного метода.Теперь фиксируем максимальный элемент aij и установим соотношение между матрицами Ak+1 иAk . Ak+1 мы задавали так:Ak+1 = Vijk+1 Ak (Vijk+1 )T .Обозначим Bk = Ak (Vijk+1 )T = (bkms ) и (Vijk+1 )T = (v k+1lm ) — обозначение для элементов матрицыk+1 T(Vij ) . Тогда, согласно определению матричного произведения, ka ,s 6= i, j;n msXakmi cos ϕ − akmj sin ϕ, s = i;bkms =akmp v k+1(2.3)ps =p=1 kkami sin ϕ + amj cos ϕ, s = j.38Глава 2. ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ— эта сложная формула вытекает из того, что среди элементов v k+1один ненулевой при p 6= i, j, иpsдва ненулевых при p = i, j.Теперь, согласно принятым обозначениям, Ak+1 = Vijk+1 Bk . Выведем отсюда формулу для элементов матрицы Ak+1 : kb ,m 6= i, j;n msXk+1 kbkis cos ϕ − bkjs sin ϕ, m = i;ak+1vmpbps =(2.4)ms =p=1 kbis sin ϕ + bkjs cos ϕ, m = j.k+1— идея построения системы та же (vmp— элементы матрицы Vijk+1 ).Напомним, что aij — максимальный по модулю элемент из A.

Положив в (2.4) m = i и s = j,получим такую формулу для него:ak+1= bkij cos ϕ − bkjj sin ϕ.ijТеперь распишем bkij и bkjj через (2.3), взяв m = i, j и s = j :ak+1ij=(akii sin ϕ + akij cos ϕ) cos ϕ − (akji sin ϕ + akjj cos ϕ) sin ϕ ==(akii − akjj ) sin ϕ cos ϕ + akij (cos2 ϕ − sin2 ϕ) ={A = AT ⇒ Ak = ATk} =(akii − akjj ) sin 2ϕ+ akij cos 2ϕ.2Напомним, мы пытаемся уменьшить внедиагональные элементы матрицы A при ортогональномпреобразовании.

Потребуем равенство нулю для ak+1ij . Таким образом, написанное выше преобразование превращается в уравнение относительно ϕ. Решая его, получим:ϕ=2akij1arctg k.2aii − akjjВычислим количественную характеристику диагональности получившейся матрицы Ak+1 :t(Ak+1 ) =n XnX2(ak+1ms ) .m=1 s=1s6=mСогласно формулам, элементы в Ak при умножении на (Vijk+1 )T изменяются только в i-м и j-мстолбцах. Аналогично, элементы в Ak+1 изменяются относительно B k только в i-й и j-й строках. Тоесть Ak+1 отлична от Ak только в i-х и j -х строках и столбцах.Выделим в сумме неменяющиеся элементы. Распишем сумму, раскроем скобки, и проведем всепреобразования:t(Ak+1 )==nXnXm=1s=1m6=i,j s6=i,j,mnnXXm=1s=1m6=i,j s6=i,j,m(akms )2 +(akms )2 +nnXX k 2 k+1 22(bmi ) + (bkmj )2 +(ais ) + (ak+1=js )m=1m6=i,jnXs=1s6=i,j k 2(ami ) cos2 ϕ + (akmj )2 sin2 ϕ − 2akmi akmj sin ϕ cos ϕ+m=1m6=i,jnX k 2(bis ) cos2 ϕ + (bkjs )2 sin2 ϕ−+(akmi )2 sin2 ϕ + (akmj )2 cos2 ϕ + 2akmi akmj sin ϕ cos ϕ +s=1s6=i,j−2bkis bkjs sin ϕ cos ϕ + (bkis )2 sin2 ϕ + (bkjs )2 cos2 ϕ + 2bkis bkjs sin ϕ cos ϕ =2.2.

СТЕПЕННОЙ МЕТОД ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ=nXnXnnXX k 2 k 2(ami ) + (akmj )2 +(ais ) + (akjs )2 + 2(akij )2 − 2(akij )2 =(akms )2 +m=1s=1m6=i,j s6=i,j,m39m=1m6=i,js=1s6=i,j= t(Ak ) − 2(akij )2 .Из всей этой последовательности формул следует, что t(Ak+1 ) < t(Ak ) — как видно, характеристикадиагональности монотонно убывает с ростом индекса, причем уменьшается каждый раз на 2(akij )2 —удвоенный квадрат максимального элемента. Из этого следует, что описанный итерационный процессприводит к диагональной матрице.Запишем оценку на скорость сходимости процесса. Пусть akij — максимальный внедиагональныйэлемент.

Простым подсчетом элементов матрицы можно получить такое неравенство:t(Ak ) 6 n(n − 1)(akij )2 .t(Ak )для n > 2. Подставим это неравенство в ранее полученноеn(n − 1)соотношение между t(Ak ) и t(Ak+1 ) :Отсюда следует, что (akij )2 >t(Ak+1 ) = t(Ak ) − 2(akij )2 6 t(Ak ) −2t(Ak ) = qt(Ak ),n(n − 1)где q = 1 −2< 1.n(n − 1)Применив эту операцию k раз, получим:t(Ak ) 6 q k t(A).Из последнего неравенства видно, что процесс нахождения матрицы Q сходится к диагональнойматрице Λ со скоростью геометрической прогрессии со знаменателем q.Можно немного оптимизировать способ выбора «плохого» элемента: сначала выбираем «плохую»nX(asi )2 имеет максимальное значение), а потом из этойстроку (например, выбираем строку s, гдеi=1строчки выбираем максимальный по модулю элемент.

При этом, немного неоптимальный выбор askкомпенсируется скоростью нахождения «плохого» элемента (2n вместо n2 сравнений).Теперь будем решать частичные проблемы. Как пример частичной проблемы, можно привестизадачу нахождения границ спектра в итерационном методе.2.2Степенной метод поиска собственных значенийБудем рассматривать задачу нахождения максимального по модулю собственного значения симметрической матрицы A.Алгоритм поискаВозьмем любой вектор x0 , отличный от нуля, и построим последовательность векторов xk такую,чтоxk+1 = Axk , k = 0, 1, . .

Характеристики

Список файлов лекций