Численные методы. Соснин (2005) (1160462), страница 2
Текст из файла (страница 2)
+ li(k−1) l(k−1)k + lik lkk = aik ⇐⇒⇐⇒ li1 lk1 + li2 lk2 + . . . + li(k−1) lk(k−1) + lik lkk = aik ,i = k, n.(1.4)Пусть i = k, тогда из (1.4) получаем:2222+ lkklk1+ lk2+ . . . + lk(k−1)= akk .2Переносим все слагаемые, кроме lkk, в правую часть и извлекаем корень:lkk =q2 − . .
. − l2akk − lk1k(k−1) .(1.5)Можно показать, что под корнем получается отношение положительного минора k-го порядка, кположительной сумме:s|...| > 0lkk =.(. . .) > 0Подставив (1.5) в (1.4), получим:lik =aik − li1 lk1 − . . . − li(k−1) lk(k−1),lkki = k + 1, n.Мы полностью определили элементы матрицы L, доказав тем самым работоспособность нашегометода. Заметим, что мы использовали симметричность матрицы A, когда брали в уравнениях длянахождения lij только те ее элементы, которые лежат на главной диагонали или ниже ее.Модифицированный метод квадратного корняТеперь несколько модифицируем метод квадратного корня, а точнее, распространим его на болееширокий класс задач. Итак, мы снова рассматриваем уравнениеAx = fпри условии, что det A 6= 0.(1.6)Предположим, что матрица A симметрична (A = AT ) — положительной определенности требоватьне будем.
Покажем, что для таких матриц справедливо представлениеA = LDLT ,(1.7)где L — нижнетреугольная матрица с единицами на главной диагонали, а D — диагональная матрица:L=10...lij... 0.. .;... 01D=d110...00............00.. . .0 dnn8Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙДля наглядности (1.7) можно представить так (LT = (lij )):10...lij... 0d11.. . 0 .... 0 ..010............010.. .... 0 .0 ..dnn0 ...lij = A.1...0Обозначим произведение матриц DLT за новую матрицу B = (bij ) и заметим, что lij = lji .Теперь найдем элементы матриц L и D.
Согласно определению произведения матриц, элемент aikпри i > k определяется так:aik = li1 b1k + li2 b2k + . . . + li(k−1) b(k−1)k + lik bkk— мы отбросили нулевые последние слагаемые, обусловленные видом матрицы L.А элемент aik при i = k определяется так:akk = lk1 b1k + lk2 b2k + . . . + lk(k−1) b(k−1)k + lkk bkk .Теперь воспользуемся тем, чтоbjk = djj ljk = djj lkj ,и перепишем два последних уравнения:aik = li1 (d11 lk1 ) + li2 (d22 lk2 ) + .
. . + li(k−1) (d(k−1)(k−1) lk(k−1) ) + lik bkk ;(1.8)akk = lk1 (d11 lk1 ) + lk2 (d22 lk2 ) + . . . + lk(k−1) (d(k−1)(k−1) lk(k−1) ) + lkk bkk .(1.9)При k = 1 (1.9) выглядит так:2l11d11 = a11 .Зная, что l11 = 1, получим формулу для d11 :d11 = a11 .Теперь рассмотрим случай k = 1, i = 2, n. Из (1.8) получим:li1 (d11 l11 ) = ai1=⇒{l11 = 1}=⇒li1 =ai1.d11Мы получили уравнения, определяющие первый столбец матрицы L.Аналогично, при k = 2 (1.9) и (1.8) перепишутся так:l21 (d11 l21 ) + l22 (d22 l22 ) = a22li1 (d11 l21 ) + li2 (d22 l22 ) = ai2=⇒=⇒{l22 = 1}li2 ==⇒2d22 = a22 − l21d11 ;ai2 − li1 d11 l21.d22И так далее. Вот как будут выглядеть общие формулы для элементов матриц D и L :22dkk = akk − (lk1d11 + .
. . + lk(k−1)d(k−1)(k−1) );lik =(aik − li1 d11 lk1 − . . . − li(k−1) d(k−1)(k−1) lk(k−1) ).dkkЗаметим, что формулы корректны — элементы dkk не могут быть равны нулю, так как в этомслучае определитель матрицы A обращался бы в ноль.1.2. ЛИНЕЙНЫЕ ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ9Итак, разложение (1.7) получено. Подставим его в исходное уравнение (1.6):LDLT x = f=⇒{обозначим DLT x = y}=⇒Ly = f.Мы получили уравнение относительно вектора y. В левой части стоит нижнетреугольная матрицаL, что позволяет нам сильно сократить вычисления, применяя только обратный ход метода Гаусса.После нахождения y вспомним, как мы его определяли:DLT x = y=⇒{обозначим LT x = z}=⇒Dz = y— это система уравнений, определяемая диагональной матрицей D.
Решая ее, получим z.Теперь у нас есть простая система для нахождения x :LT x = z.Наличие в левой части верхнетреугольной матрицы L, опять же, сильно упрощает вычисления.Сложив количество операций, необходимых для получения разложения (1.7) и для поиска x, по3лучим сложность данного метода, равную ≈ n6 .Итерационные методы решения СЛАУЭто следующая группа методов поиска решений систем линейных уравнений, особенно эффективнаяпри решении больших систем (с числом неизвестных порядка тысячи и более).В общем случае сначала задаются некоторым вектором x0 , называемым начальным приближением. От него строится последовательность x1 , x2 , .
. . , xk и так далее, где число k называютномером итерации. В общем случае (k + 1)-е приближение зависит от всех предыдущих:xk+1 = Fk+1 (x0 , x1 , . . . , xk ).От последовательности, естественно, ожидается сходимость к вектору x, который будет являтьсярешением исходной системы.Определение. Итерационный метод называется m−шаговым, если каждое последующее итерационное приближение строится лишь по m предыдущим:xk+1 = Fk+1 (xk−m+1 , . .
. , xk−1 , xk )— на практике наиболее часто используется m = 1 и m = 2.Определение. Если Fk+1 — линейная функция, то соответствующий итерационный метод называется линейным.1.2Линейные одношаговые итерационные методыСогласно определению, общее выражение для xk+1 в линейных одношаговых итерационных методахтаково:xk+1 = S k+1 xk + ψk+1 ,(1.10)где S k+1 — некоторая матрица, а ψk+1 — некоторый вектор.Логично потребовать, чтобы вектор x = A−1 f (то есть искомое решение) при подстановке вместоxk+1 и xk обращал (1.10) в тождество:A−1 f = S k+1 A−1 f + ψk+1 ,10Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙверное для всех k.Преобразовав, получим:(A−1 − S k+1 A−1 )f = ψk+1=⇒ψk+1 = Qk+1 f,где Qk+1 = A−1 − S k+1 A−1 .
Потребуем, чтобы эта матрица была обратима, и формулу для Qk+1домножим на A :S k+1 = E − Qk+1 A.Это выражение подставим в (1.10):xk+1 = xk − Qk+1 Axk + Qk+1 f.Домножим это выражение на (Qk+1 )−1 и перенесем слагаемые:Qk+1−1(xk+1 − xk ) + Axk = f.Теперь домножим и разделим первое слагаемое на некоторую константу τk+1 :Qk+1−1τk+1xk+1 − xk+ Axk = f.τk+1−1Обозначив Bk+1 = Qk+1τk+1 , получим так называемую каноническую форму записи одношагового итерационного метода:xk+1 − xk+ Axk = f.(1.11)Bk+1τk+1Примечание. Мы могли бы и не вводить константу τk+1 , однако впоследствии будет видно, чтос матрицей Bk+1 и константой τk+1 удобнее работать по отдельности.Примечание.
Далее мы иногда будем отождествлять итерационный метод и его каноническуюформу — так что не пугайтесь.Примечание. Матрица B обязана быть обратимой — иначе вывод канонической формы записиитерационного метода будет некорректен.Определение. Если Bk+1 = E, то соответствующий метод называется явным (находим xk+1 , нерешая системы уравнений), в противном случае — неявным.Определение. Если Bk+1 = B и τk+1 = τ, то метод называется стационарным, в противномслучае — нестационарным.Определение. Вектор z k = xk − x (отклонение от точного решения) называется погрешностьюна k-й итерации.k→∞Определение. Итерационный метод называется сходящимся, если kz k k −→ 0 для некоторойвыбранной нормы || · ||.Понятно, что предел может не быть достигнут за конечное число шагов.
В таком случае мы задаемся некоторой точностью ε — достаточно малым числом ( ∼ 10−3 , 10−5 ) — и будем останавливатьитерационный процесс при выполнении неравенства||xk − x|| 6||x0 − x||1ε,или ||xk − x|| 6 ε||x0 − x||для некоторого k. В этом случае мы будем говорить, что решение получено с точностью ε.1.2.
ЛИНЕЙНЫЕ ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ11Кроме того, естественно потребовать, что если это неравенство выполняется для k0 = k0 (ε), то онодолжно выполняться и для любого k > k0 (ε).Число k0 (ε) называется минимальным числом итераций, необходимых для достижениизаданной точности ε. Понятно, что чем меньше k0 (ε), тем лучше метод.Примеры одношаговых итерационных методовДля иллюстрации нижеприведенных методов сначала более детально распишем уравнение Ax = f :a11 x1 + a12 x2 + . . .
+ a1n xn = f1 ; a21 x1 + a22 x2 + . . . + a2n xn = f2 ;. . . . . . . . . . . . . . . . . . . . .an1 x1 + an2 x2 + . . . + ann xn = fn .Расставляя итерационные индексы над компонентами вектора X, можно получить различные итерационные процессы.Метод ЯкобиРасставим индекс следующего приближения (k + 1) над диагональными элементами, а индекс k —над остальными:+ a12 xk2 + . . . + a1n xkn = f1 ; a11 xk+11 a21 xk + a22 xk+1 + .
. . + a2n xk = f2 ;1n2. . . . . . . . .an1 xk1+ an2 xk2= fn .+ . . . + ann xk+1nИз этой системы легко получить общее уравнение для xk+1:ifi −xk+1=ii−1Xaij xkj −j=1nXaij xkjj=i+1,i = 1, n.(1.12)aiiЧтобы от этой расчетной формулы перейти к канонической форме записи метода, представим матрицу A в виде суммы трех матриц — нижнетреугольной, диагональной, и верхнетреугольной:A1 = 0...0aij..;.0A = A1 + D + A2 ,a11... 0;D=...0annA2 = 0...0aij....0(1.13)Из (1.12) легко вывести, чтоA1 xk + Dxk+1 + A2 xk = f.Для приведения к канонической форме (1.11) прибавим и вычтем Dxk :D(xk+1 − xk ) + Axk = f.Отсюда получаем формальное определение метода Якоби через матрицу Bk+1 и константу τk+1 :Bk+1 = D;τk+1 = 1.Легко видеть, что он неявный и стационарный.12Глава 1.
РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙМетод Зейделя2Пусть теперь итерационные индексы (k + 1) стоят не только на диагонали, но и под ней:a11 xk+1+ a12 xk2+ . . . + a1n xkn= f1 ;1 a21 xk+1 + a22 xk+1 + . . . + a2n xk= f2 ;n12. . . . . . . . . .an1 xk+11+ an2 xk+12+ . . . + ann xk+1n= fn .Относительно xk+1 эта система решается так: сначала из первого уравнения находится xk+1, потом1из второго — xk+1ит.д.2Для записи метода в канонической форме представим матрицу A в виде суммы матриц A1 , D, A2(как и в методе Якоби):A = A1 + D + A2 ,где A1 — нижнетреугольная, D — диагональная, а A2 — верхнетреугольная матрица.
Тогда написанная выше итерационная схема будет выглядеть так:A1 xk+1 + Dxk+1 + A2 xk = f.Преобразуя это уравнение, приведем его к каноническому виду:(A1 + D)xk+1 + A2 xk = f⇐⇒⇐⇒(A1 + D)(xk+1 − xk ) + Axk = f.Отсюда получаем формальное определение метода Зейделя:Bk+1 = A1 + D;τk+1 = 1.Как несложно заметить, он является неявным и стационарным.Метод релаксацииФормально определим метод следующим образом:Bk+1 = D + ωA1 ;τk+1 = ω,где ω — некий числовой параметр. Нетрудно заметить, что метод Зейделя является частным случаемэтого метода при ω = 1.Кроме того, метод релаксации является стационарным и неявным.Метод простой итерацииИ здесь никакой подробной системы не будет:Bk+1τk+1= E;= τ.Подставив эти равенства в (1.11), получим каноническую форму метода простой итерации:xk+1 − xk+ Axk = f.τНа этот раз метод является явным и стационарным.2 ПредложенЛ.
Зейделем в 1874 г.1.3. СХОДИМОСТЬ ОДНОШАГОВЫХ СТАЦИОНАРНЫХ МЕТОДОВ13Метод РичардсонаВ данном случае матрица Bk+1 — единичная, а константа τk+1 определяется по некоторым расчетным формулам. Выражение для τk+1 приводить не будем, запишем только каноническую формузаписи метода:xk+1 − xk+ Axk = f.τk+11.3Сходимость одношаговых стационарных методовТеорема 1.1 (Достаточное условие сходимости одношагового стационарного итерационного метода).Пусть матрица A — симметрическая положительно определенная, B — положительно определена,(xk+1 − xk )и τ > 0, тогда итерационный процесс B+ Axk = f сходится для любого начальногоττприближения x0 , если B − A > 0.2Доказательство этой теоремы было дано в курсе «Введение в численные методы», кроме того, егоможно найти в [1].Утверждение 1.1. Пусть матрица A — симметрическая положительно определенная, и выполненоnXусловие диагонального преобладания: aii >|aij |.
Тогда метод Якоби является сходящимся дляj=1i6=jлюбого начального приближения.Доказательство. Так как матрица A — положительно определена, ∀ x 6= 0Распишем скалярное произведение поподробнее:hAx, xi =n XnXaij xi xj 6i=1 j=1612n XnXn XnXi=1 j=1|aij |x2i +i=1 j=1|aij ||xi xj | 6hAx, xi > 0.x2i + x2j|xi xj | 662nnn XnnnXXX1 XX|aji |x2j = {i ↔ j} =|aij |x2i =x2i |aii | +|aij | <2 i=1 j=1i=1 j=1i=1j=1i6=j< {так как матрица A — матрица с диагональным преобладанием} <nnXX< 2x2i |aii | = 2|aii |xi xi = 2 hDx, xi .i=1i=1Таким образом, получаем:0 < hAx, xi < 2 hDx, xi .Откуда следует, что2 hDx, xi − hAx, xi > 0=⇒2D − A > 0.То есть матрица (2D − A) положительно определена. В этом случае выполняются все условиятеоремы 1.1 при B = D и τ = 1 — это и дает сходимость метода Якоби.Утверждение 1.2 (Сходимость метода релаксации).