Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности (1051191), страница 3
Текст из файла (страница 3)
Так как j < k, то все коэффициенты матрицы U, присутствующие в (2.14),уже определены.8Теперь для i > k элементы нижнетреугольной матрицы L равны нулю.16П РИМЕР. Для матрицы (1.2) приведенный алгоритм дает следующие матрицы L и U (с точностью три знака):L =U =10100.09110.222 0.091 0.18510.111000.08510000010.222 0.182000.235 0 1901103219.818 1.9097.8891000.77811.823;(2.17)0102 00 00 .0 0.889 80 7.205(2.18)Соответствующая им матрица ошибки разложения R равнаR = A − LU = 00000000000000000000000000 −0.364 −0.848000000000000 −0.182 0 −0.404 .000000Замечание.
Элементы матриц L и U, полученных в результатеILU-разложения, не совпадают с соответствующими элементами матриц LA и UA , полученных при полной LU-факторизации. В этомлегко убедиться, сравнивая (2.8) и (2.9) с (2.17) и (2.18).Очевидно, что матрица LU удовлетворяет всем трем требованиям к матрице предобусловливателя, сформулированным в §2.1. Действительно, она приближает матрицу A, так как на множестве индексов PA точно воспроизводит ее; она легко вычисляется по приведенному выше несложному алогоритму; наконец, она легко обратима, так как является произведением двух треугольных матриц.Таким образом, выбор M = LU является достаточно хорошимспособом предобусловливания.2.3.Симметричный случайВ случае, когда матрица A симметрична, из равенств A = LU и A == AT для треугольных матриц L и U следует, что L = U T .
Фактори17зация (2.10) в этом случае принимает видA = LLT + R = UTU + R .(2.19)Очевидно, однако, что для нахождения диагональных элементовтакой матрицы L требуются операции извлечения квадратного корня. Чтобы избежать этого, вместо факторизации (2.19) используетсянесколько отличающийся вариантA = LDLT + R ,(2.20)где главная диагональ L состоит из единичных элементов, а D —диагональная матрица.Для нахождения коэффициентов L и D применяется такой жеподход, как и в §2.2, поэтому сразу приведем алгоритм без пояснений:Для k = 1, nДля j = 1, k − 1Если (k, j) ∈/ PAто lkj := 0иначе lkjувеличить jdk := akk −увеличить kk−1Pi=1j−1X1 :=akj −di lik ljk dji=12 d2lkiiП РИМЕР.
Для симметричной матрицыA=90300 3 08 0 10 11 11 1 9приведенный алгоритм дает (с точностью три знака):L=1010.3330100.125 0.091 1;D = diag (9, 8, 11, 8.784) .Соответствующая им матрица ошибки разложения R = A − LDLTравна нулю.182.4.О программной реализацииILU-предобусловливанияОдним из условий при построении ILU-разложения было соответствие портретов PL ⊂ PA и PU ⊂ PA .
При этом одна из матрицявляется верхне-, а другая нижнетреугольной, так что их портретыне совпадают, за исключением диагональных элементов. Однако алгоритм разложения построен так, что l kk ≡ 1 и, следовательно, явнохранить диагональные элементы L необходимости нет.Таким образом, PL+U ⊂ PA и разложение лучше всего хранитькак еще одну матрицу с тем же портретом что и A — вследствиесовпадения портретных массивов они хранятся только один раз, ак массиву элементов aelem (для CSlR — altr, autr, adiag) добавляетсямассив luelem (для CSlR — массивы ltr, utr, ludiag).
Это эквивалентнохранению объединенной матрицы разложенияB = L + U − I,диагональ которой считается отнесенной к матрице U.Для нахождения ILU-факторизации формат хранения CSlR оказывается более удобным по сравнению с CSR. Как можно видеть из §2.2,в алгоритме для вычисления коэффициентов l kj и ukj используютсяскалярные произведения «строка-столбец» видаmax{i,k}−1Xlki uij ,(2.21)i=1в которых доступ к элементам нижнетреугольной матрицы L производится по строкам, а к элементам верхнетреугольной матриицыU — по столбцам. Для CSlR порядок доступа совпадает с порядкомхранения для обеих матриц, тогда как для CSR — только для матрицы L.
Для постолбцового доступа к элементам U необходимы операции поиска, что существенно замедляет расчеты 9 .При использовании разреженных строчных форматов внешниециклы по перебору индексов элементов матрицы в алгоритме из§2.2, очевидно, должны быть переписаны. Построчный доступ кэлементам матрицы L выполняется непосредственно, для построчного доступа к элементам U необходимы дополнительные операции9С другой стороны, как будет показано ниже, из-за построчной генерации матриц операции поиска необходимы и при использовании CSlR, однако в этом случаена каждую сумму вида (2.21) нужен только один поиск, тогда как для CSR операция поиска должна производиться для каждого слагаемого суммы. Можно показать, что вычислительная сложность для CSR на порядок выше, чем для CSlR.19поиска.
Общая схема нахождения ILU-факторизации выглядит следующим образом10 :для k от 1 до n {для j от iptr[k] до iptr[k+1]-1 {вычислить ltr[j] по формуле для lk,jptr[j]}вычислить ludiag[k] по формуле для u k,kдля j от k+1 до n {найти такое q (iptr[j]<=q<iptr[j+1] что jptr[q]=kесли найденото вычислить utr[q] по формуле для u k,j}}Для решения треугольных СЛАУ с матрицами L и U (см. формулу (2.4)) могут быть использованы алгоритмы из §1.4.10Снова (ср. §1.3) упорядоченность элементов внутри строки способна дать некоторый выигрыш при поиске.20Глава 3Классические итерационныеметоды и релаксацияИсторически первые итерационные методы основывались на циклическом покомпонентном изменении вектора решения, осуществляемом таким образом, чтобы обнулить соответствующий коэффициент вектора невязки и тем самым уменьшить его норму.
Подобнаяметодика уточнения решения получила название релаксации.Хотя в настоящее время такие методы в их классической формулировке уже практически не применяются, существуют определенные классы задач [4], для которых разработаны их модификации,хорошо себя зарекомендовавшие.
Кроме того, как будет показано в§3.3, эти методы могут быть применены (и применяются) не в качестве самостоятельного средства решения СЛАУ, а для предобусловливания1 .3.1.Методы Якоби и Гаусса-ЗейделяПусть матрица A системы (1) такова, что ее главная диагональ несодержит нулевых элементов. Представим ее в виде разности трехматриц:A = D − E − F,(3.1)где матрица D содержит только диагональные элементы A; матрицаE — только поддиагональные; матрица F — только наддиагональные (см. рис. 3.1)Система (1) тогда может быть записана в видеDx − Ex − Fx = b .1См. также §4.4, где идея методов, основанных на расщеплении, используетсядля построения подпространств Крылова.21Рис. 3.1.
Представление матрицы в виде разностиЕсли имеется некоторое приближение x k к точному решениюСЛАУ x∗ , то при xk 6= x∗ это соотношение не выполняется. Однако,если в выраженииDxk − Exk − Fxk = b(3.2)одно или два из вхождений вектора xk заменить на xk+1 и потребовать, чтобы равенство имело место, можно получить некоторую вычислительную схему для уточнения решения.Наиболее простой с точки зрения объема вычислительной работы вариант получается при замене в (3.2) Dx k на Dxk+1 . При этомполучается схемаxk+1 = D−1 (E + F)xk + D−1b ,(3.3)известная как метод Якоби2 .Выражение (3.3) в скалярной форме имеет вид(k+1)xii−1nXX1 (k)(k)=bi −aij xj −aij xj ,aiij=1j=i+1i = 1, n ,(3.4)откуда хорошо видна основная идея метода: на (k + 1)-й итерацииi-й компонент вектора решения изменяется по сравнению с k-й итерацией так, чтобы i-й компонент вектора невязки r k+1 стал нулевым(при условии отсутствия изменений в других компонентах вектораx).Очевидным недостатком схемы (3.3)–(3.4) является то, что при(k+1)нахождении компонента xiникак не используется информация(k+1)(k+1)о уже пересчитанных компонентах x1, .
. . , xi−1 . Исправить этотнедостаток можно, переписав (3.4) в виде(k+1)xi2=i−1XnX1 (k+1)(k)bi −aij xj−aij xj aiij=1j=i+1i = 1, n ,(3.5)Такая форма записи подразумевает, что главная диагональ матрицы не содержит нулевых элементов. Если это не так, то при невырожденности матрицы выполнения условия можно добиться перестановкой строк.22что в векторной форме эквивалентноxk+1 = (D − E)−1 Fxk + (D − E)−1b .(3.6)Такая схема называется методом Гаусса-Зейделя.Выражение (3.5) получается из (3.2) заменой x k на xk+1 при матрицах D и E. Если вместо этой пары матриц взять D и F, то получится похожая схемаxk+1 = (D − F)−1 Exk + (D − F)−1b ,(3.7)которая называется обратным методом Гаусса-Зейделя.Еще одной очевидной модификацией является симметричныйметод Гаусса-Зейделя, который заключается в циклическом чередовании (3.6) и (3.7) на соседних итерациях.Нетрудно заметить, что выражения (3.3), (3.6), (3.7) могут бытьзаписаны в видеKxk+1 = Rxk + b ,(3.8)где матрицы K и R связаны соотношениемA = K − R.(3.9)Подобное представление матрицы A называется расщеплением, аметоды вида (3.8) — методами, основанными на расщеплении.
Очевидно, матрица K должна быть невырожденной и легко обратимой.Условия сходимости методов (3.8) устанавливает следующая теорема:Теорема 3.1.1 Вычислительная схема (3.8) сходится при любомначальном приближении x0 тогда и только тогда, когда матрицыK и R удовлетворяют условию maxj |λj (K−1 R)| < 1.Доказательство. В силу соотношения (3.9) для точного решения x ∗системы (1) справедливоKx∗ = Rx∗ + b .Почленно вычитая из этого равенства (3.8), получимK(x∗ − xk+1 ) = R(x∗ − xk ) ,(3.10)где векторы в левой и правой частях есть ошибки решения ϑ k+1 и ϑkсоответственно.Очевидно, сходимость схемы (3.8) к точному решению эквивалентна сходимости ϑk → 0 при k → ∞. Из (3.10) имеемϑk+1 = K−1 Rϑk = (K−1 R)k+1 ϑ0 .23(3.11)Следовательно, для сходимости (3.8) необходимо и достаточновыполнения условия (K−1 R)k → 0 при k → ∞, откуда и вытекаетутверждение теоремы.2Замечание 1. Из (3.11) видно, что чем меньше присутствующая вусловии величина maxj |λj (K−1 R)|, тем быстрее сходимость метода.Замечание 2.
Матрицы K и R в (3.8) могут, вообще говоря, зависеть от номера итерации (K = Kk , R = Rk )3 . В этом случае условияmaxj |λj (K−1k Rk )| < 1 достаточно для гарантированной сходимости,однако оно перестает быть необходимым.3.2.Ускорение сходимостирелаксационных методовИз теоремы 3.1.1 следует, что скорость сходимости методов, основанных на расщеплении, непосредственно связана со спектральнымрадиусом матрицы K−1 R; с другой стороны, выбор K ограничен требованием легкой обратимости.Одним из распространенных способов улучшения сходимости является введение параметра.
Пусть ω — некоторое вещественное число. Рассмотрим вместо (1) масштабированную системуωAx = ωb ,(3.12)и вместо (3.1) воспользуемся представлениемωA = (D − ωE) − (ωF + (1 − ω)D) ,(3.13)где матрицы D, E и F имеют тот же смысл, что и в (3.1).Тогда на основании (3.12) и (3.13) можно построить итерационную схему, похожую на метод Гаусса-Зейделя,(D − ωE)xk+1 = [ωF + (1 − ω)D]xk + ωb .(3.14)Схема (3.14) называется методом последовательной верхней релаксации (SOR4 ) Для нееKSOR (ω) = D − ωE ;RSOR (ω) = ωF + (1 − ω)D .Выбор параметра ω, минимизирующего спектральный радиусρ(ω) = maxj |λj (K−1SOR (ω)RSOR (ω))| является, вообще говоря, достаточно сложной проблемой.
Однако для многих классов матриц (например, связанных с конечно-разностными аппроксимациями уравнений в частных производных [4, 5]) такая задача исследована иоптимальные значения ω известны.34Такие итерационные схемы называются нестационарными.Successive OverRelaxation.24Замечание. Записав (3.14) в скалярной форме, нетрудно получить, что в SOR решение удовлетворяет соотношению(k+1)xi(k)= ωxGSi + (1 − ω)xii = 1, n ,где xGSi — итерация метода Гаусса-Зейделя, определяемая формулой(3.5). Действительно, при ω = 1 (3.13) совпадает с (3.1).Выражение (3.13) остается тождеством, если в нем поменять местами матрицы E и F.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
