Fletcher-1-rus (1185917), страница 42
Текст из файла (страница 42)
Мга(и+1! о. 1Г(2*ИВ .ЕО. М)В(2,М+1) О. 50 б 1 1,2 ВО 5 Э 1,ИН ЭЛ = 2"3 ООАК 1,1 ЗВ Л "1+К 4 А(ЭВ) В(ЭВ! "В(1,ЭВ)'В(ЭВ-1) 5 СОЯТ1ИОЕ б СОИТ1ЯОЕ ИКИ ИМ вЂ” 1ВС Х(И) В(М)/В(З,И) 1РПВС .ВЭ. ЫХ(И-и = СА(и-и - В(б,н-п*Х(И))/В(Э,ИОО 7 Э 1.МЕМ ЗА И-2*4+1-1ВС Х(ЗА) (Е(ЭА) - В(4,ЭЫ*Х(Эх+1) - В(5/ОА)*Х(Л+2))/В( Х(Л-1) = (В(Л-Ы - В(4,Л-1)*Х(Л))/В(Э,Л-1) 7 СОЯТ1ИОЕ ВЕТОВИ ЕЯО Рнс. 6.)9. Распечатка подпрограммы ВАР)6О).. $6.2.
Прямые методы для линейных систем 243 ния Ч заполнить левые части уравнений (6.23) и убедиться в равенстве левых частей правым. Это, однако, не гарантирует точности построенного решения Ч, если матрица А плохо обусловленная, иначе говоря, если уравнения, составляющие систему (6.23), близки к линейным. Плохо обусловленные системы и их влияние на точность решения обсуждаются в книгах [Оега1д, 1987; Ра!т!с(ц(з1, В)огск, 1974].
б.2.4. Обобщенный алгоритм Томаса Ьс с, Ье се а, Ь, с, е, а, ее сс сс ан , Ьн , ен с сн с ен ан Первый этап состоит в исключении всех элементов еь что дает сс с, а' Ь' с' а' Ь' с' с с с с е е а,, Ьн, с'„, а'„Ь', Применение конечно-разностных или конечно-элементных схем повышенного порядка приводит к матрицам с более широкими лентами, чем трехдиагональная матрица.
Однако не. трудно провести обобщение алгоритма Томаса. Допустим, что нам необходимо решить следующую пятидиагональную систему: Гл. 6. Стациоиариые задачи где Ь, '= Ь, — е,с',,/а',, се = сс еД ~/ц еД ~~ц (6.32) Следующий этап по существу совпадает с первым этапом обыч- ного алгоритма Томаса. Иначе говоря, все элементы а', исклю- чаются, что дает 1с",), 1 с" ~' 1 с" ,)", (а с (6.33) а 1 с,, где I е а с~ — аД а с,= е а Ь~ — асск-1 (6.34) е е а Ьз — а~с~ е а лс аА-~ е I Ф Ь~ — а~с~ Очевидно, что для решения пятидиагональной системы потребуются две прямые и одна обратная прогонки. Различные этапы обобщенного алгоритма Томаса можно интерпретировать как серии операций (прямых прогонок), требуемых для приведения матрицы А к верхнетреугольному виду, После этого для решения уравнений (6.33) требуется обратная подстановка следующего общего вида: о, = с(" ,— с",с, +, — 1",о, (6.35) э 6.2.
Прямые методы для линейных систем 24б после чего с помощью обратных подстановок можно рассчитать Ч. Ясно, что чем более широкую ленточную структуру будет иметь матрица А, тем менее экономичным будет вышеописанный алгоритм. Если внутри самой ленты вкраплено достаточно большое число нулей, то, как правило, более экономичными оказываются некие альтернативные приемы обработки разреженных матриц.
(Зепи(паз, 1977а], хотя они вызывают больше трудностей при программировании. 5.2.5. Блочные трехдиагональиые системы Чт Ь, с, ,ь, (6.36) ас Ьс сс Ч ам 1 Ьм ~ сч 1 ам Ьм где аь Ьь с; — субматрицы размера М м,М, а Ч; и б; суть М- компонентные субвекторы. Число М соответсвует числу уравнений, записываемых в каждой точке сетки; например, для трехмерного течения вязкой сжимаемой жидкости М =Ь (п. 1!.6.3). Отсюда следует, что Ч; — субвектор решения, связанный с конкретной точкой сетки. Уравнение (6.36) — система„ состоящая из 1Ч блоков уравнений, причем каждый блок связанный с конкретной точкой сетки, включает в себя М уравнений.
Решение системы (6.36) идет по методике, очень близкой методике решения (6.28). Сначала трехдиагональная матрица Описание алгоритма Томаса для решения трехдиагональных систем дается в п. 6.2.2. Как правило, этот алгоритм следует применять в тех случаях, когда единственное исходное уравнение подвергается дискретизации с помощью неявного алгоритма ($7.2). Однако многие задачи гидроаэродинамики описываются системами уравнений ($ 10.2 и гл. 11). Попытки реализовать неявные алгоритмы приводят тогда, как правило, не к скалярным трехдиагональным системам уравнений,а к блочным трехдиагональным структурам. Однако алгоритм Томаса без труда распространяется и на действия в таких ситуациях. Блочно-трехдиагональная система уравнений, эквивалентная системе (6.28), может быть представлена в виде Гл.
6. Стационарные задачи 246 блоков из (6.36) преобразуется к верхнетреугольной форме за счет исключения субматриц аь По аналогии с (6.29) первый блок уравнений дает ч ч с1 =(Ь1) сь 41=(Ь!) Ви (6.37) тогда как для блока общего вида имеем Ь',=Ь, — а,с,'. „с,'=(Ь,') 1с,, д! =(Ь,'.) '(а,. — а,д,'. 1). (6.38) В выражениях (6.37) и (6.38) фигурируют явные выражения обратных матриц.
На практике более экономичный путь связан с решением М-компонентной субсистемы; например, уравнение Ь,'.сз = сз решается относительно с,'. После выполнения операций (6.37) и (6.38) уравнение (6.36) приобретает верхне- треугольную форму прн замене с; и д; нас„и бз и при замене Ь| на единичную матрицу 1. Второй этап, эквивалентный применению формул (6.31), требует обратных подстановок по формулам Как показывает подсчет, для реализации блочного алгоритма Томаса требуется 0(5АГМз/3) операций, что явно предпочтительно числу О((/чМ)з/3), требуемому для реализации полного исключения по Гауссу. Однако если бы оказалось возможным расщепить блочно-трехдиагональную систему (6.36) на М скалярных трехдиагональных систем, то полный подсчет показал бы 0(бл/М) операций, т.
е. грубо говоря, уменьшение времени счета в М'/3 раз при /ч' » М. В результате возникает стимул к тому, чтобы построить неявные алгоритмы, которые позволили бы части уравнений отделиться от системы в целом. Вариант блочного алгоритма Томаса, описанный выше, излагается по книге 11заасзоп, Ке!!ег, 1966). 6.2.6. Прямые алгоритмы решения уравнения Пуассона Уравнение Пуассона, так же как и связанное с ним уравнение Лапласа, достаточно часто встречается в гидроаэродинамике, чтобы имело смысл исследовать возможности разработки специальных процедур для решения дискретизированного уравнения Пуассона.
Трехточечная конечно-разностная дискретизация двумерного уравнения Пуассона приводит к системе уравнений типа (6.23), для которой матрица А имеет форму, показанную на рис. 6.1. Для случая однородной сетки (бх = бу), покрывающей прямоугольную область О < х < 1, О < у < 1, один из блоков 4 6.2. Прямые методы для линейных систем 24Т системы уравнений может быть записан в виде 1У»-~+ аЧ»+1У» -~ =Ьы й=1 М (640).
где М вЂ” число активных точек в направлении оси у. Каждый вектор Ч» содержит У неизвестных значений опм связанных с одной из линий сетки (у = у»), направленной вдоль оси х. Матрица д является трехдиагоиальной, имеет порядок У и отличные от нуля элементы 1, — 4, 1, относяшиеся соответственно к точкам (/ — 1, (с), (/, й) и (/+ 1, А).
Вектор Ь» содержит все те вклады в В из уравнения (6.23), которые вносит в правую часть единственная линия сетки, идущая параллельно оси х, т. е. у =у». При решении уравнения (6.40) особенно эффективными оказываются два метода — метод циклической редукции и метод разложения в ряд Фурье. Предпочтительная стратегия [Зхчаг(х(гацЬег, 1977] сводится к тому, чтобы вначале с помощью циклической редукции уменьшить размер М в уравнениях (6.40), после чего строить решение редуцированной системы с помошью разложения в ряд Фурье.
В процессе циклической редукции блок, расположенный в А-й (четной) строке уравнения (6.40), умножается на — ц и складывается с блоками в (й — !)-й и (й+ 1)-й строках с целью исключения У», и Ч»эн В результате получим 1Ч»»+ м Ч»+!У»»»=Ь», (6.417 где дш= 21 — цн = ( ~/21+ д) ( т/2! — и), Ь~»п= Ь», — аЬ»+ Ь,+г Этот процесс можно повторить требуемое число раз. После осуществления 1 редукций уравнение (6.40) приобретает форму !Ч»»+ а~пЧ+ 1У»,»» = Ьг, (6.42) где ы д' = — ~ (д — ~»!), Р» = 2 соз [(2й — 1) и/2'+']. » ! В принципе редукция может быть определена до тех пор, пока граничные условия не позволят определить У» и и Ч»+,т [хогг, 1970]. Затем с помощью обратных подстановок будут определены все промежуточные векторы Ч».
На практике, однако, оказывается более эффективным завершить данный процесс после 1 редукций, переключаясь на применение рядов Фурье. Гл. В. Стаииоиарине задачи Если на всех границах ставятся граничные условия Дирихле, то надлежащее представление в виде ряда Фурье будет следующим: о! а — — ~~ У, а в(п( + ) прн я=1, ..., М!, (6.43) тде М! — число линий сетки, параллельных оси х и оставшихся в уравнениях (6.42) после 1 редукций. Коэффициенты У,,» определяются путем решения трехдиагональных систем вида У, а, +»,У„»+У, »,— — Н, а при в=1, ..., М, (6.44) хде Н, »= + ~~ й(х;, у») и!п ( + ), (6.45) :величина й(хь у») — элемент вектора Ь!'! в уравнении (6.42) и Д,= — 2~1+ 2 в!па( — '.)] (6.46) .Быстрое преобразование Фурье [Соо1еу е! а1., 1970; Вг!дйаш, 1974] применяется для нахождения неизвестных величин в соотношениях (6.43) и (6.45).