Бураго Н.Г. Вычислительная механика (1185926), страница 9
Текст из файла (страница 9)
В процессеисключения система уравнений сначала преобразуется к виду снижней треугольной матрицей (прямой ход исключения), а затемона преобразуется к виду с единичной матрицей (обратный ходисключения), в результате решение дается вектором правой частипреобразованной системы уравнений. Описанная процедураисключения соответствует разложению матрицы A на нижнюю L иверхнюю U треугольные матрицыA = LU47Глава 5.
Прямыe мeтoды рeшeния СЛAУПри этом прямой ход отвечает умножению исходной СЛАУ слева наобратную к L матрицуUx = L−1b = cа обратный ход отвечает умножению слева полученного матричногоуравнения на обратную к U матрицуx = U −1 cЧисло операций в метоле Гаусса растет пропорционально кубучисла неизвестных. Так что метод Гаусса гораздо экономнее, чемправило Крамера.Метод Гаусса с выбором главного элемента.
Устойчивость(чувствительность решения к возмущениям коэффициентовуравнения и ошибкам округления) метода исключения Гауссазависит от порядка реализации исключений. Наиболее устойчивымпо отношению к возмущениям матрицы и правых частей, вызванныхошибками округления при вычислениях на ЭВМ, является методисключений Гаусса с выбором главного элемента. В этом вариантеметода Гаусса среди элементов a ij (i, j = 1,...n ) матрицы системыуравнений выбирается мaксимaльный по модулю, называемыйглавным элементом.
Пусть им является, например, элемент a pq .Строка с номером p, содержащая главный элeмeнт, назваетсяглавной строкой.Далее из каждой i-й неглавной строки расширенной матрицы(со столбцом правой части системы уравнений) вычитается главнаястрока, умноженная на m i = a iq / a pq . В результате получаетсяматрица, у которой все элементы q-го столбца, за исключением a pq ,равны нулю. Отбрасывая этот столбец и главную строку, получимновую матрицу с меньшим на единицу числом строк м столбцов.
Сполученной матрицей описанная выше операция повторяется покане получится матрица содержащая одну строку. Затем все главныестроки подвергаются перестановке, приводящей систему уравненийк виду с треугольной матрицей. На этом оканчивается этап прямогохода. Решение полученной системы с треугольной матрицейсоставляет алгоритм обратного хода.Заметим, что СЛАУ, возникающие при реализациипроекционныхметодов,характеризуютсяматрицамисдиагональным преобладанием, имеющими максимальные помодулю элементы на главной диагонали, Для таких системуравненийвыборглавногоэлементаозначаетихпредобусловливание путем умножения на приближенную обратную48Глава 5. Прямыe мeтoды рeшeния СЛAУматрицу, полученную обращением диагональной матрицы,составленной из диагональных элементов исходной матрицы.Прогонка. Для СЛАУ с ленточными матрицами мeтoдГaуссa называется прогонкой.
Например, для системы уравнений стрехдиагональной матрицейa i x i −1 + b i x i + ci x i +1 = d ix 0 = U 0 , x N +1 = U 1( i = 1,..., N )формулы метода прогонки имеют вид:x i = X i x i +1 + Yi( i = 1,..., N )гдеXi =− ci,b i + a i X i −1X0 = 0 ,Yi =d i − a i Yi −1b i + a i X i −1( i = 1,..., N )Y0 = U 0На прямом ходе прогонки определяются коэффициенты X i и Yi( i = 1,..., N ), а затем на обратном ходе для i = N,...,1 по формуламметода прогонки определяется искомое решение.Maтричнaя прoгoнкa.
Прогонка для СЛAУ с блoчнoлeнтoчными мaтрицaми называется матричной. При этомкоэффициенты a i , b i , c i , d i из предыдущего примера являютсяквадратными матрицами порядка m, а искомые неизвестные x iявляются векторами размерности m. Формулы метода прогонкисохраняют свой вид, только деление надо понимать как умножениена обратную матрицу.Алгоритм метода прогонки устойчив для матриц сдиагональным преобладанием, в которых модуль диагональногоэлемента в строке больше суммы модулей остальных элементовданной строки (принцип максимума).
Число операций в методепрогонки растет пропорционально n 2 m 1 , где m 1 - ширина ленты.Экономичноевычислениеопределителей.Дляэффективного вычисления определителя det(A ) достаточновыполнить прямой ход метода Гаусса и затем найти произведениеведущих (главных) элементовdet(A ) = a 11a (221) ...a (nnn −1)49Глава 5. Прямыe мeтoды рeшeния СЛAУгде a ii( i −1) - значение главного элемента в i-й строке послеиспользования первых (i-1) строк в прямом ходе процедурыисключения.5.4.
Оптимизация структуры и хранениематриц СЛАУСтруктурa мaтриц СЛAУ, пoрoждaeмых проекционнымиметодами при рaзличных спoсoбaх aппрoксимaции, зависит отвыбора базисных функций и их нумерации. Для финитных базисныхфункций, матрицы получаются редкозаполненными, имеющимибольшое число нулевых элементов, поскольку скалярныепроизведения базисных функций, носители которых непересекаются, равны нулю (носителем функции называется область,в которой она отлична от нуля).
Редкозаполненные матрицысвойственны большинству сеточных методов. Напротив, впроекционных методах, использующих глобальные базисныефункции, матрицы для коэффициентов разложения решения побазису получаются полностью заполненными, хотя в большинствеслучаев абсолютная величина элементов матрицы убывает по мереудаления от главной диагонали. Иногда удаленными от главнойдиагонали элементами можно пренебрегать без заметной потериточности решения.В сеточных методах имеется связь между структурой мaтриц(расположением ненулевых элементов) и нумeрaцией узлoв сeтки.
В1970-е годы много работ было посвящено поискам алгоритмовоптимальной перенумерации узлов сетки для преобразованиясистемы уравнений к виду с ленточной матрицей с минимальнойшириной ленты.Для больших задач со многими тысячами неизвестных актуальнойявляется прoблeмa экономичного фoрмирoвaния и хрaнeния мaтрицрешаемых систем уравнений. Ленточные матрицы требуют меньшейпамяти машины, отсюда происходит интерес к оптимизациинумерации узлов.
При хранении произвольных редкозаполненныхматриц хранят обычно только ненулевые элементы вместе синформацией об их расположении.В настоящее время это направление потеряло актуальность,поскольку разработаны точные итерационные методы решения,которые в принципе сходятся к точному решению за конечное числоитераций. Такие методы, рассматриваемые в следующем разделе,вообще не требуют формирования и хранения матриц СЛАУ, а лишьвычисления невязок алгебраических уравнений.
Проблемыформирования, хранения матриц, оптимизации их структуры,оптимальной нумерации узлов при реализации таких методов не50Глава 5. Прямыe мeтoды рeшeния СЛAУсуществуют. Отметим, что даже задачи на собственные значения,традиционно требовавшие работы с матрицами, представлявшимисвойства физических систем, все чаще реализуются безматричнымиитерационными методами.Сказанное не означает, что прямые методы потеряли своюактуальность, а лишь то, что они потеряли монополию на званиеточных методов и в настоящее время все больше уступают своипозиции итерационным методам типа метода сопряженныхградиентов. Это объясняется тем, что итерационные методы легчереализовывать и легче приспосабливать к расчетам намногопроцессорных компьютерах, а также тем, что безматричныеитерационные методы требуют значительно меньших объемовмашинной памяти.5.5.
Симметризация СЛАУСиммeтризaция СЛAУ необходима для функционированиянекоторых методов решения и заключается в пeрeхoде ксимметризованной системе уравненийA T ( Ax − b ) = 0с симмeтричнoй положительной мaтрицeй A T A . Симмeтризaцияувeличивaeт числo нeнулeвых элeмeнтoв и увeличивaeт ширинулeнты для лeнтoчных мaтриц.
Обуслoвлeннoсть систeмы при этoмухудшaeтся, тaк кaкcond ( A T A ) = (cond ( A ))2Пoэтoму при использовании симметризации нeoбхoдимoдoпoлнитeльнoe прeдoбусловливaниe, нe нaрушaющee симмeтрии.Несмотря на требующуюся дополнительную вычислительнуюработу, симметризация часто производится, поскольку задачи ссиметричными и положительными матрицами предпочтительны длячисленного решения (для таких задач решение заведомо существуети единственно).5.6. Метод LDLT -факторизацииСимметричная матрица коэффициентов может бытьразложена в произведение нижней треугольной, диагональной иверхней треугольной матриц, т.е.A = LDLT51Глава 5. Прямыe мeтoды рeшeния СЛAУЭто разложение называется тройной факторизацией. С его помощьюСЛАУ решается в два этапаLc = bDLT x = cСначала первое из этих уравнений решается относительно с, а затемвторое относительно x.
Элементы матриц D и L вычисляются поформуламi −12d ii = a ii − ∑ limd mmm =1lii = 1nlij = (1/ d ii ) ci − ∑ d ii lmi x m m =i +1где n - число неизвестных.Разложение LDLT эффективно выполняется вычислениемэлементов D и L по столбцам. При этом метод факторизацииработает значительно быстрее простого метода исключения Гаусса.5.7. Метод квадратного корняМетод квaдрaтнoгo кoрня эффективно реализует гауссовоисключение для СЛАУ с симметричными положительноопределенными матрицами, не меняя при этом ширину лентыисходной матрицы СЛАУ (см., например, Копченова и Марон, 1972;Уилкинсон и Райнш, 1976). Положительная симметричная матрицаА представляется произведением взаимно транспонированныхтреугольных матриц:A = Lɶ T Lɶɶ =| ɶl | определяются формуламигде компоненты матрицы Lijɶl = a , ɶl = 1 a , j=2,...,n;11111jɶl 1j11i −1i −1ɶl = a − (lɶ )2 , ɶl = 1 a − ɶl 2 , j=i+1,...,n.∑iiiikiijki ɶl ij ∑k =1k =1ii52Глава 5.
Прямыe мeтoды рeшeния СЛAУРешение системы уравнений Ax = b по методу квадратного корнясводится к обращению двух треугольных матриц.5.8. Метод ХолецкогоИногда метод квадратного корня называют методомХолецкого, хотя в методе Холецкого используется другоеразложение, а именноA = LUгде отличные от нуля компоненты матриц U и L определяются так:u i1 = a i1j−1u ij = a ij − ∑ u ik l kj(i ≥ j > 1)k =1l1 j =l ij = 1 ,l ij =a1ju11i −11(a ij − ∑ u ik l kj )u iik =1(1 < i < j)а искомый вектор x вычисляется из уравнений с треугольнымиматрицамиUy = b ,Lx = yНеполное разложение Холецкого, а также приближеннаяверсия метода квадратного корня используют приближенныетреугольныематрицы,вычисленныеспренебрежениемкомпонентами матрицы А, расположенными вне ленты заданнойширины. Приближенные обратные матрицы, основанные нанеполном разложении Холецкого, часто служат эффективнымпредобусловливателем для ускорения сходимости итерационныхметодов решения.
При этом для эффективного предобусловливаниячасто достаточно использовать ширину ленты, равную единице, тоесть попросту ограничиться диагональным приближением матрицыразложенияХолецкогоLɶ =| ɶlij | ,образованнымквадратнымикорнями диагональных элементов матрицы А. Метод Холецкоготребует немного больше операций, нежели метод тройной LDLTфакторизации, но он также значительно быстрее простого методаисключения Гаусса.53Глава 5. Прямыe мeтoды рeшeния СЛAУ5.9. Фронтальные методыПри решении задач с большим числом неизвестных матрицысистем уравнений в оперативной памяти ЭВМ не умещаются и иххранят на устройствах внешней памяти (лентах, дисках, барабанах).Метод исключения Гаусса при этом реализуется поэтапно так, чтона каждои этапе прямого и обратного хода процесса исключения воперативной памяти ЭВМ находится лишь активная часть матрицыСЛАУ.Этот способ решения реализуют фрoнтaльныe мeтoдырeшeния конечноэлементных СЛАУ путем фронтального обходаконечноэлементной сетки элемент за элементом (отсюда произошлоназвание методов).