Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 9
Текст из файла (страница 9)
А для его вычисления нам нужно будетвспомнить явную формулу для многочленов Чебышёва:n n ipp1 hTn (x) =x + x2 + 1 + x − x2 + 1.(27)2λk =Вычислим tn по этой формуле:1tn =2" √√√n n √ !nM+ µM + µ + 2 MµM + µ − 2 Mµ1√++=√M −µM −µ2M− µ√√ !n #M− µ√.√M+ µ(28)Далее у Кобелькова была ошибка! Он не в ту сторону оценил, поэтому получался полный бред. А потом он ещё раз ошибся, исработало свойство чётности ошибок. А надо-то вот как:Второе слагаемое мы просто откинем (оно мало́ по сравнению с первым).
Обозначив t :=162tnИтак, мы считаем по формулеполучаем√ n1+ t√.1− t(29)√ n√ n1− t√≈2 1−2 t .1+ t(30)tn >ОтсюдаµM,1rn =· Tntn12(M + µ)I − 2AM −µr0 ,то есть в среднем на n шагах знаменатель геометрической прогрессии составляетно вывести, что для достижения точности ε нужно сделать порядка√n ≈ t ln ε−131(31)√√n2 1 − 2 t . Отсюда неслож(32)итераций. При этом общее число шагов должно быть кратно n.
В этом большой дефект метода. В реальностиметод работает плохо, для улучшения его нужно переставлять корни.Можно использовать следующий алгоритм, предложенный в 1969 году Лебедевым (учеником Соболева) дляперестановки корней при n = 2m .
Пусть имеется некоторая нумерация корней на (m−1)-м шаге (здесь s := 2m−1 )(b0 , b1 , . . . , bs , s + b0 , . . . , s + bs ).(33)Эта последовательность переходит в последовательность(s + 1 + b0 , b0 , s + 1 − b1 , b1 , s + 1 − b2 , b2 , . . . ).(34)Оказывается, в этом случае нормы промежуточных многочленов при вычислении ограничены константой 1.4.2.4. Линейный оптимальный процессНедостаток метода Ричардсона заключался в том, что нужно было идти шагами по n. А мы хотим метод,который давал бы оптимальную аппроксимацию на каждом шаге при любом значении n. Оптимальным мыназовём такой процесс, при котором на каждом шаге имеет место равенство(M + µ)I − 2A1M +µrn =· Tn.(35)r0 , tn := TntnM −µM −µ|{z}pНапишем три последовательных шага метода.
Получим соотношенияtn−1 rn−1 = Tn−1 (p)r0 ,tn rn= Tn (p)r0 ,tn+1 rn+1 = Tn+1 (p)r0 ,(36)Вспомним, что для Tn выполнены рекуррентные соотношенияTn−1 (p) − 2pTn (p) + Tn+1 (p) = 0.(37)Умножим второе уравнение на −2p и сложим. Получаем соотношениеtn+1 rn+1 − 2tn4tnM +µrn +Arn + tn−1 rn−1 = 0.M −µM −µ(38)Так как rn = xn − x, то, в силу равенстваtn+1 − 2tnполучаемtn+1 xn+1 − 2tnM +µ+ tn−1 = 0,M −µM +µ4tnxn +(Axn − b) + tn−1 xn−1 = 0.M −µM −µ(39)(40)Отсюда получаем процедуру оптимального поиска решения, которая использует только три «слоя».Осталось сделать данную процедуру устойчивой. В силу равенства (39) можно переписать соотношение так:tn+1 xn+1 − (tn+1 + tn−1 )xn +2(tn+1 + tn )(Axn − b) + tn−1 xn−1 = 0.M +µ(41)tnПоделим всё на tn+1 и введём обозначение ωn := − tn+1. Тогдаxn+1 − (1 + ωn−1 ωn )xn +2(1 + ωn−1 ωn )(Axn − b) + ωn−1 ωn xn−1 = 0.M +µ(42)Имеем ω0 = − tt10 = − M−µM+µ .
Очевидно, |ω0 | < 1 по построению.Разделив соотношение (39) на tn+1 , получаем1 + 2ωnM +µ+ ωn−1 ωn = 0,M −µ32(43)откудаωn = −M+µ2 M−µ1,+ ωn−1(44)а значит, |ωn | < δ(M, µ) < 1, и потому процесс сходится.Чтобы сделать первый шаг, нужно вычислить x1 по формулеx1 = x0 −2(Ax0 − b).M +µ(45)Среди недостатков этого процесса можно отметить то, что в нём используется информация о спектре матрицы.
В реальном мире спектр далеко не всегда известен даже приблизительно.4.3. Другие методы4.3.1. Метод скорейшего спускаМы будем рассматривать симметричные положительно определённые матрицы A. Рассмотрим функционал(46)F (x) := (Ax, x) − (2b, x).Пусть x — точное решение, то есть Ax = b. Заметим, чтоA(x − x), x − x = (Ax, x) − 2(b, x) + (Ax, x) = F (x) + (Ax, x).(47)Поскольку A — положительно определённая матрица, имеем A(x − x), x − x > 0, причём равенство нулюдостигается тогда и только тогда, когда x = x.
Значит,arg min F (x) = x.(48)kxk2A := (Ax, x).(49)xn+1 = xn − αn grad F (xn ).(50)Введём норму, задаваемую матрицей A:Метод скорейшего спуска устроен следующим образом:Несложно видеть, что (здесь верхние индексы обозначают координаты векторов)∂F= 2(Ax − b)k .∂xk(51)Обозначим ξn := Axn − b, δn := 2αn . С учётом этого обозначения, задача переписывается в виде(52)xn+1 = xn − δn ξn .Итак, нам нужно найти min функционала F по δn :F (xn+1 ) = A(xn − δn ξn ), xn − δn ξn − 2(b, xn − δn ξn ) == (Axn , xn ) − 2δn (Axn , ξn ) + δn2 (Aξn , ξn ) − 2(b, xn ) + 2δn (b, ξn ).(53)Дифференцируем по δn , получаем222F ′ (xn+1 ) = 2δn kξn kA − 2(Axn − b, ξn ) = 2δ kξn kA − 2 kξn k .(54)Приравнивая эту производную, как и положено, к нулю, получаемРассмотримy := xn −kξn k2.(55)2(Axn − b).M +µ(56)δn =2kξn kAЭто шаг в направлении градиента F .
Поскольку при вычислении xn+1 мы выбирали параметр оптимально, апри вычислении y — рассмотрели некоторое фиксированное его значение, то во всяком случае имеемF (xn+1 ) 6 F (y).33(57)Теперь можно оценить погрешность, причём нам будет удобно оценивать её в норме k·kA . Из определения Fследует, что22F (x) = kx − xkA − kxkA .(58)Заметим, что имеет место равенствоy − x = xn − x −Отсюдаr := y − x =2A(xn − x).M +µ(59)2E−A rn .M +µ(60)2A Arn .M +µ(61)Применим к последнему равенству матрицу A:Ar =E−Применим известное утверждение из линейной алгебры: для всякого самосопряжённогооператора существуPет ортогональный собственный базис {e1 , .
. . , em }, и пусть ei ∼ λi . Пусть rn = ci ei , тогдаXX2λi2λir=1−ci e i ,Ar =1−ci λi ei ,(62)M +µM +µоткуда2krkA = (Ar, r) =X1−2λiM +µ2c2i λi .(63)Заменяя λi в скобках на самое маленькое, что там может быть, а именно µ, получаем оценку2krkA 6следовательно,M −µM +µ2 XkrkA 6c2i λi =M −µM +µ22krn kA ,M −µkrn kA .M +µ(64)(65)Из формул (57) и (58) следует, чтоkxn+1 − xkA 6 ky − xkA ,(66)krn+1 kA 6 krkA .(67)то естьКомбинируя это неравенство с неравенством (65), получаемkrn+1 kA 6M −µkrn kA ,M +µ(68)то есть имеет место сходимость со скоростью геометрической прогрессии.
То, что мы оценили сходимость вкакой-то левой норме, нас не смущает, потому что пространство конечномерно и все нормы эквивалентны.У этого метода имеется один недостаток: на каждом шаге приходится вычислять ξn , то есть умножать матрицу на вектор. Кроме того, при вычислении δn опять придётся умножать матрицу на вектор. Этот недостатокможно вылечить следующим образом: вычислять ξn+1 = ξn − δn Aξn , но при этом тащить за собой не толькосам вектор невязки ξn , но и его произведение на матрицу, то есть вектор Aξn .
Тащить ещё один вектор — этонебольшая потеря памяти, поскольку памяти под саму матрицу A всё равно нужно на порядок больше.4.3.2. Метод Ричардсона для несимметричных матрицЗдесь мы рассмотрим итерационный метод, который нам уже встречался. НО теперь мы будем решатьуравнение Ax = b с несимметричной матрицей.Теорема 4.3. Пусть As := 21 (A + At ) > 0.
Тогда существует τ0 > 0, такое что для всех τ ∈ (0, τ0 ) методxn+1 − xn+ Axn = bτсходится.34(69) Мы уже знаем, что со всякой симметричной положительно определённой матрицей B можно связать2норму kxkB := (Bx, x).Пусть x — точное решение, rn — погрешность на n-м шаге. Тогда наш метод перепишется в видеrn+1 − rn+ Arn = 0.τ(70)Введём обозначения r := rn , rb := rn+1 , rt := rb−rb = r + τ rt .τ . Тогда имеет место тождество rВ новых обозначениях метод имеет вид rt + Ar = 0. Пусть Aa — кососимметрическая часть матрицы, то естьAa := 21 (A − At ). Тогда A = As + Aa , следовательно, метод можно переписать так:rt + As r + Aa r = 0,(71)(E − τ As ) rt + As rb + Aa r = 0.| {z }(72)или, подставляя r во втором слагаемом,BОчевидно, матрица B симметрична. Кроме того, ясно, что при всех достаточно малых τ матрица B будетположительно определённой. Будем считать, что τ0 уже выбрано так, чтобы при τ ∈ (0, τ0 ) имеем B > 0.τb = τ (br + r) + τ 2 rt .
ОтсюдаИмеет место тождество rb = rb+r2 + 2 rt , поэтому 2τ r(Brt , 2τ rb) = Brt , τ (br + r) + (Brt , τ 2 rt ) = B(br − r), rb + r + τ 2 (Brt , rt ) =222= (Bbr , rb) − (Br, rb) + (Bbr , r) − (Br, r) + τ 2 (Brt , rt ) = kbr kB − rbB+ τ 2 krt kB . (73)Теперь умножим скалярно равенство (72) на 2τ rb, получим2222kbrkB − krkB + τ 2 krt kB + 2τ kbr kAs + 2τ (Aa r, rb) = 0.(74)(Aa r, rb) = (Aa rb, rb) − τ (Aa rt , rb) = −τ (Aa rt , rb).(75)В силу кососимметричности матрицы Aa имеем (Aa rb, rb) = 0, поэтомуПоэтому можно переписать равенство, написанное выше, следующим образом:2222kbr kB − krkB + τ 2 krt kB + 2τ kbr kAs − 2τ 2 (Aa rt , rb) = 0.(76)Все операторы у нас ограниченные, поэтому найдутся такие константы C1 и C2 , что22kAa yk 6 C1 kykB ,22kyk 6 C2 kykAs .(77)Нам потребуется неравенство 2 |ab| 6 εa2 + 1ε b2 . Оценим скалярное произведение с помощью этого неравенства:2поэтому2τ 2 |(Aa rt , rb)| 6 τ 2 ε kAa rt k +2kbr kB−2krkB2+ τ (1 −τ2C2 τ 2222kbr k 6 C1 ετ 2 krt kB +kbr k As ,εε2C1 ε) krt kBC2 τ2+τ 2−kbrkAs 6 0.ε(78)(79)Поскольку ε можно брать любым, положим ε := C11 , тогда третье слагаемое умрёт.
Кроме того, уменьшая принеобходимости значение τ0 , можно добиться того, что 2 − Cε2 τ = 2 − C1 C2 τ > 1 при τ ∈ (0, τ0 ). После всего этогошаманства получаем неравенство222kbr kB + τ kbr kAs 6 krkB .(80)В силу эквивалентности норм найдётся константа C3 , для которой2Следовательно, имеем неравенството есть, в старых обозначениях,2kbr kAs > C3 kbr kB .(81)(1 + C3 τ ) kbr k2B 6 krk2B ,(82)12krn kB ,1 + C3 τоткуда следует сходимость метода со скоростью геометрической прогрессии. 2krn+1 kB 6(83)Тут ещё были какие-то слова про чебышёвское ускорение и про то, что если спектр лежит в эллипсе, но всё равно ничего непонятно.354.3.3. Метод решения симметричных плохо обусловленных системНазвание придумано самостоятельно. Суть метода состоит в введении переобуславливателя. Во всяком случае лучше соответствует тексту, чем то, что написано в билетах.
Судя по всему, эквивалентно билету про «итерационные методы со спектральноэквивалентными операторами».А ещё в книжке на странице примерно 301 транспонирование у матрицы D (см. ниже) всё-таки забыто не было. Если оно иправда должно там быть (а это правдоподобно), тогда всё становится на свои места.Пусть A — симметрическая положительно определённая (но плохая) матрица, и пусть мы решаем системуметодом итераций:xn+1 − xn+ Axn = b.(84)τМы знаем, что если Mµ ≫ 1, тогда при реальных вычислениях всё разлетится. Будем это лечить. Пусть B —симметрическая положительно определённая матрица, причём такая, чтобы система By = c легко решалась.Заменим исходную систему на такую систему:xn+1 − xn+ Axn = b.(85)Bτ√Большой беды в этом нет, потому что мы легко сможем вернуться к старой системе.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
