Численные методы. Ионкин (2012) (неоффициальные) (косяки есть) (1160437), страница 2
Текст из файла (страница 2)
А так как необходимо решить m систем, то на решение всех систем понадоm3 − m3бится m действий. Факторизация, проведенная один раз, требуетдействий.34mОткуда получаем, что общее количество действий равно m3 − .33Однако, если воспользоваться спецификой видов матриц, то число действийможно уменьшить до m3 . Докажем данное утверждение.Рассмотрим систему (3). Учитывая нижнетреугольную форму матрицы B, получим:(j)b11 y1 = 0y1 = 0b21 y1 + b22 y2 = 0 =⇒···y2 = 0(j)(j)(j)(j)=⇒(j)jbj−1,1 y1 + bj−1,2 y2 + · · · + bj−1,j−1 yj−1=0=⇒(j)(j)yj−1 = 010(j)Откуда получаем, что yi = 0, где i 6 j − 1Следующее j-ое уравнение даст нам (учитывая, что bjj 6= 0):(j)bjj yj = 1(j)=⇒yj =1bjj(∗)Оставшиеся уравнения системы имеют вид:−(j)bi,j yj+(j)bi,j+1 yj+1+ ··· +bi,i yij= 0, i = j + 1, m, bii 6= 0=⇒(j)yi=i−1Pbil yljl=jbiiФиксируем все индексы, т.е.
i и j и считаем число действий:1 деление(i − j) умножений.+Теперь отпускаем один из индексов, например, i. Тогда получим(m − j) + (m − j − 1) + · · · + 2 + 1 =(m − j + 1)(m − j)2умножении при фиксированном j.Вдобавок к этому (m − j) делений + 1 деление от (*), откуда получаем, что прификсированном j всего(m − j + 1)(m − j + 2)⇒2умножений и делений.Теперь, отпустив j, получимmX(m − j + 1)(m − j + 2)(5)2j=1делений и умножений при решении системы вида B · Y = f .Задача. Доказать, что для реализации (5) необходимоm(m + 1)(m + 2)операций.6Решение:mX(m − j + 1)(m − j + 2)j=12= { проведем замену k = m − j + 1} =mXk(k + 1)k=1=mXkk=12+2=mXk2k=12k(k + 1)k(k + 1)(2k + 1), а вторая —.
Сло412m(m + 1)(m + 2)жив полученные результаты, получим требуемое значение.6Откуда следует, что первая сумма равна11m(m − 1)действий (обратный ход метода Гаус2m2 (m − 1)са), а учитывая, что нам надо решить m уравнений, получимумножений2и делений.Таким образом, весь процесс Гаусса-Жордана обращения матрицы требует m3действий (факторизация + решение системы (3) + решение системы (4)):Для реализации (4) необходимоm3 − m m(m + 1)(m + 2) m3 − m2++= m3 .362§4 Метод квадратного корняРассмотрим матричное уравнение видаAx = f,(1)где A — матрица размера (m × m), |A| =6 0, A = A∗ (эрмитова, т.е. aij = aji ).Представим матрицу A в виде A = S ∗ DS, гдеd11 0 .
. . 0 0 0 d22 . . . 0 0 D = ...... .... , ... .. 00 . . . 0 dmms11 s21 . . . 0 s22 . . .S = ...... ...00 ...s1ms2m .. ,. smmгде sii > 0, dii = ±1, i = 1, m.Следовательно S ∗ будет нижнетреугольная.Покажем для симметричной матрицы второго порядка, что такое разложениевозможно. Для простоты потребуем, чтобы она была вещественной.a11 a12s11 s12d11 0s11 0∗TA=, S=, D=, S =S =a21 a220 s220 d22s12 s22Найдем в начале произведение DS. d11 0s11 s12d11 s11 d11 s12DS =·=0 d220 s220d22 s22Теперь домножим на S ∗ слева. s11 0d11 s11 d11 s12d11 s211s11 d11 s12∗S DS =·==As12 s220d22 s22s11 d11 s12 d22 s222 + s212 d1112Откуда получим, чтоa11a12a21a22= d11 s211= s11 d11 s12= a12= d22 s222 + s212 d11Для того чтобы факторизация была возможна, необходимо, чтобы система быларазрешима. Ее решение находится по формулам:d11 = sign(a11 )p|a11 |s=11a12s12 =s11 d11 d22 = sign(a22 − s212 d11 )q s22 = |a22 − s2 d11 |12Таким образом показано, что для матрицы A (вещественной и симметричной)возможна факторизация и все находится по явным формулам.
Рассмотрим общийслучай. Элементы матрицы DS находятся по формуле:(DS)ij =mXnodil slj = с учетом свойств матриц = dii sijl=1Заметим, что(S ∗ )ij = S jiУмножим слева матрицу S ∗ на DS:∗(S DS)ij =mXsli dll sljl=1Распишем сумму, приравняв ее к aij элементуi−1XmXsli dll slj + sii dii sij +l=1sli dll slj = aij ,i, j = 1, ml=i+1Учтем специфику матрицы S: sli = 0,mXl > i. Следовательно,sli dll slj = 0l=i+1Тогда получимsii dii sij = aij −i−1Xl=1sli dll slj ,i6j(3)13Рассмотрим два случая (i = j и i < j).Положим i = j. Получим:sii sii dii = aii −i−1Xsll sli dlll=1Так как zz = |z|2 , то2|sii | dii = aii −i−1X|sli |2 dlll=1Таким образом, мы получили формулы, из которых можем найти диагональныеэлементы матрицы D:i−1X|sli |2 dll )dii = sign(aii −l=1Далее понятно, чтоvui−1Xutaii −|sli |2 dll sii =l=1Осталось вновь вернуться к исходной формуле, частный случай которой мырассмотрели.
Из формулы (3) можно найти элементы sij , считая, что dii sii 6= 0:sij =aij −Pi−1l=1sli dll sljsii dii,i < j.Система перепишется в виде:S ∗ DSx = fОбозначая DSx = Y , получим:(S ∗Y = fDSx = Y(4)(5)Матрица S ∗ нижнетреугольная и из (4) легко находится вектор Y . Затем, решаяуравнение (5), находим вектор x.Был рассмотрен метод квадратного корня. Этот метод имеет преимущество поколичеству арифметических действий (умножений и делений) перед методом Гаусса,m3m3который пропорционалендействий, а рассмотренный нами способ , так как36расчеты проводились только для части матрицы (но здесь есть еще m извлеченийкорня).Следовательно, если матрица эрмитова или самосопряженная, то один из наиболее эффективных прямых методов - метод квадратного корня.14§5 Примеры и канонический вид итерационных методов решения СЛАУРассмотрим матричное уравнение видаAx = f,(1)где A — матрица размера (m × m), |A| =6 0,x = (x1 , x2 , .
. . , xm )T ,f = (f1 , f2 , . . . , fm )T .Когда решается линейная система (1), то в правой части, как правило, стоитфункция наблюдения. Ясно, что прямой метод дает точное (числовое) решение, абстрагируясь от округления - в силу конечной разрядной сетки.Также, если m достаточно велико, то количество действий даже для суперкомпьютера может оказаться очень большим. А использование итерационных методовпозволит решить систему (1) быстрее. Более того, в некоторых случаях возможнооценить число итераций, необходимых для достижения заданной точности.В качестве примера рассмотрим методы Якоби и Зейделя. Запишем систему (1)покоординатно:mX(2)aij xj = fi , i = 1, mj=1Рассмотрим метод Якоби (МЯ).
Перепишем равенство (2) следующим образом:i−1Xaij xj + aii xi +j=1mXaij xj = fi ,i = 1, m.j=i+1Тогда, отсюда можно выразить xi :fi −xi =i−1PmPaij xj −j=1aij xjj=i+1aii,c предположением, что aii 6= 0.Обозначим xni — n-я итерация i-й координаты.Чтобы организовать метод Якоби, слева «навешивают» итерацию n+1, а справасобирают все с n-ой итерацией.fi −xn+1=ii−1Paij xnj −j=1mPj=i+1aiiaij xnj,(3)где n = 0, 1, 2, . .
. , x0 — задано.В каждом конкретном итерационном методе выбор начального приближенияявляется самостоятельной задачей.15Понятно, что никакой процесс не может продолжаться бесконечно долго, нужнобудет где-то оборвать это действие. Заканчиваем итерационный процесс тогда, когдадостигается оценка:kxn − xk < εЗапишем метод Зейделя:fi −i−1Paij xn+1−jj=1xn+1=imPaij xnjj=i+1aii,(4)где n = 0, 1, 2, . . . , x0 - задано.Формулу (4) можно так же реализовать по явным методам.
Рассмотрим этотспособ.Положим i = 1. Найдем первую координату (n + 1)-ой итерации:f1 −xn+1=1mPa1j xnjj=2a11Найдем вторую координату:f2 − a21 xn+1−1xn+1=2mPa2j xnjj=3a22И так далее. То есть, организуя вычислительный процесс с первой координаты,можно найти все значения вектора x по явным формулам для (n + 1)-й итерации.Таким образом, метод Зеделя реализуется по явным формулам. При изученииитерационных методов важно не упускать из виду два вопроса: первый - сходимостьметода и условия этой сходимости, второй - получение оценки скорости сходимостиметода.Для исследования этих вопросов удобно рассматривать системы в матричномвиде. Любую матрицу мы можем представить в виде суммы трех матриц R1 +D+R2 ,где00..R1 = .
— нижнетреугольная матрица c нулевой диагональю,aij0a110..D= — диагональная матрица,.0amm0aijR2 = . . . — верхнетреугольная матрица с нулевой диагональю.0016Тогда перепишем уравнение (1):R1 x + Dx + R2 x = f⇓Dx = f − R1 x − R2 xЕсли предположить, что матрица D имеет обратную, а это как раз и требуетвыполнения условия aii 6= 0, то тогда, домножив слева на матрицу D−1 , получаемx = D−1 f − D−1 R1 x − D−1 R2 x.Итак, организуем итерационный процесс для метода Якоби и Зейделя.МЯ:xn+1 = D−1 f − D−1 R1 xn − D−1 R2 xn(5)МЗ:xn+1 = D−1 f − D−1 R1 xn+1 − D−1 R2 xn(6)Dxn+1 = f − R1 xn − R2 xn(7)(D + R1 )xn+1 = f − R2 xn(8)D(xn+1 − xn ) + Axn = f(9)(D + R1 )(xn+1 − xn ) + Axn = f(10)МЯ:МЗ:МЯ:МЗ:где n = 0, 1, 2, .
. . , x0 — задано.Видно, что если есть сходимость, то сходится к решению нашей системы. Такматематики и пришли к каноническому виду двухслойного итерационного процесса.Определение. Канонической формой записи двухслойного итерационного методарешения системы (1) называется запись видаBn+1xn+1 − xn+ Axn = f,τn+1(11)−1где n = 0, 1, 2, . . . , x0 — задано, τn+1 > 0 (итерационный параметр) и ∃Bn+1. ЕслиBn+1 = E, то метод (11) называется явным, в противном случае – неявным.Замечание. Если параметр τ и матрица B зависят от итерации, то метод называется нестационарным, иначе стационарным.
Далее мы будем рассматриватьтолько стационарные методы.17Отметим, что если Bn+1 – диагональная матрица, то метод реализуется по явнымформулам, хотя формально метод неявный.Рассмотрим еще один метод — метод простой итерации (метод релаксации).xn+1 − xn+ Axn = f,ττ >0n = 0, 1, 2, .
. . , x0 – задано.Если τ зависит от n, то получающийся методxn+1 − xn+ Axn = f,τn+1называется методом Ричардсона.Попеременно треугольный итерационный метод (метод Самарского)Представим матрицу A в виде A = R1 + R2 , где0.5a110..R1 = — нижнетреугольная матрица,.aij0.5amm0.5a11aij..R2 = — верхнетреугольная матрица..00.5ammТогда попеременно треугольный итерационный метод имеет следующий вид:(E + wR1 )(E + wR2 )xn+1 − xn+ Axn = f,τ(12)где n = 0, 1, .
. . , x0 — задано, τ > 0, w > 0 – итерационные переметры.В данном методе имеется два итерационных параметра. С точки зрения алгоритма реализации, этот метод не сложнее, чем предыдущие, а по сходимости - напорядок лучше.Введем векторxn+1 − xn(E + wR2 )= W n+1τОпределение.
Разность между правой и левой частями решаемой системы называется невязкой.В нашем случае невязка имеет видf − Axn = rn18Тогда видно, что на первом этапе можем решать уравнение(E + wR1 )W n+1 = rnПравая часть известна, а (E + wR1 ) — нижнетреугольная матрица. Нахождение вектора системы с нижнетреугольной матрицей осуществляется по явным формулам,начиная с первой компоненты вектора W .Теперь, на втором этапе решаем уравнение(E + wR2 )V n+1 = W n+1 ,гдеxn+1 − xn.V=τНа третьем этапе (n + 1)-ю итерацию находим по формуле:n+1xn + τ V n+1 = xn+1 .Таким образом, несмотря на то, что метод Самарского неявный с двумя итерационными параметрами, его реализация не представляет никакой трудности.§6 Теоремы о сходимости итерационных методовРассмотрим матричное уравнение видаAx = f,(1)где A — матрица размера (m × m), |A| =6 0.Рассмотрим также двухслойный стационарный метод решения уравнения (1):Bxn+1 − xn+ Axn = f,τ(2)где τ > 0, существует обратная матрица B −1 , n = 0, 1, 2, .
. . , и задано начальноеусловие x0 .Когда мы говорим, что начальное условие задано, то это значит, что либо егонадо выбирать исходя из каких-то жестких условий, либо есть свобода выбора.Говоря о сходимости нужно четко понимать, по какой норме эта сходимостьполучается, поэтому нам необходимо ввести линейное нормированное пространство.Пуcть H - линейное пространство размерности m:dim H = mВозьмем два произвольных вектора x и y из этого пространства:x ∈ H,x = (x1 , x2 , . . . , xm );y ∈ H,y = (y1 , y2 , . . .