Fletcher-1-rus (1185917), страница 41
Текст из файла (страница 41)
б.!5. Распечатка подпрограммы РАСТ. Гл. 6. Стационарные задачи ЗОВЙОВТ1ИЕ ЗОЬЧЕ(И,А,ЗРЧТ,В) В1НЕИЗ10И А(50,50),дРЧТ(50),6[50) ГОЙЧАЙВ ЕЫН1ИАТХОИ ИИ И-1 002 Е 1,НН ЕР к+1 Ь дРЧТ(Ю 3 В(Ы В(Ы В(Е) В(Е) 3 ВО 1 1 = ЕР,И В(1) В(1) + А(1.Е)ак СОИТХИ0Е СОНТ1кок ВАСЕ ЗОВЗТ1ТОТ10И ВО 4 ЕА 1,ИН ЕИ И-ЕА Е ЕН+1 В(х) В(Е)/А[к,Е) 3 — В (Е) ВО 3 1 1,П( ВГЙ) В(1) + А(1,Е)*Е СОИТ1ИОК СОИТ1ИНК ВГА) ВП)/А[1,1) Йетоки ЕИВ Рнс. 6.16. Распечатка подпрограммы ЯО[.ЧЕ. наибольший ведущий элемент аьа (РАСТ, строки 17 — 24).
Если наибольший ведущий элемент равен нулю, то программа выдает сигнал ошибки [ЛРЧТ(Х) = — Ц и осуществляется возврат к основной программе. В предположении о том, что наибольший ведущий элемент не равен нулю, осуществляется перестановка рядов (РАСТ, строки 39 — 41), чтобы переместить наибольший ведущий элемент аь а в положение (й, й). В каждом столбце 1 (1= А+ 1, У) коэффициенты (-й строки ((=й+ 1, У) модифицируютси (РАСТ, строки 38 — 46) в соответствии с формулой го(,ат а, ) —— а, [ — [ч — '~аа Р ч а,а) (6.25)) 1 г с 3 С 4 С 5 с 6 С 7 з с 9 С 10 С 11 Ьг 13 1 ° 15 16 17 16 19 ЗО 21 2 22 С 23 С 24 С 25 26 27 2В 29 30 31 32 3 33 4 34 35 36 ЗОЬЧЕЗ ВТЕКАЙ Зтктзн, Аак В АЗЗОНЕЗ А 16 ГАСТОЙ1ЗЕВ 1ИТО ЬЛ ГОЙИ (ВТ ГЛСТ) ЙЕТОЙИЗ ЗОЬОТ10И.
Х, 1И В $6.2. Прямые методы для линейных систем 237 (6.26) Массивы хранятся в памяти вычислительных машин по столбцам (индекс 1 в аь;). Поэтому наиболее целесообразно сделать так, чтобы самый внутренний цикл подпрограммы РАСТ (строки 43 — 45) реализовался по первому индексу. Подпрограмма БОьУЕ прежде всего преобразует элементы вектора В (строки 12 — 21) с целью взаимодействия с ведущей строкой и модификации элемента Ьь 1= й+ 1, )т', путем вычитания из него ведущей строки с соответствующим множителем, т. е. 5,=5,— ( — "") 5„ что эквивалентно операции (6.
25) в подпрограмме РАСТ. Окончательное решение получается с помощью последовательных подстановок, реализуемых посредством подпрограммы 50)Л'Е (строки 25 — 34). Это дает .,=(ь, т .,ю.)1'„. мм е т+1 Для экономии памяти решение У записывается иа место использованного массива и выдается уже как массив В. Реализация исключения по Гауссу в два этапа, соответствующих подпрограммам РАСТ и 50ьЧЕ, полезна в тех случаях, когда требуется строить решения для многих вариантов величины В в правой части уравнения (6.23).
Это связано с тем, что для выполнения подпрограммы РАСТ требуется 0(У') операций, тогда как 50ьЪЕ требует 0(АГ') операций. Обе эти подпрограммы, представленные в несколько расширенной форме, подвергаются более тщательному обсуждению в книге [Рогзурпе е1 а!., 1977, равд. 3.3). Как правило, при использовании спектральных методов (2 5.6) и панельных методов ($14.1) получаются матрицы А с плотной структурой. В противоположность этому методы конечных элементов и конечных разностей приводят к матрицам А разреженной структуры и, если используется методика расщепления (факторизации, $ 8.2), матрицам узколенточной структуры.
Если матрица А разреженная, но не обязательно ленточная, то исключение по Гауссу часто осуществляется при введении А в память в форме одного массива, содержащего значения ац вместе со связанным с иим индикаторным массивом 1А, обеспечивающим информацию о положении точки (й 1). Основная трудность, возникающая при программировании процесса гауссовского исключения с разреженными матрицами, связана с вероятностью заполнения матрицы А. Иначе говоря, в процессе Га.
б. Стационарные задачи исключения элементы матрицы А, первоначально равные нулю, становятся отличными от нуля. Различные методы надлежащего представления А и !А обсуждаются в книгах (3епп(пдз, 1977а, гл. 5; 1)цП, 1981], последняя из которых предлагает более высокий уровень анализа. Существуют пакеты программ для реализации гауссовского исключения при разреженных матрицах; можно рекомендовать разработанные в Харуэлле программы МА28 и т. д., рассмотренные в вышеназванной книге Даффа.
б.2.2. Трехдиагональные системы: алгоритмы Томаса Использование трехточечных конечно-разностных формул или конечных элементов с линейной интерполяцией приводит после расщепления ($8.2) к появлению матриц А в уравнении (6.23), имеющих трехдиагональную структуру. Использование конечно-разностных схем более высокого порядка или конечных элементов более высокого порядка приводит к ленточной структуре А, где ширина области ненулевых элементов будет больше. Алгоритм Томаса пригоден для решения уравнения (6.23) в случае трехдиагональной матрицы А.
Обобщение алгоритма Томаса на случай пятидиагональной матрицы А описывается в п. 6.2.4. Если речь идет о системах уравнений, то матрица А имеет, как правило, блочную (трехдиагональную) структуру. Порядок действий в этом случае рассматривается в п. 6.2.5. Если ненулевые элементы располагаются вблизи главной диагонали, то полезно рассмотреть те варианты исключения по Гауссу, которые используют преимущества ленточной структуры А. Один из примеров такого варианта соответствует рассмотренной в 9 9.3 задаче о конвекции — диффузии. Если воспользоваться формулами с центральными разностями, то при обозначениях данного параграфа получается следующий алгоритм: — (1+ 0.цыц) о,, + 2о; — (1 — 0.5Я,.п) от,, —— О, (6.27) при повторении которого по отношению к каждому узлу по- лучим Ь, с~ а, Ь, с, (6.28) ас Ьс с~ а„ , Ь„ , с„ , а„ Ьн $ б.2.
Прямые методы для линейных систем хх ххх Прогонка аиерел х х х х х х хх рт + 1+ 1+ 1+ Прогонка Рис. 6.!7. Алгоритм Томаса для решения трехдиагональной системы урав- нений. преобразуются к виду — (/ ! 1 с', 1 с' 1 с', Ог 1 с'„, 1 и и р г(' т. е. из них исключаются коэффициенты аь тогда как коэффициенты Ь; нормализуются к единичным значениям. Для первого уравнения этой цели служат формулы с с !Ьи г( г2 /Ьр (6.29) а для произвольно расположенного уравнения— / сг, «1г — агог с', = йг — агг,, Ь; — агс; (6.30) где аг = — (1+ 0.5 1тсеи), Ьг = 2, с; = — (1 — 0.5 ггееи) Ненулевые значения 4 связаны с членами типа источника или, как с 4 и г(и, с граничными условиями.
Все элементы матрицы А, кроме показанных выше, равны нулю. Можно отметить изменение некоторых обозначений по сравнению с п. 6.2.1, в частности относящихся к Ьь Алгоритм Томаса для решения уравнений (6.28) состоит из двух частей (см. схему на рис. 6.17). Сначала уравнения (6.28) Гл. 6. Стациоиариые зааачи Модификация уравнений по формулам типа (6.30) проводится в процессе прямой прогонки (рис. 6.17). Второй этап алгоритма сводится к обратной подстановке (обратная прогонка на рис. 6.17), т. е.
к использованию формул он=де, о,=дс — о!,.1си (6.31) Алгоритм Томаса чрезвычайно экономичен; он требует для своего исполнения всего 5У вЂ” 4 операций (умножений и делений). Однако во избежание плохой обусловленности (и, следовательно, накопления ошибок округления) необходимо, чтобы (о;1,)' )а;! + (с;!. Как правило, использование расщепления при решении многомерных задач Я 8.2) приводит к трехдиагональным системам уравнений, которые могут эффективно решаться при помощи алгоритма Томаса. 6.2.3. ВАХРАС/ВАЯЯОЬ: исключение ло Гауссу для узколенточных матриц В тех случаях, когда матрица А имеет узколенточную структуру, для реализации исключения по Гауссу подходят подпрограммы ВАрчРАС и ВА)ч(ЯОЬ.
Подпрограмма ВА(ч(РАС (рис. 6.18) выполняет прямую прогонку или несколько повторяющихся прямых прогонок, чтобы придать матрице А верхне- треугольную форму. Подпрограмма ВА)ч8ОЬ (рис. 6.19) преобразует правую часть (6.23), а именно вектор В, и полученная система уравнения решается путем обратной подстановки. Для случая 1)чТ = 2 подпрограммы ВАХРАС и ВАЫЯОЬ составлены так, чтобы воспользоваться особым видом структуры матрицы А, полученным при использовании одномерных квадратичных элементов.
А именно матрица А оказывается попеременно то «трехдиагональной» (три соседних ненулевых элемента в одной строке), то «пятидиагональной» (пять соседних ненулевых элементов в одной строке), что соответствует уравнениям, связанным с узлами в серединах или в углах элементов соответственно (см. рис. 5.!О). В этом случае процесс факторизации (ВАХРАС) состоит нз двух прогонок. Первая прогонка воздействует только на «пяти- диагональные» строки с целью исключения лишнего коэффициента аа, е,. Вторая прогонка приводит матрицу к верхне- треугольной форме.
Конкретный вид реализации зависит от того, связаны ли первое и последнее уравнения из системы (6.23) с узлами в серединах или в углах элементов. Граничные условия Дирихле соответствуют срединным узлам; граничные же условия Неймана соответствуют угловым узлам. Читателю может оказаться полезным изучить часть программы, относящуюся ОПЫИ$10И В(5,65) 1Р(1ИТ .ЕО. 2)ООТО 2 1ИТ 1 ИР И-1 0013 1,ИР ЗР 3+1 ВИ,ЗЮ В(2,39)(В(3,3] В(э,зю В(3,39] - В(2,39]еВ(4,3] 1 соитпвк АЕПЫИ г ии и(г ВО 5 1 1,2 33 *3 -1 ВО 4 д дЗ,ИВ ЗА 2*ы-ц 1Р(1 жо.
цоото з зв ЗА+ г Вн„)В) - ВСЦЗВ)(Ва.зВ-Ц В(2,3В) В(2,3В] - В(1,3$]*В(З,ЮВ-Ц 9(з.зв) ° в(з,зв) - ви,зв]*в(а,эв-ц оото а к высказанным положениям и начинающуюся с 1].'1Т = 2, после прочтения п. 6.2.4. Примером использования подпрограмм ВАРГАС/ВА11801 как в случае трехдиагональных, так и попеременно трехдиаго- 16 К Флетчер, т. ! 2 зс 4 С 5 С ЗС 7 в 9 С 10 С 31 С 12 13 14 15 16 17 за 19 С 2О С 21 С 22 с 23 24 25 26 27 2В 29 С зо с 31 С Зг зз за З5 36 37 С зв с 39 С 40 41 аг аэ 44 45 46 ат 4$ 49 50 6 6тд Прямые методы для линейных систем Звввоотзие вм(РЯС (В, и, зит] РдстОА1$ез ВАЯВ ивта11 АА1$1ио Ряои ьгикда ОА ооввадт1с еьеиеитз 1ито ьщ А!икая кькпкитз ТА191АО(ывь Зтзткк 1 1, Р1АЗТ РАЗЗ, АЕВОСЕ ТО ТЯ101АООИАЕ 1 2. зксоив РА$$, Аквоск то ВРРВА ТА13иепдА гит 2, аозвявтзс кьппитв Рп(тввздооидь зтвтки Аззоиез Р1азт ЕООА710и Ровней Ат И10$1ве иове З 3$ ЗА + З В(2,3В-Ц ° В(2,ЮВ-1)/В(э,дВ-2) В(3,3В-Ц В(З,ЗВ-Ц вЂ” В(2,3В-Ц *В(4,3В-2) 1Г ЫВ .ОТ.
И)СОТО ° В(г.ав) В(2,39)/В(3,3В-Ц В(3,3В) е В(3,3В) - В(2,39)"В(4,3В-ц В(4,3$) В(а,дв) - В(2,3$)*В(5,3В-Ц 4 СОИТ1ИОЕ 5 СОИТ1ИПЕ ЯЕТОЯИ ЕИВ Рис. 6.18. Распечатка подпрограммы ВА]ЧРАС. Гл. 6. Стационарные задачи нальных/'пятидиагональных матриц является решение задачи Штурма — Лиувилля (п. 5.4.2).
Полезная проверка правильности составления подпрограмм решения уравнений (6.23) состоит в том, чтобы после определе- 1 г з с 4 С 5 С 6 5 С ВС )О С 11 эг 13 14 15 16 17 15 19 го гэ с 22 С гз с 24 С 25 С 26 С 27 2В 29 ЗО Эг 32 33 Зб Эб 36 Зт ЗВ 39 40 41 42 43 44 45 46 47 45 эовхонтэме вамэоь(м,х,в,и,эит) шкэ ь.п гастоиэхтэои то эоьтк ток х, Озгки и 01МЕИ510И Е(65),Х(65),В(5,65) 1Р(1МТ .ЕО. 2)СОТО 3 1ИТ = 1, ТЫШАООИАЬ 575ТЕН ИР И-1 ОО 1 3 = 1,ИГ дг Э+1 1 В(ЭР) В(ЭР) — В(г,зг)*В(д) Х(М) = В(и)/В(з,и! ОО 2 д 1,ИР ЭА = М " д Х(Л) ЬЧ(Л! - В[4,Л) Х(ЭЬ1))/В(Э,ЭА) 2 СОМТ1ИОЕ ВЕТОАИ 1МТ 2, РЕИТАВ1ЛСОИАЬ 515ТЕН 155ОНЕ5 11ВЗТ ЕООАТ10И ГОВНКО АТ Н1551ОЕ ИОВЕ 5ет 1Вс ° 0 эг ( А5т еонлт10и гокмен Ат нэоээое иоэе эет 1Вс 1 19 ьаэт еООАт10Я РОВнеп Ат соемеа иоое 31ВС 1 ИН М/2 эг(г ин .ко.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
