Содержание к диссертации
Введение
Глава 1. Методы неполной факторизации для сеточных систем дифференциальных уравнений 14
1.1. Алгебраические особенности сеточных систем уравнений и постановка задачи 14
1.1.1. Основные понятия и обозначения 14
1.1.2. Особенности сеточных систем при неявных аппроксимациях параболических уравнений 15
1.1.3. Постановка двумерной параболической краевой задачи с постоянными коэффициентами 19
1.1.4. Трехмерная эллиптическая краевая задача 20
1.2. Применение итерационных методов для решения сеточных систем 24
1.2.1. Общие свойства и анализ итерационных алгоритмов . 24
1.2.2. Рассмотрение метода сопряженных градиентов 32
1.2.3. Предобусловленный метод сопряженных градиентов . 35
1.2.4. Применение метода неполной факторизации 39
1.2.5. Корректность и устойчивость неполной факторизации . 44
Глава 2. Применение методов неполной факторизации к решению двумерных параболических краевых задач и их матричный анализ 47
2.1. Анализ параболической краевой задачи с постоянными коэффициентами 47
2.1.1. Оценки элементов точных факторизации 47
2.1.2. Получение улучшенных оценок элементов точных факторизации 57
2.2. Оценки неполных факторизации модельной параболической краевой задачи 63
2.2.1. Теорема об оценках норм матриц G 63
2.2.2. Вспомогательные результаты для неполной блочной факторизации 65
2.2.3. Вспомогательные результаты для неполной факторизации Холесского 68
2.2.4. Теоремы об оценках чисел обусловленности предобуславливателей матрицы В 71
2.3. Методы неполной факторизации для параболических краевых задач с переменными коэффициентами 73
Глава 3. Решение трехмерных эллиптических краевых задач методом неполной факторизации 76
3.1. Постановка эллиптической краевой задачи и оценки блочных факторизации 76
3.1.1. Постановка задачи 76
3.1.2. Оценки элементов матрицы Gjj 77
3.1.3. Оценки элементов матрицы компенсации S(j 84
3.2. Оценки неполных факторизации трехмерных эллиптических краевых задач 96
3.2.1. Теорема об оценках норм матриц Gin 96
3.2.2. Вспомогательные результаты для неполной блочной факторизации 97
3.2.3. Теорема об оценках предобуславливателей типа неполной блочной факторизации для эллиптических краевых задач 100
Глава 4. Численное решение эллиптических и параболических краевых задач 101
4.1. Реализация методов решения на основе алгоритмов работы с разреженными матрицами специального вида 101
4.1.1. Схемы хранения разреженных матриц
4.1.2. Алгоритмы работы с разреженными матрицами 105
4.1.2. Реализация методов неполной факторизации 111
4.2. Численный эксперимент 114
Литература 119
Приложение 1 127
- Особенности сеточных систем при неявных аппроксимациях параболических уравнений
- Получение улучшенных оценок элементов точных факторизации
- Теорема об оценках предобуславливателей типа неполной блочной факторизации для эллиптических краевых задач
- Реализация методов решения на основе алгоритмов работы с разреженными матрицами специального вида
Введение к работе
При моделировании явлений или процессов, как природно-естественного, так и технического характера часто возникает необходимость решения сложных задач математической физики эллиптического и параболического типа с использованием сеточных методов.
Зачастую при рассмотрении таких задач возникает необходимость решать их в областях достаточно сложной формы. Именно поэтому в последнее время для их решения наиболее эффективным средством стали итерационные методы неполной факторизации. Эти методы хорошо зарекомендовали себя на практике и легли в основу ряда пакетов прикладных программ. Их главными достоинствами являются рекордная практическая экономичность, высокие теоретические оценки скорости сходимости итераций и широкие возможности конструирования адаптивных алгоритмов для разных классов задач.
Большое значение при решении таких задач математической физики придается градиентным итерационным методам, в частности, методу сопряженных градиентов [39],[56]. Этот метод пригоден для областей достаточно сложной формы, когда алгоритмы типа быстрого преобразования Фурье [45], переменных направлений для прямоугольных областей [66] становятся не эффективными. Он является быстросходящимся и применяется для решения возникающих при аппроксимации исходной дифференциальной задачи систем линейных алгебраических уравнений (СЛАУ) произвольного вида. Достаточно только, чтобы матрица этой системы была положительно определенной и симметричной. Его привлекательность обусловлена, однако, не только его весьма высокой скоростью сходимости (это достигается и у Чебышевского метода с оптимальным набором параметров [39]), но и тем, что метод не требует дорогостоящей априорной информации о матрице системы (граница спектра, число обусловленности).
Вместе с тем, для систем высокого порядка, метод сопряженных градиентов сходится достаточно медленно. Поэтому применяются различные способы ускорения метода сопряженных градиентов. В настоящей работе используется ускорение метода сопряженных градиентов с помощью неполных факторизации.
Метод сопряженных градиентов в сочетании с алгоритмами работы с разреженными матрицами [57], [81] позволяет решать дифференциальные задачи достаточно большой размерности, порядка десятков тысяч. За последние десятилетия итерационные методы неполной факторизации, применяемые совместно со спектральными или вариационными принципами ускорения сходимости последовательных приближений, стали наиболее эффективным средством для решения многих дифференциальных задач.
Методы неполной факторизации были открыты в 50-е годы Н.И. Булеевым [25]. Методы Н.И. Булеева предназначались для итерационного решения пятиточечных и семиточечных конечно-разностных уравнений, аппроксимирующих двумерные и трехмерные краевые задачи математической физики. Принципы их конструирования можно считать формальным обобщением безытерационного метода прогонки для решения одномерных краевых задач, основанного на разложении трехдиагональной матрицы в произведение нижней и верхней треугольных матриц. Существенно новым моментом для итерационных процессов явилось введение приема, названного Н.И. Булеевым принципом компенсации. Он заключается в добавлении к уравнению таких членов, которые приводят к взаимному исключению ошибок итерационных приближений, если они являются гладкими функциями координат. А это как раз представляется естественным для задач математической физики.
Несмотря на отсутствие теоретических обоснований в рассматриваемых работах, в многочисленных экспериментах по численному решению практических задач методы неполной факторизации обнаруживали высокую эффективность. Особые преимущества эти алгоритмы имели при их использовании с учетом свойств решаемых задач. Тем не менее в связи с недостатком теоретических выкладок полученные результаты еще долгое время не привлекали пристального внимания отечественных специалистов, а исследования в России долгое время продолжались исключительно Н.И. Булеевым и его учениками(В.П. Гинкиным [32],[33], В.К. Артемьевым [16], [17] и некоторыми другими).
Зарубежными учеными методы неполной факторизации неоднократно переоткрывались в последующие годы и получили дальнейшее развитие в работах О. Axelsson'a [1]-[6], R. Beauwens'a [7], X. Van der Vorst'a [14], J. Gol-ub'a [10], P.Concus'a, J. Meurant'a [9], T. Manteuffel'a [12], Y. Notay'a [13], T. Chan'a [8], P. Vassilevski [15] и многих других авторов. Они стимулировали новые исследования и привнесли ряд интересных результатов.
Одно из важных направлений работ связано с исследованием итерационных процессов на основе использования предобуславливающих матриц в виде приближенного разложения матрицы исходной системы на треугольные множители. Для симметричных уравнений это дает семейство алгоритмов, называемых разложением Холесского [39],[56]. Известно, что вычислительная сложность разложения Холесского для разреженных матриц высокого порядка определяется главным образом степенью заполнения вычисляемых треугольных матричных множителей ненулевыми элементами. Отсюда появляются различные возможности для конструирования неполных факторизации, основанных на принудительном аннулировании некоторых новых матричных элементов в соответствии с какой-либо выбранной стратегией. В результате формируется класс итерационных процессов, отличающихся как скоростью сходимости, так и трудоемкостью реализации отдельного итерационного шага.
Другим важным моментом явилось построение алгоритмов, основанных на вычислении ленточной части матрицы, обратной к ленточной матрице со свойством диагонального преобладания. Это позволило получать хорошие приближенные факторизации для М-матриц и на их основе конструировать эффективные неявные итерационные методы. Такие подходы были рассмотрены впервые в работах В.П. Ильина [44],[41], Н.С. Бахвалова [18], С.А. Сандер [61], O.Axelsson'a [5], G. Golub'a и G. Meurant'a [9].
К настоящему времени исследования по решению дифференциальных уравнений в частных производных с помощью методов неполной факторизации пока не получили полного освещения в учебной и монографической литературе российских специалистов. Однако, следует отметить, что одним из первых систематических изложений современных результатов в данной области является монография В.П. Ильина [39]. Этой тематике посвящены также последние работы И.А. Блатова [20]-[24], А.Ю. Еремина [38], И.Е. Капорина [46], Л.Ю. Колотилиной [48].
Главными теоретическими вопросами решения краевых задач методами неполной факторизации являются обоснование самого существования и устойчивости применяемых матричных преобразований, определение условий сходимости итерационных процессов и получение оценок скорости сходимости итераций. Первые вопросы за последние годы уже практически решены, а в последних имеется еще немало открытых аспектов. Эксперименты в этой области зачастую опережают теорию и ставят ряд таких исследовательских проблем, для решения которых не хватает существующего теоретического аппарата.
Особое место здесь занимают проблемы с разреженными матрицами специальной структуры, возникающими из сеточных аппроксимаций краевых задач и предъявляющими особенно жесткие требования к быстродействию и оперативной памяти компьютера. По данной тематике изучались работы В. Базова и Дж. Форсайта [29], Г.И. Марчука [53], [54], И.А. Блатова [21], А.А. Самарского и Е.С. Николаева [59], О.А. Ладыженской [50], В.Б. Андреева [60], А.Н. Тихонова [64].
Одной из основных трудностей при разработке и, особенно, при обосновании вычислительных алгоритмов с разреженными матрицами является то, что многие основные вычислительные операции над разреженными матрицами (обращение, переход к треугольным и ортогональным факторизациям) не сохраняют разреженность. Но вычислителями-практиками давно было замечено, что, несмотря на потерю разреженности, большинство вновь возникающих ненулевых элементов часто оказываются малыми по абсолютной величине, и замена их нулями приводит к хорошим разреженным аппроксимациям заполненных матриц, к которым применимы и эффективны вычислительные алгоритмы типа методов неполной факторизации. Теоретическое обоснование этого феномена (за исключением некоторых специальных случаев) практически отсутствовало, что не позволяло обосновать ряд существующих алгоритмов и затрудняло разработку новых. В связи с этим актуальным является любой вклад в построение соответствующей теории. В настоящем времени этой тематикой стали заниматься Y. Notay [13], И.А. Блатов [21]-[23], М. Р. Ларин [11], чьи последние работы были внимательно изучены и учтены при исследовании.
Важнейшим направлением исследований решения задач математической физики при помощи методов неполной факторизации является расширение типов систем уравнений, для которых обосновываются и эффективно применяются специальные или универсальные алгоритмы.
Одними из таких интересных в математическом плане задач являются дифференциальные уравнения в частных производных эллиптического и параболического типа, аппроксимации которых приводят к алгебраическим системам. Здесь оказываются актуальными алгоритмы, учитывающие свойства исходных дифференциальных задач.
Актуальность темы обусловлена также приложениями полученных теоретических результатов к важным классам дифференциальных задач.
В постановке численного эксперимента помогли работы Х.Гулда и
Я.Тобочника [36], А.Н. Бубенникова [28], Н.Н. Лебедева, И.П. Скальской, Я.С. Уфлянд [51], Б.С. Польского [58], С.Г. Мулярчика [55], А.В. Лыкова [52].
Таким образом, получение и обоснование эффективных алгоритмов решения параболических и эллиптических краевых задач на основе методов неполной факторизации представляет большой как теоретический так и практический интерес. Поскольку, во-первых, такие задачи большой размерности возникают во многих отраслях науки и техники, а итерационные методы неполной факторизации позволяют достаточно быстро находить их решение в условиях ограниченных вычислительных ресурсов. Во-вторых, вопросы теоретического обоснования применения данных методов для дифференциальных уравнений эллиптического и параболического типа еще недостаточно хорошо изучены.
Целью работы является изучение численного решения дифференциальных уравнений эллиптического и параболического вида с помощью методов неполной факторизации с различными видами предобуславливаниелей.
В соответствии с поставленной целью в работе анализируются разреженные матрицы сеточных систем уравнений, соответствующих исходным дифференциальным задачам, строится теория оценок их элементов, оценок элементов обратных матриц и элементов точных и неполных факторизации. Для соответствующих задач доказываются оценки спектральных чисел обусловленности, характеризующих скорость сходимости итерационного процесса. Полученная теория применяется для решения эллиптических и параболических краевых задач в областях достаточно сложной формы.
В работе используется аппарат теории обыкновенных дифференциальных уравнений, асимптотических методов, элементы анализа, линейной алгебры и функционального анализа, теории итерационных методов.
Диссертационная работа состоит из введения, четырех глав, четырех приложений и списка цитированной литературы, содержащего 81 наименование. Главы разделены на параграфы, параграфы - на подпараграфы. Нумерация подпараграфов сквозная. Формулы, замечания, предложения, леммы и теоремы нумеруются двойными числами, первое из которых - номер главы, а второе - номер формулы или утверждения внутри параграфа.
В первой главе приводится постановка изучаемых краевых задач эллиптического и параболического вида и рассматриваются особенности аппроксимации этих задач сеточными системами уравнений. Обосновывается применение итерационных методов, сравнительный анализ которых показывает эффективность использования предобусловленного метода сопряженных градиентов. В качестве предобуславливателей рассматриваются точечные и блочные неполные факторизации.
Во второй главе изучается матрица системы, получаемой в результате разностной аппроксимации двумерной параболической краевой задачи с постоянными коэффициентами. На основе анализа данной матрицы получены оценки элементов неполной факторизации Холесского и неполной блочной факторизации в зависимости от структуры допустимого заполнения предобуславливателя, при помощи которых доказаны теоремы об оценках чисел обусловленности, оказывающих решающее влияние на скорость сходимости итерационного процесса. С целью дальнейшего уточнения результатов проведен дополнительный анализ и получены улучшенные оценки элементов факторизации. Рассмотрена параболическая краевая задача с переменными коэффициентами и обосновано применение методов неполной факторизации, для которых справедливы полученные оценки.
В третьей главе рассматривается решение трехмерной эллиптической краевой задачи с помощью методов неполных факторизации. Доказаны оценки элементов соответствующей матрицы сеточной системы уравнений и обратной к ней. На основе этих оценок получены оценки числа обусловленности предобуславливателя типа неполной блочной факторизации, подтверждающие эффективность решения исходной задачи выбранным методом. В целях дальнейшего ускорения сходимости рассмотрена неполная блочная факторизация матрицы с диагональной компенсацией. Для элементов матрицы компенсации получены соответствующие оценки.
В четвертой главе предложены эффективные по скорости и требованиям к машинной памяти схемы хранения и алгоритмы работы с разреженными матрицами большой размерности, ориентированные на матрицы специального вида, изученные во второй и третьей главах. Приведены формулы, позволяющие оценить размерность эффективно решаемых задач в зависимости от доступной памяти. Проведены численные эксперименты по решению тестовых эллиптических и параболических краевых задач с известным решением. Результаты численных экспериментов подтверждают эффективность применения методов неполной факторизации и справедливость полученных во второй и третьей главах оценок.
Работа имеет как теоретическую, так и практическую направленность. Ее результаты могут быть использованы для дальнейшего исследования дифференциальных уравнений в частных производных с граничными условиями, являющихся моделями различных процессов, а также для написания ряда прикладных программ при решении задач математической физики. Полученные выводы позволят детально проанализировать имеющиеся плюсы и минусы широкого класса итерационных алгоритмов при решении рассмотренных дифференциальных задач, актуальных в многочисленных приложениях.
Таким образом на защиту выносятся:
1. Оценки точечных и блочных факторизации, отличающиеся от ранее известных существенно более высокой точностью.
2. Оценки числа обусловленности предобусловленной матрицы, позволяющие оценить скорость сходимости итерационного процесса в зависимости от структуры допустимого заполнения.
3. Метод решения параболических краевых задач с переменными коэффициентами с применением неполных факторизации.
4. Эффективные схемы хранения и алгоритмы обработки разреженных матриц специального вида большой размерности, ориентированные на матрицы, получаемые при решении рассматриваемых задач математической физики.
Работа над диссертацией выполнялась на кафедре "Вычислительная математика" Воронежского государственного университета.
Результаты работы докладывались на весенних Воронежских математических школах "Понтрягинские чтения" (1998, 1999гг.), на конференции "Математическое моделирование систем. Методы, приложения и средства"(Воронеж, 1998г), на конференции "Нелинейный анализ и функционально-дифференциальные уравнения" (МНК АДМ-2000)(Воронеж, 2000г),на конференции "Математика. Образование. Экология. Тендерные проблемы. "(Воронеж, 2000г), на семинарах профессоров И.А. Блатова(Самара) и В.Г. Задорожнего, на ежегодных научных конференциях профессорско-преподавательского состава Воронежского государственного университета.
Результаты диссертации опубликованы в 12 работах [69]-[80].
Автор выражает искреннюю благодарность своему научному руководителю, доктору физ.-мат. наук, профессору В.В. Стрыгину за постоянное внимание, обсуждение и ряд ценных замечаний в процессе исследования, а также доктору физ.-мат. наук, профессору И.А. Блатову.
Особенности сеточных систем при неявных аппроксимациях параболических уравнений
Таким образом, получение и обоснование эффективных алгоритмов решения параболических и эллиптических краевых задач на основе методов неполной факторизации представляет большой как теоретический так и практический интерес. Поскольку, во-первых, такие задачи большой размерности возникают во многих отраслях науки и техники, а итерационные методы неполной факторизации позволяют достаточно быстро находить их решение в условиях ограниченных вычислительных ресурсов. Во-вторых, вопросы теоретического обоснования применения данных методов для дифференциальных уравнений эллиптического и параболического типа еще недостаточно хорошо изучены.
Целью работы является изучение численного решения дифференциальных уравнений эллиптического и параболического вида с помощью методов неполной факторизации с различными видами предобуславливаниелей.
В соответствии с поставленной целью в работе анализируются разреженные матрицы сеточных систем уравнений, соответствующих исходным дифференциальным задачам, строится теория оценок их элементов, оценок элементов обратных матриц и элементов точных и неполных факторизации. Для соответствующих задач доказываются оценки спектральных чисел обусловленности, характеризующих скорость сходимости итерационного процесса. Полученная теория применяется для решения эллиптических и параболических краевых задач в областях достаточно сложной формы.
В работе используется аппарат теории обыкновенных дифференциальных уравнений, асимптотических методов, элементы анализа, линейной алгебры и функционального анализа, теории итерационных методов.
Диссертационная работа состоит из введения, четырех глав, четырех приложений и списка цитированной литературы, содержащего 81 наименование. Главы разделены на параграфы, параграфы - на подпараграфы. Нумерация подпараграфов сквозная. Формулы, замечания, предложения, леммы и теоремы нумеруются двойными числами, первое из которых - номер главы, а второе - номер формулы или утверждения внутри параграфа.
В первой главе приводится постановка изучаемых краевых задач эллиптического и параболического вида и рассматриваются особенности аппроксимации этих задач сеточными системами уравнений. Обосновывается применение итерационных методов, сравнительный анализ которых показывает эффективность использования предобусловленного метода сопряженных градиентов. В качестве предобуславливателей рассматриваются точечные и блочные неполные факторизации.
Во второй главе изучается матрица системы, получаемой в результате разностной аппроксимации двумерной параболической краевой задачи с постоянными коэффициентами. На основе анализа данной матрицы получены оценки элементов неполной факторизации Холесского и неполной блочной факторизации в зависимости от структуры допустимого заполнения предобуславливателя, при помощи которых доказаны теоремы об оценках чисел обусловленности, оказывающих решающее влияние на скорость сходимости итерационного процесса. С целью дальнейшего уточнения результатов проведен дополнительный анализ и получены улучшенные оценки элементов факторизации. Рассмотрена параболическая краевая задача с переменными коэффициентами и обосновано применение методов неполной факторизации, для которых справедливы полученные оценки.
В третьей главе рассматривается решение трехмерной эллиптической краевой задачи с помощью методов неполных факторизации. Доказаны оценки элементов соответствующей матрицы сеточной системы уравнений и обратной к ней. На основе этих оценок получены оценки числа обусловленности предобуславливателя типа неполной блочной факторизации, подтверждающие эффективность решения исходной задачи выбранным методом. В целях дальнейшего ускорения сходимости рассмотрена неполная блочная факторизация матрицы с диагональной компенсацией. Для элементов матрицы компенсации получены соответствующие оценки.
В четвертой главе предложены эффективные по скорости и требованиям к машинной памяти схемы хранения и алгоритмы работы с разреженными матрицами большой размерности, ориентированные на матрицы специального вида, изученные во второй и третьей главах. Приведены формулы, позволяющие оценить размерность эффективно решаемых задач в зависимости от доступной памяти. Проведены численные эксперименты по решению тестовых эллиптических и параболических краевых задач с известным решением. Результаты численных экспериментов подтверждают эффективность применения методов неполной факторизации и справедливость полученных во второй и третьей главах оценок.
Работа имеет как теоретическую, так и практическую направленность. Ее результаты могут быть использованы для дальнейшего исследования дифференциальных уравнений в частных производных с граничными условиями, являющихся моделями различных процессов, а также для написания ряда прикладных программ при решении задач математической физики. Полученные выводы позволят детально проанализировать имеющиеся плюсы и минусы широкого класса итерационных алгоритмов при решении рассмотренных дифференциальных задач, актуальных в многочисленных приложениях.
Получение улучшенных оценок элементов точных факторизации
Для эффективности алгоритма (1.60) существенным является требование простоты решения системы Mr = г. С другой стороны, для того, чтобы предобуславливание было эффективным, желательно, чтобы матрица М хорошо аппроксимировала матрицу А. Ясно, что эти два требования вступают в противоречие.
Критерий близости матрицы М к А можно определить, рассмотрев соотношение которое следует из (1.56) и (1.59). Видно, что матрица М 1А подобна матрице А, и поэтому число обусловленности предобусловленной матрицы равно отношению наибольшего собственного значения матрицы М 1А к ее наименьшему собственному значению:
Таким образом, поскольку целью предобуславливания является уменьшение числа обусловленности матрицы А настолько, насколько это возможно, критерием близости матрицы М к А будет служить малость отношения наибольшего собственного значения матрицы М 1А к наименьшему.
Рассмотрим теперь способ проверки сходимости итераций. Естественный критерий rfc+12 є требует дополнительного вычисления на каждой итерации скалярного произведения (rfc+1,rfc+1), так как в алгоритме оно не фигурирует. Вместо него появляется величина (f 1 4 1) = (M_1rfc+1,rfc+1), которая представляет собой скалярное произведение, определенное при помощи матрицы М-1. Эта величина может быть использована для проверки сходимости, однако при этом условие сходимости может оказаться выполненным, в то время как норма rfc+1J2 еще не достигает желаемой малости. Поэтому в случае положительного результата проверки ( + 7- +1) є необходимо выполнение вспомогательной проверки значения rfc+12.
Применение метода неполной факторизации Важным подходом к предобусловливанию матрицы А является использование неполных факторизации. Если матрица А представлена в обычном точечном виде, будем рассматривать неполную факторизацию Холесского где L - нижнетреугольная матрица, R - некоторая ненулевая матрица. Если выполнить полную факторизацию Холесского A = LLT [31], множитель L оказывается гораздо менее разреженным, чем матрица А, из-за возникающего заполнения. Неполная факторизация позволяет сохранить структуру разреженности матрицы А или добавить незначительное количество ненулевых элементов по некоторому принципу, что в несколько раз сокращает время ее выполнения (полная факторизация представляет собой прямой метод решения системы, и если бы мы могли ее выполнить, не было бы смысла рассматривать итерационные методы и, в частности, метод сопряженных градиентов).
Пусть S = {(г,і)} - заданное множество индексов, где і и j принимают значения между 1 и JV. Данное множество и определяет структуру разреженности матрицы L. На его основе можно сформулировать алгоритм неполной факторизации Холесского: hi — va« 5Zfc=i hk В связи с необходимостью вычисления квадратных корней в последней строке алгоритма (1.64), неполная факторизация возможна только для симметричных положительно определенных матриц, которые являются //"-матрицами [56]. Если это условие не выполнено, используются дополнительные приемы. Пусть на г-м шаге алгоритма (1.64) получилось Ц% : 0) т-е- либо требуется извлечение квадратного корня из отрицательного числа, либо матрица L вырождена. Тогда можно заменить \\ на некоторое положительное число, например, на значение предыдущего диагонального элемента lf_u_i, или на величину (Е г)2, что позволяет обеспечить диагональное преобладание в новой г -й строке матрицы L. Вместо матрицы А можно также использовать матрицу А = А + а/, где а 0 выбирается так, чтобы матрица А имела строгое диагональное преобладание. Тогда А является і7-матрицей, и неполная факторизация (1.64) выполнима. В зависимости от выбора множества S можно получить конкретную реализацию неполной факторизации Холесского. Если S = {(i,j) : ( ф 0}, то матрица L имеет точно такую же структуру разреженности, как и матрица А (принцип отсутствия заполнения). Этот принцип можно ослабить. Например, если матрица А является пятидиагональной (это соответствует решению эллиптических и параболических задач на прямоугольной области), можно использовать IC(h, &)-принцип, допускающий заполнение h — 1 дополнительных диагоналей, смежных с первой поддиагоналыо, и к — 1 дополнительных диагоналей, смежных со второй поддиагоналыо. Таким образом, допускается наличие h+k диагоналей, расположенных ниже главной в матрице L. Если матрица А имеет блочно-трехдиагональный вид А = tridiag{—Dp,Epi—Fp}, 1 р п, что имеет место при решении параболических и эллиптических краевых задач на произвольных областях, будем рассматривать неполную блочную факторизацию, основанную на полной факторизации
Теорема об оценках предобуславливателей типа неполной блочной факторизации для эллиптических краевых задач
Рассмотрим матрицы систем, возникающих при разностной аппроксимации двумерных параболических уравнений (1.2)-(1.12) и трехмерных эллиптических уравнений (1.13). Составим для них наиболее эффективные схемы хранения элементов, характеризующиеся наименьшим объемом требуемой памяти. Это возможно сделать, учитывая специальный вид этих матриц. В целях достижения наибольшей экономии памяти будем отдельно рассматривать матрицы для различных задач параболическиого и эллиптического типа. Соответственно будем выделять схему хранения для матрицы каждой задачи.
Двумерная параболическая краевая задача с постоянными коэффициентами {к1(х,у) = 1, к2(х,у) = 1 на квадратной области.Для этой задачи матрица В имеет предопределенный вид В = tridiag{—Di, НІ, — i }, где Hi = tridiag{—l, 4+ г, —1}, Di = Fi = I, а I - единичная матрица. Таким образом, для решения этой задачи фактически вообще нет необходимости хранить элементы матрицы, так как они для любой строки являются постоянными. В результате экономится минимум 5iV ячеек памяти (N -количество неизвестных, может исчисляться десятками тысяч), необходимых для хранения элементов по диагоналям.
Двумерная параболическая краевая задача с постоянными коэффициентами на произвольной области. В отличие от предыдущей задачи матрицы Di и Fi уже не являются диагональными, но для них только один элемент в каждой строке или столбце соответственно может равняться 1. Поэтому, учитывая симметричность матрицы 5, для нее достаточно хранить один целочисленный вектор Д размерности N — щ (щ - размерность блока Hi), элементы которого содержат расстояние от главной диагонали столбца с единичным элементом для каждой строки блоков Di и одновременно строки с единичным элементом для каждого столбца блоков Fj_i, г = 2,..., п. Это расстояние равно разности последовательных номеров узлов сетки, находящихся на соседних линиях г и одной линии j. Дополнительно необходимо хранить вектор {щ}, содержащий размерности каждого блока Щ. Учитывая, что для хранения целого положительного числа, большего 65535, требуется 4 байта машинной памяти, необходимый для хранения матрицы объем памяти составляет
Схема 3. Двумерная параболическая краевая задача с переменными коэффициентами на квадратной области.Для этой задачи блоки D{ и F{ диагональные, однако значения всех ненулевых элементов блоков JD;, НІ, Fi зависят от строки и столбца i,j и вычисляются по формулам (1.5). В принципе элементы матрицы В можно было бы не хранить, как и в схеме 1, однако при этом сильно возрастает сложность вычисления элементов (которые в схеме 1 суть константы), и с точки зрения вычислительной эффективности эти элементы лучше рассчитать один раз и хранить по схеме типа диагональной. Для этого необходимо 3 вектора: d размерности N для элементов главной диагонали, / размерности N — п для нижней поддиагонали и a (N — п) для еще одной ненулевой поддиагонали. Данная схема отличается от диагональной тем, что в векторе I не хранятся нулевые элементы, возникающие на стыке блоков Н{, поэтому его размерность равна N — п, а не N—1. Элементы верхних ненулевых поддиагоналей можно определить исходя из симметричности матрицы В. Учитывая, что для хранения вещественного числа двойной точности требуется 8 байт машинной памяти, необходимый для хранения матрицы объем памяти составляет байт.
Схема 4- Двумерная параболическая краевая задача с переменными коэффициентами на произвольной области. Это самый сложный случай, отличающийся от предыдущего тем, что матрицы D{ и F; не являются диагональными. Поэтому необходимо применить смешанную схему хранения: элементы главной диагонали и смежной с ней нижней (верхней) поддиагонали хранятся в векторах d размерности N и N — п соответственно, дополнительно хранится целочисленный вектор А размерности N — п\, элементы которого содержат расстояние от главной диагонали столбца с ненулевым элементом для каждой строки блоков Di и одновременно для каждого столбца блоков Fj_i, і = 2,..., п. Сами ненулевые элементы блоков D{ (Fi) хранятся еще в одном векторе а размерности N — п\. Аналогично схеме 2 хранится также целочисленный вектор {щ} размерности п. Тогда необходимый для хранения матрицы объем памяти составляет байт.
Трехмерная эллиптическая краевая задача с постоянными коэффициентами на кубе.Для этой задачи матрица А имеет предопределенный вид A = fivediag{—Qij,—Dij,EiJT—Fij,—Rij} Eij = tridiag{—1,6,—1}, D(j = Рц = Qij = R%j = I- Таким образом, для решения этой задачи, так же как и в первом случае, нет необходимости хранить элементы матрицы, так как они для любой строки являются постоянными.
Трехмерная эллиптическая краевая задача с постоянными коэффициентами на произвольной области. В отличие от предыдущей задачи матрицы Dij и Fij, Qij и Щ уже не являются диагональными, но для них только один элемент в каждой строке или столбце соответственно может равняться 1. Поэтому, учитывая симметричность матрицы А, для нее достаточно хранить один целочисленный вектор А размерности N — E"=i j=2 lij элементы которого содержат расстояние от главной диагонали столбца с единичным элементом для каждой строки блоков Dij и одновременно столбца блоков F -\. Аналогично необходим еще один вектор Д(2) размерности N — s\ для хранения индексов единичных элементов в блоках Qi и Ri-i. Необходимый для хранения матрицы объем памяти составляет
Реализация методов решения на основе алгоритмов работы с разреженными матрицами специального вида
Методы неполной факторизации для решения эллиптических и параболических краевых задач реализованы на языке программирования Си++ с использованием технологии объектно-ориентированного программирования [63] в среде Windows NT 4.0.
Реализация основана на двух базовых классах: матрица (CMatrix) и вектор (СVector), для которых определены операции доступа к элементам, присвоения, умножения на число, умножения матрицы на вектор. Для векторов определена операция скалярного произведения, а для матриц -произведение матрицы на матрицу и обращение. Кроме того, для матриц имеется возможность получения столбцовых индексов и значений ненулевых элементов каждой строки.
Все перечисленные операции с матрицами определены как виртуальные методы класса CMatrix, что позволяет в конкретной реализации схемы хранения разреженной матрицы использовать соответствующий этой схеме алгоритм. С другой стороны, программный код нахождения неполных факторизации и решения СЛАУ предобусловлениым методом сопряженных градиентов не зависит от конкретной схемы хранения. В нем используются вызовы соответствующих методов класса CMatrix, которые выполняют требуемый алгоритм в зависимости от конкретного объекта-матрицы. Всего от класса CMatrix произведено шесть классов, реализующих рассмотренные схемы хранения и алгоритмы работы с разреженными матрицами.
Неполная факторизация в свою очередь производится виртуальным методом основного класса CIFSolver, выполняющего решение соответствующей СЛАУ. От этого класса произведены классы, выполняющие неполную факторизацию Холесского или неполную блочную факторизацию для рассмотренных блочно-трехдиагональных и блочно-пятидиагональных (в трехмерном случае) матриц.
Основная процедура решения, таким образом, получает на входе указатель на объект-матрицу системы и на объект-вектор свободных членов. Далее выполняется неполная факторизация матрицы и итерационный алгоритм решения с заданной извне точностью є. Параметры неполной факторизации задаются в конкретном объекте, реализующем соответствующий алгоритм.
Дополнительный виртуальный метод основного класса используется для сравнения результатов с задаваемым извне точным решением (эталоном). Если этот метод перекрыт в производном классе, считается, что точное решение известно и решение методом неполных факторизации производится для проверки функционирования алгоритма. В этом случае в выходной файл кроме результатов решения и количества итераций выводится максимальная абсолютная и относительная невязка.
Для решения краевых задач методами неполной факторизации предназначены отдельные классы CParabPDE и CEllipticPDE, обеспечивающие непосредственно формирование матрицы (объект класса, производного от CMatrix) и вектора свободных членов СЛАУ (объект класса, производного от CVector) в результате разностной аппроксимации исходной дифференциальной задачи с учетом вида используемых в задачах (1.2), (1.13) функций f(x,y), f(x,y,z), p(x,y,t), tp(x,y,z), а также с учетом вида области П.
Функции / и р вычисляются при помощи виртуальных методов классов CParabPDE и CEllipticPDE с соответствующими параметрами. Для задания вида области используется виртуальная функция nodeType{i,j) (или nodeType(i,j, к)), которая возвращает одно из трех значений: узел является внутренним, граничным или расположен вне области. После того, как матрица и вектор свободных членов СЛАУ сформированы, выполняется ее решение при помощи соответствующего объекта, производного от класса CIFSolver и реализующего конкретный тип неполной факторизации.
Все эти операции выполняет основной метод solve класса CParabPDE и CEllipticPDE, параметры которого задают величину шагов h и г (для параболических уравнений), точность вычислений, используемая при проверке условия останова итераций, тип неполной факторизации (блочная или Холесского), параметр к, задающий структуру допустимого заполнения в предобуславливателе.
Фрагмент исходного текста программы, содержащий описание перечисленных классов, приведен в приложении 2.
Для проверки полученных в работе алгоритмов численного решения параболических и эллиптических краевых задач методами неполной факторизации было поставлено два численных эксперимента на задачах, для которых известен точный аналитический вид решения, что позволило кроме практической проверки скорости сходимости метода вычислить также невязку по отношению к точному решению. Кроме того, вычислялось также число обусловленности матрицы М гА для экспериментальной проверки полученных в работе теоретических результатов (время вычисления числа обусловленности не принималось во внимание и не фиксировалось). Все вычисления проводились на компьютере с процессором типа Celeron с тактовой частотой 333 МГц и 64 МБайт оперативной памяти (24 свободно). В первом эксперименте решалась параболическая краевая задача с переменными коэффициентами (1.2), в которой соответствующие функции имели вид
