Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 8
Текст из файла (страница 8)
Остаётсяоценить погрешность. ИмеемZRn (f ) 6|g(x) − Ln (x)| dx.(45)Стандартным образом оценивая подынтегральное выражение и домножая на длину отрезка, получаем оценкудля Rn (f ).Определение. В случае, когда многочлен Лагранжа берётся по трём узлам, эта формула называется формулой Филона.3.5.4. Оптимальные квадратурыТут были слова про поиск оптимальной квадратуры на заданном классе. Подробности см. БЖК., а в программе экзамена этого, судя по всему, нет.Пусть F — какой-то класс функций, а K некоторый класс квадратур. Рассмотрим Ek (f ) := |I(f ) − K(f )|,где k — квадратура класса K, а f ∈ F .
Наилучшей квадратурой на заданном классе называется квадратураk ∗ := argmin inf sup Ek (f ).k∈K f ∈F(46)Оценки погрешностей квадратур выражаются через производные функций данного класса. Так, например,формула прямоугольников является оптимальной для функций с ограниченной второй производной (теоремаНикольского).4. Численные методы линейной алгебрыВ этой главе мы будем пытаться решать линейные системы вида Ax = b.4.1. Точные методыИмеет место очевидная оценка для сложности алгоритмов решения систем и нахождения собственных значений — O(m2 ), так как даже самый оптимальный алгоритм должен задействовать все элементы матрицы.Имеется и оценка сверху — алгоритм Гаусса даёт оценку 23 m3 операций для решения линейных систем.Классический алгоритм Гаусса ломается довольно быстро, потому что набегает погрешность.
Поэтому обычно используют его модификации (метод Гаусса с выбором главного элемента и т. п.).Кроме метода Гаусса, имеются более качественные методы (например, метод отражений, метод вращений,метод Холецкого). Многие из них лучше обычного метода Гаусса тем, что используют ортогональные преобразования, а это уменьшает погрешность. Естественно, что трудоёмкость при этом возрастает (за всё приходитсяплатить).Про эти методы можно почитать, например, в книжке Богачёва, но мы всё же опишем некоторые из них.4.1.1. Метод отраженийПусть w — вектор единичной длины. Рассмотрим матрицу Uw = E − 2wwt . Изучим её свойства:1.
Uw симметрична.2. Uw2 = E. В самом деле,Uw2 = I − 4wwt + 4w(wt w)wt = E.3. Из свойств 1 и 2 следует, что Uw Uwt = E, то есть матрица Uw ортогональна.27(1)4. Из свойства 2 следует, что Spec Uw ⊂ {1, −1}.5. Uw w = w − 2w(wt w) = −w.6. Если x⊥w, то Uw x = x − wwt x = x − w(w, x) = x, ибо (w, x) = 0.Таким образом, матрица Uw является матрицей отражения относительно гиперплоскости с нормалью w.Пусть заданы два вектора e и y. Найдём вектор w, такой что Uw y = αe (α ∈ R).
Запишем вектор y в видеy = γw + v, где v⊥w. Поскольку Uw не меняет длины векторов, получаем, что α = ± |y||e| . После примененияматрицы отражения получаем систему(−γw + v = αe,(2)γw + v = y,значит,2γw = y − αe.(3)Поскольку w должен быть единичным, получаем, чтоw=y − αe.|y − αe|(4)Теперь уже можно построить алгоритм.
Подбираем вектор w1 так, чтобы после отражения с помощью матрицы U1 := Uw1 первый столбец матрицы A перешёл бы в вектор, пропорциональный e1 (здесь {ei } — стандартныйбазис). После этого переходим к системе U1 Ax = U1 b, и так далее.Тогда в процессе работы алгоритма получаем система будет преобразовываться так:Um−1 Um−2 · . . . · U1 Ax = Um−1 Um−2 · . .
. · U1 b.(5)Отметим, что вычисление числа α требует порядка m операций, вычисление w — ещё порядка m операций.Что касается умножения на матрицу Uw слева, то оно требует порядка m2 операций, если его производитьпоследовательно, то есть дважды умножая матрицу на вектор (а не m3 , как может показаться, если сначалаявно вычислить матрицу Uw , а потом умножать её на матрицу системы!)Преимущество этого метода: его устойчивость. Мы всякий раз совершаем ортогональные преобразования,значит, погрешность будет накапливаться не очень сильно.4.1.2.
Метод ХолецкогоМетод Холецкого, он же метод квадратного корня, основан на построении разложения матрицы в произведение A = SDS t , где D — диагональная матрица, а S — верхнетреугольная.На этом краткий обзор точных методов заканчивается, и мы переходим к итерационным методам.4.2. Итерационные методы4.2.1. Метод простой итерацииОдним из наиболее примитивных является метод простой итерации, который состоит в преобразованиисистемы Ax = b к виду x = Bx + c и применении итераций xn+1 = Bxn + c.Если kBk < 1, то данный процесс, очевидно, сходится со скоростью геометрической прогрессии.
В самомделе, пусть x — точное решение то есть x = xn − rn . Имеем(x = Bx + c,(6)xn+1 = Bxn + c.Вычитая уравнения друг из друга, получаем krn+1 = Brn k → 0 при n → ∞.В дальнейшем мы будем рассматривать следующие матричные нормы:P• kAk∞ = max |aij |,jPi• kAk1 = max |aij |,ijp• kAk2 = max λ(At A).Чтобы доказать сходимость метода, хотя бы в одной из этих норм должно быть выполнено kBk < 1.Теорема 4.1 (Критерий сходимости метода простой итерации). Пусть λi — собственные значенияматрицы B. Метод простой итерации сходится тогда и только тогда, когда |λi | < 1 при всех i.
Докажем сначала обратное утверждение. Нам потребуется28Лемма 4.2. Для всякой матрицы B и всякого η > 0 существует невырожденная матрица T , такая чтоB = T −1 ΛT,(7)где Λ имеет вид матрицы из жордановых клеток, у которых над главной диагональю вместо единиц стоитчисло η (а на диагонали, как обычно, стоят собственные значения матрицы B).
Потом напишем. Применить теорему о ЖНФ и подкрутить. Пусть |λi | < 1 при всех i. Их конечное число, поэтому можно выбрать столь малое η > 0, что |λi | + η < 1 привсех i. Тогда имеем kΛk1 < 1. Покажем, что метод сойдётся:rn = B n r0 = (T −1 ΛT )n r0 = T −1 Λn T r0 ,(8)но kΛn k1 → 0, поэтому и rn → 0.Теперь докажем прямое утверждение. От противного. Допустим, что существует собственное значение λ,такое что |λ| > 1, и пусть вектор e отвечает собственному значению λ. Выберем подходящее начальное условие,такое что x − x0 = e, то есть r0 = e.
Тогда rn = B n r0 = λn e 9 0, что и требовалось доказать. Пример 2.1. Это будет не просто пример, а даже контрпример, который показывает, что в общем случаевсё плохо. Пусть n < m (m судя по всему, размерность пространства).Пусть B = αE + βJ(0), где J(λ) — жорданова клетка с собственным значением λ, а α и β — некоторые числа.Тогда kB n k1 = (|α| + |β|)n . В самом деле, при n < m матрица будет иметь вид nαC1n αn−1 β .
. . β n 0 . . .00αn. . . . . . . . . . . . . . . . . . . . .Bn = (9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ,n0 ....................... 0 αто есть она будет верхнетреугольной, а по строкам будут биномиальные слагаемые. Ясно, что максимум приподсчёте нормы будет достигаться в первой строке, а сумма модулей чисел в ней — это как раз (|α| + |β|)n .Пусть c := (0, 0, . . . , 0, 1). Несложно видеть, что итерации отображения B имеют вид!n−1 n−2βββ1,1xn =0, . . . , 0,,,...,(10)1−α1−α1−α1−α(это проще всего проверить по индукции). β Пусть числа α и β таковы, что |α| < 1, α < 0, |β| > 1, и при этом |α| + |β| > 1 и 1−α < 1.
Отметим, что притаком выборе kBk1 растёт с экспоненциальной скоростью. Оценим норму итерации решения:kxn k1 =βn(1−α)nβ− 1−α1−161−α.1 − (α + β)(11)Казалось бы, ничего особенного. Положим r0 := (0, 0, . . . , 0, 1). А теперь начинаем считать погрешность:rn = B n r0 = (0, . . .
, 0, β n , . . . , C1n αn−1 β, αn ).(12)Если размерность достаточно большая, то погрешность сначала будет экспоненциально расти, так как β n → ∞.Через несколько итераций это прекратится, так как β n «выползет за пределы» вектора (количество нулевыхкоординат у вектора xn уменьшается с ростом n). Но та ошибка, которая успеет набежать за это время, ужеуспеет всё испортить.Чтобы поверить в то, что ни одна машина с такой системойне справится, рекомендуется прогнать этот метод β 5 2157на таких параметрах: α = − 2 , β = 4 . Тогда |α| + |β| = 4 , а 1−α = 4 · 3 = 56 < 1, то есть все условия выполнены.Отметим, что сама норма итераций решения, как видно из оценки, не превосходитm ∼ 100, n ∼ 50, чтобы всё сломалось.321+ 21 − 54= 6. Достаточно взятьТут был замечен такой бред: на 50 итерациях β n вырастет сильно больше, чем оценка нормы решения...294.2.2.
Модификация метода простой итерации (метод Ричардсона)Из приведённого в предыдущем разделе примера видно, что просто так применять описанный метод простойитерации нельзя. Будем пытаться его улучшить, исходя из того, что мы обладаем некоторой информацией оспектре исходного оператора A.Пусть мы знаем, что A = At и Spec(A) ⊂ [µ, M ], причём µ > 0.Перепишем наше уравнение в видеx = x − τ (Ax − b).(13)Здесь τ 6= 0 — некоторое число, которое мы сейчас будем подбирать. Это самый обычный метод простой итерации: достаточно переписать это уравнение в видеx = (E − τ A) x + |{z}τb .| {z }(14)xn+1 = xn − τ (Axn − b).(15)cBТеперь сделаем из нашего уравнения итерационный процесс:Если переписать его в видеxn+1 − xn+ Axn = b,(16)τполучим, что это разностный аналог параболического уравнения.Когда он сходится? Имеем λB = 1 − τ λA , поэтому, применяя теорему о сходимости метода простой итерации,получаем, что сходимость будет тогда и только тогда, когда(τ λA > 0,2⇔ 0<τ <|1 − τ λA | < 1 ⇔.(17)2Mτ< λ .A2А теперь поставим себе целью найти оптимальное τopt ∈ 0, M, для которого метод сходится быстрее всего.Фактически, нас интересует τ = arg min kBkевкл , а в силу симметричности матрицы B получаем, чтоτ(18)τopt = arg min max |1 − τ λ| .τλ∈[µ,M]Поскольку функции кусочно-линейные, максимум достигается в концах отрезка, стало быть,(19)τopt = arg min max {|1 − τ µ| , |1 − τ M |} .τДля наглядности построим графики функций f1 (τ ) := |1 − τ µ| и f2 (τ ) := |1 − τ M |.1Mτopt1µτ2При τopt = M+µимеем минимум нормы kBk = M−µM+µ .
Предположим, что t :=вать как малый параметр и раскладывать по нему в ряд. ОтсюдаµM≪ 1. Тогда t можно рассматри-M −µ1−t== (1 − t)(1 − t + t2 − t3 + . . . ) = 1 − 2(t − t2 + . . . )M +µ1+t(20)Мы хотим, чтобы |rn | 6 kBkn |r0 | 6 ε|r0 |, то есть n ln kBk 6 ln ε. Таким образом, число итераций длядостижения точности ε должно быть не меньше, чемn>ln ε−1ln ε−11M≈=ln ε−1 .− ln kBk2t2 µНа лекциях тут, судя по всему, была допущена ошибка: потерян множитель301,2который возникает из разложения по t.(21)4.2.3.
Upgrade метода Ричардсона, или чебышевское ускорениеСудя по всему, этот текст соответствует билету про чебышевское ускорение. На лекциях этих слов, кажется, не произносилось,но многочлены Чебышёва присутствуют достаточно явно.µНо нам этого мало. Мы хотим уменьшить зависимость параметра n от отношения M. Будем старатьсяподбирать τ так, чтобы норма оператора перехода от 1 шага к n-му была минимальной. На каждом шаге будемиспользовать поправку τk . Имеем(22)rn = (E − τn−1 A) · . .
. · (E − τ0 A) r0 .{z}|Pn (A)На r0 действует некоторый многочлен Pn от оператора A. Имеем Pn (0) = 1. Так как матрица A симметрична,то по лемме об отображении спектра имеемkPn (A)k =(23)max |Pn (λ)|.λ∈Spec AВпрочем, у нас нет полной информации о спектре, а мы знаем только, на каком отрезке он лежит.
Поэтому мысможем лишь оценить сверху норму этого оператора:(24)max |Pn (λ)| 6 kPn (λ)kC[µ,M] .λ∈Spec AСреди всех многочленов найдём многочлен с минимальной нормой. Разумеется, это многочлен Чебышёвастепени n на отрезке [µ, M ]. Автоматически получаем n корней и n значений параметров τk . Пусть tk — нулистандартного многочлена Чебышёва. Тогда имеемtk :=откуда(2k + 1)πM + µ − 2λk= cos,M −µ2n(25)M +µ M −µ(2k + 1)π1−cos, τk =.(26)222nλkОценим сверху нормировочный множитель t1n , где tn := Tn M+µM−µ .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
