Численные методы. Соснин (2005) (1160462), страница 5
Текст из файла (страница 5)
ЧЕБЫШЕВСКИЙ НАБОР ИТЕРАЦИОННЫХ ПАРАМЕТРОВ25Оценим скорость сходимости этого метода при q 6= 0. По определениюqrrδδ1+34∆1δδ∆qqln = ln= ln 1 + 4+≈4.q∆∆1− δ1− δ∆∆Применив попеременно-треугольный итерационный метод к решению системы уравнений (1.29) стрехдиагональной матрицей A, получим оценку на k0 (ε) :ln 1k0 (ε) = qε .δ4 ∆δλmin (B −1 A)Как уже показывалось для метода простой итерации,== tg2∆λmax (B −1 A)откуда получаем оценку на k0 (ε) :k0 (ε) ≈ln 1ε4 πh2=πh2≈π 2 h2,4ln 1εN = {ε = 5 · 10−5 ≈ e−10 } ≈ 1, 6N.2πЭто уже неплохой результат (гораздо меньше, чем в случае метода простой итерации), например,взяв N = 100, мы должны будем сделать 160 итераций вместо 20000 в МПИ.1.6Чебышевский набор итерационных параметровВ данном разделе мы будем применять для уже порядком заколебавшего нас уравнения Ax = f такуюитерационную схему:xl − xl−1+ Axl−1 = f.(1.33)BτlТеперь фиксируем число итераций (скажем, k, а то m — маловато...), и постараемся выбрать τlтак, чтобы погрешность на k-й итерации была минимальной:||z k || −→ inf .τlЧтобы получить оценку на погрешность, перейдем от l-го приближения к погрешности z l = xl − x,как это уже делали раньше сто раз.
Тогда итерационный процесс будет выглядеть просто прекрасно:Bz l − z l−1+ Az l−1 = 0.τlОтсюда без проблем получается выражение для z l :z l = (E − τl B −1 A)z l−1 ,l = 1, k.Рекурсивно применяя эту формулу для z k , получим очень-очень длинную формулу:z k = (E − τk B −1 A)(E − τk−1 B −1 A) . . . (E − τ1 B −1 A)z 0 .(1.34)Для тех, кто в танке, повторяем: необходимо так подобрать τl , чтобы минимизировать z k . Оказывается, что эта задача не только имеет решение, но нам даже разрешили его опубликовать (тожемне, ценность великая...). Это решение, а, точнее, набор τl называется чебышевским наборомитерационных параметров.
Прямо корпускулярно-волновая теория света, мдаа...А теперь приколитесь — следующая теорема пойдет без доказательства (кому неймется, можетпосмотреть в [1]).26Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙТеорема 1.6. Пусть матрицы A и B симметричны и положительно определены, τl считаютсяпо следующей формуле (k — фиксировано):2; τ0 = λ−1A) + λmax (B −1 A)min (Bτ01−ξλmin (B −1 A)τl =, гдеρ0 =, ξ=;1 + ρ0 tl1+ξλmax (B −1 A) tl = cos (2l − 1)π , l = 1, k,2kа xl вычисляется по формуле (1.33).Тогда погрешность ||xk − x||A будет минимальной среди всех возможных, и для нее справедливонеравенство:||xk − x||A 6 qk ||x0 − x||A ,√ρk11− ξ√ .где qk =,ρ=11 + ρ2k1+ ξ1Замечание 1.
Подробную инструкцию о том, что надо выкурить, чтобы догадаться до такихформул, можно найти в [2].Замечание 2. Найдем, при каких k выполняется условие остановки итерационного процесса:||xk − x||A 6 ε||x0 − x||A .(1.35)Из теоремы 1.6 следует, что оно заведомо (т.е. всегда (т.е. даже ночью)) выполняется, если qk 6 ε.Это, в свою очередь, эквивалентно тому, чтоρk16 ε ⇐⇒ ε(ρk1 )2 − ρk1 + ε > 0.1 + ρ2k1Корнями этого квадратного трехчлена будутρk1 =1±√ kρ1 = ε;1 ± (1 − 2ε2 )1 − 4ε2≈=⇒ 12ε2ερk1 = + ε.εОтсюда следует, что qk 6 ε приρk1 6 ε;ρk1 >1+ ε.εВ теореме 1.6 мы выбирали ρk1 так, что оно было меньше единицы. Поэтому второй случай нам неподойдет ( 1ε 1), из чего следует, чтоqk 6 ε ⇐⇒ ρk1 6 ε.Таким образом, (1.35) верно при&k > k0 (ε) =ln 1εln ρ11'.В данном случае за скорость сходимости принимается ln ρ11 .
Оценим этот параметр:√pp11+ ξ√ = ln 1 + 2 ξ + O(ξ) ≈ 2 ξ.ln= lnρ11− ξ(1.36)1.6. ЧЕБЫШЕВСКИЙ НАБОР ИТЕРАЦИОННЫХ ПАРАМЕТРОВ27Посмотрим, как работает этот метод (1.33) (для упрощения сделаем матрицу B единичной) применительно к нашей модельной задаче (1.28). Напомним, матрица системы в примере была вида:2 −1−1 2 .
. . 01 ......A= 2....h......0−1−1 2Собственные значения для нее выглядят так:λmin =4πhsin2 ( ),h22λmax =4πhcos2 ( ),h22minпоэтому ξ = λλmax= tg2 ( πh2 ).И снова, в который раз (и не в последний, по-моему) зададимся ε = 0, 5 · 10−4 ≈ e−10 . В этом случаеиз (1.36) можно поиметь следующую формулу:k0 (ε) =1011010= {h = } =N ≈ 3, 5N.≈πhπhNπ2 tg( 2 )Тогда для N = 10 нам понадобится 35 итераций, а для N = 100 — 350 итераций, то есть ужеможно себя не убивать...
Откровенно говоря, этот метод немеряно круче (на порядок) метода простойитерации, да и то, что схема явная, не может не радовать.С другой стороны, этот метод медленнее, чем попеременно-треугольный, но там-то надо было матрицы обращать, а здесь этого иногда можно и не делать...Упорядоченный набор чебышевских параметровВообще говоря, в описанном в предыдущем разделе процессе не предполагалось монотонности убывания погрешности.
Обычно ее и нет, что на практике при реализации метода с чебышевским наборомпараметров часто приводит к возникновению неприятных ситуаций (переполнение регистров и прочийотстой). Для предотвращения таких дел набор параметров можно особым образом упорядочить, чтобыкак-то скомпенсировать увеличение погрешности за счет ее уменьшения на соседних итерациях.К примеру, в нашем методе отличия на разных итерациях заключались в изменяющемся параметреτl , который в свою очередь зависит от tl = cos( (2l−1)π). По-другому это можно записать так:2kπ(2l − 1)πtl = cos= cosθlk , θlk = {1, 3, . .
. , 2k − 1}.2k2kОдним из способов упорядочивания параметров τl является простая перестановка значений θlk . Мырассмотрим пример такой перестановки, когда k является степенью двойки: k = 2p . Тогда значениябудут считаться по рекуррентным формулам:θ11 = 1;2mθ= θim ,i = 1, m;(1.37) 2i−12m2mθ2i= 4m − θ2i−1, i = 1, m.— последние две формулы необходимо применять для всех m = 1, 2, 4, . . . , 2p−1 .
Проще это будетпонять на примере:k = 24 = 16 :28Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙθ112θ1,24θ1,...,48θ1,...,8m=1:m=2:m=4:1;1, 4-1=3;1, 8-1=7, 3, 8-3=5;1, 15, 7, 9, 3, 13, 5,11;— потом по θlk строятся tl и τl , при этом получается более устойчивый метод.В следующей теореме, которую мы опять примем без доказательства, будет показано, как можно соединить попеременно-треугольный итерационный метод с упорядоченным набором чебышевскихпараметров.Теорема 1.7. Пусть уравнение Ax = f решается по следующей итерационной схеме (число шагов —k — фиксировано):xl − xl−1B+ Axl−1 = f, l = 1, k,τlгде B = (E + ωR1 )(E + ωR2 ), R1 + R2 = A (матрицы R1 , R2 имеют тот же смысл, что и в теореме1.4).Пусть также существуют такие положительные константы δ и ∆, чтоA > δE,Тогда, если выбиратьτ0 =τl = tl =∆A > R 1 R2 .4итерационные параметры τl так:2;γ1 + γ2√δ ∆√ ;√=γ1τ01−ξ2( ∆ + δ), ξ= ,, ρ0 =√1 + ρ0 tl1+ξγ2δ∆ γ2 =,4πcosθk , θlk считаются по формулам вида (1.37).2k l γ1то для достижения точности||A(xk − x)||B −1 6 ε||A(x0 − x)||B −1ln 2εитераций.достаточно выполнить k0 (ε) = √ qδ2 2 4 ∆Доказательство.
Доказательство можно найти в [1].Теперь посмотрим, какие улучшения даст этот метод в модельной задаче (1.28). Напомним, чтоконстанты δ, ∆ можно взять такими:42 πh δ = λmin (A) = 2 sin ( );h24πh ∆ = λmax (A) =cos2 ( ).h22Тогда при ε = 0, 5 · 10−4 ≈ e−10 получим, что√√ln 2 + 10ln 2 + 101(ln 2 + 10) N√k0 (ε) = √ q≈ √ q= {h = } =≈ 3 N.N2 π2 2 4 tg2 ( πh2 2 πh2 )21.7. ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ ВАРИАЦИОННОГО ТИПА29Таким образом, для N = 10 нам понадобится 10 итераций, а для N = 100 — 30 итераций. Превосходство в скорости над ранее рассмотренными методами очевидно, однако недостатки тоже налицо:необходимо уметь обращать матрицу B и много знать о спектре матрицы A.
Заметим, что это общиетребования у быстрых итерационных методов — но, как мы увидим далее, они не являются необходимыми.1.7Одношаговые итерационные методы вариационного типаВ данном разделе мы вновь будем работать с одношаговыми методами. Запишем итерационную схемуметода:xk − xk−1B+ Axk−1 = f.τkКак и раньше, перейдем от xk к погрешности z k = xk − x :Bz k − z k−1+ Az k−1 = 0.τkОтсюда выразим (k + 1)-ю погрешность:z k+1 = (E − τk+1 B −1 A)z k .(1.38)Теперь зафиксируем некоторую матрицу D : D = DT > 0. Как уже говорилось, в этом случае1существует такая матрица D 2 , что11D 2 D 2 = D,11D 2 = (D 2 )T .Далее мы будем подбирать τk+1 так, чтобы минимизировать ||z k+1 ||D (считаем, что z k уже задана) — такой способ построения итерационного процесса называется локальной минимизацией.Согласно определению ||z k+1 ||D , имеем:rDqE||z k+1 ||D =hDz k+1 , z k+1 i =11D 2 z k+1 , D 2 z k+1 = ||y k+1 ||,1где y k+1 = D 2 z k+1 .То есть, минимизация ||z k+1 ||D эквивалентна минимизации ||y k+1 ||.
Преобразуем немного выражение для z k+1 в (1.38):11z k+1 = (E − τk+1 B −1 A)D− 2 D 2 z k .1Домножим обе части равенства на D 2 :1D 2 z k+11⇐⇒ D 2 z k+1111= D 2 (E − τk+1 B −1 A)D− 2 D 2 z k ⇐⇒=111(E − τk+1 D 2 B −1 AD− 2 )D 2 z k .Согласно принятым обозначениям, имеем:11y k+1 = (E − τk+1 D 2 B −1 AD− 2 )y k .11Переобозначив C = D 2 B −1 AD− 2 , получим более короткую запись:y k+1 = (E − τk+1 C)y k .Вспомним, что нам надо уменьшать ||y k+1 ||. Более удобно работать с квадратом этого выражения: ||y k+1 ||2 = y k+1 , y k+1 = (E − τk+1 C)y k , (E − τk+1 C)y k .30Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙИз линейности скалярного произведения получаем:2||y k+1 ||2 = ||y k ||2 + τk+1||Cy k ||2 − 2τk+1 y k , Cy k .Потребуем дополнительное условие: положительную определенность C.
В этом случае квадратнормы ||y k+1 || примет удобный для нас вид: #Cy k , y k=||y|| = ||y || + Cy , Cy− 2τk+1hCy k , Cy k i" k k # 2 k k 2 kCy , yCy , yk 2k−.= ||y || + Cy , Cyτk+1 −kkhCy , Cy ihCy k , Cy k i"k+1 2k 2kk2τk+1Очевидно, что наименьшее значение ||y k+1 || достигается, когда второе слагаемое обращается в нуль,то есть, когда: k kCy , yτk+1 =hCy k , Cy k i— это выражение больше нуля, так как C > 0. Теперь преобразуем эту формулу, используя ранние111обозначения (y k = D 2 z k , C = D 2 B −1 AD− 2 ):Dτk+11111D 2 B −1 AD− 2 D 2 z k , D 2 z kEDB −1 Az k , z kE==D 1.11111hDB −1 Az k , B −1 Az k iD 2 B −1 AD− 2 D 2 z k , D 2 B −1 AD− 2 D 2 z k(1.39)Получилось, что для получения оптимального выражения для τk+1 необходимо знать z k . Далеес помощью варьирования матрицы D сведем это к вычислению известных нам величин.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
