Многосеточные структурно-алгебраические алгоритмы

Министерство образования и науки Российской Федерации 
Сибирский федеральный университет 

Федеральное агентство научных организаций 
Институт вычислительного моделирования СО РАН

В. В. Ефремов
В. В. Шайдуров
Л. В. Гилева

МНОГОСЕТОЧНЫЕ

СТРУКТУРНО-АЛГЕБРАИЧЕСКИЕ

АЛГОРИТМЫ

Монография

Красноярск
СФУ
2016

УДК 519.712.2+004.421.2 

ББК 22.127+22.18

Е924

Р е ц е н з е н т ы:
Г. В. Муратова, доктор физико-математических наук, профессор 
кафедры 
высокопроизводительных 
вычислений 
и 
информационнокоммуникационных технологий Института математики, механики и компьютерных наук ФГАОУ ВПО "Южный федеральный университет";

А. В. Лапин, доктор физико-математических наук, профессор кафедры

математической статистики Института вычислительной математики ФГАОУ
ВПО "Казанский (Приволжский) федеральный университет"

Ефремов, В. В.

Е924

Многосеточные структурно-алгебраические алгоритмы : монография / 
В. В. Ефремов, В. В. Шайдуров, Л. В. Гилева. – Красноярск : 
Сиб. федер. ун-т, 2016. – 154 с.

ISBN 978-5-7638-3575-5

Представлены результаты исследований в области создания эффективных вычислительных алгоритмов для решения задач математической физики
многосеточными методами. Теоретическое обоснование подкреплено численными расчетами.
Предназначена для научных работников, преподавателей, студентов
старших курсов, магистрантов и аспирантов вузов, занимающихся численным
решением задач математической физики.

УДК 519.712.2+004.421.2

ББК 22.127+22.18

ISBN 978-5-7638-3575-5
c⃝ Сибирский федеральный университет, 2016

Содержание

Предисловие . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Глава 1. Элементы теории многосеточных алгоритмов
и используемые сеточные аппроксимации . . . . . . . . . . . . . . . 16
1.1. Общее описание многосеточных алгоритмов . . . . . . . . . . . . . . . . . . . .16
1.2. Некоторые свойства итерационных процедур . . . . . . . . . . . . . . . . . . .21
1.3. Дискретизация уравнения диффузии . . . . . . . . . . . . . . . . . . . . . . . . . . .25
1.4. Дискретизация уравнения теплопроводности . . . . . . . . . . . . . . . . . . . 35
1.5. Учет краевых условий на криволинейной границе . . . . . . . . . . . . . .40
1.6. Аппроксимация на прямоугольных адаптивных сетках . . . . . . . . 44

Глава 2. Обоснование сходимости двух вариантов
алгебраического многосеточного алгоритма . . . . . . . . . . . . .52
2.1. Поточечная и матричная формулировки основного
многосеточного алгоритма . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .52
2.2. Общая оценка сходимости многосеточного алгоритма . . . . . . . . . . 61
2.3. Оценка сходимости продольно-поперечной редукции . . . . . . . . . . . 70
2.4. Оценка сходимости квадратно-гнездовой редукции . . . . . . . . . . . . .83

Глава 3. Обоснование сходимости и вычислительные
эксперименты для некоторых задач . . . . . . . . . . . . . . . . . . . . . .94
3.1. Сходимость алгоритма на прямоугольной
составной сетке . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.2. Сходимость алгоритма для аппроксимации
уравнения теплопроводности . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.3. Сходимость в области с криволинейной границей . . . . . . . . . . . . . 107
3.4. Формулировка алгоритма и обоснование его сходимости
для задачи с переменными коэффициентами . . . . . . . . . . . . . . . . . . 117
3.5. Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Глава 4. Многосеточный метод с шахматным
исключением неизвестных . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.1. Шахматное исключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.2. Анализ Фурье . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .130
4.3. Вычислительный эксперимент для двухсеточного
варианта . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.4. Оценка трудоемкости многосеточного метода . . . . . . . . . . . . . . . . . 138
Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .144
Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

Предисловие

В первой главе настоящей монографии изложены основные сведения
о структуре многосеточных итерационных алгоритмов и системах линейных алгебраических уравнений, получаемых при дискретизации эллиптических уравнений методом конечных элементов.
В параграфе 1.1 дано описание структуры и основанной на ней классификации многосеточных итерационных алгоритмов на абстрактном
уровне. В параграфе 1.2 изложены вспомогательные сведения по ускорению сходимости итерационных процессов с применением свойств многочленов Чебышева, а также по спектральным свойствам энергетически
эквивалентных операторов.
В параграфе 1.3 изложен вывод элементарных матриц жесткости метода конечных элементов для билинейных элементов на квадратных ячейках и линейных элементов на треугольниках с последующим применением
простых квадратурных формул. Спектральная эквивалентность сеточных
операторов в последующем базируется на свойствах элементарных матриц
жесткости, поэтому изложение их конструкции в каждом случае проводится довольно подробно.
В параграфе 1.4 рассмотрены виды дискретизации по времени двумерного уравнения теплопроводности, сводящие исходную задачу к последовательности эллиптических уравнений, к которым затем применимы подходы, данные в предыдущем параграфе, для дискретизации по пространству.
В параграфе 1.5 изучена ситуация с аппроксимацией краевых условий
Дирихле для задачи с криволинейной границей. После дискретизации билинейными конечными элементами на регулярной квадратной сетке, покрывающей область решения, полученные узловые значения вне области
исключаются с помощью линейной аппроксимации значений в гранич

Предисловие

ных и внутренних узлах. Здесь также главным результатом являются вид
и свойства получаемых элементарных матриц жесткости приграничных
элементов.
В параграфе 1.6 рассмотрено применение метода конечных элементов при адаптации квадратной сетки к неоднородной гладкости решения
путем деления части квадратных ячеек сетки на четыре равных части.
Получающиеся на такой составной сетке элементарные матрицы жесткости имеют специфический вид, который в последующем используется для
обоснования сходимости.
В главе 2 приведено обоснование сходимости двух вариантов алгебраического многосеточного алгоритма и сопоставлены их скорости сходимости. В параграфе 2.1 сначала детально изложена поточечная формулировка основного многосеточного алгоритма, используемого далее в различных
модификациях, и дана его иллюстрация на примере уравнения Пуассона,
дискретизированного на равномерной квадратной сетке. Затем изложение
идет в более общей векторно-матричной форме.
Именно эта форма используется в последующих параграфах гл. 2 и 3
для формулировки многосеточных алгоритмов для разных приложений.
В том числе в параграфе 2.2 она использована на абстрактном уровне
для обоснования скорости сходимости двухсеточного варианта основного
алгоритма. Центральным пунктом обоснования является анализ и оценка собственных чисел операторов перехода нескольких матричных преобразований в алгоритме. Способ обоснования сходимости в этом разделе
довольно близок к идеям и результатам работ [59, 62], посвященным аналогичной тематике.
В параграфе 2.3 проведено полное и детальное обоснование сходимости многосеточного алгоритма для продольно-поперечной редукции неизвестных – простейшего варианта реализации. В нем сначала приближенно
прореживается каждый второй слой в направлении x, а затем – каждый
второй в направлении y. Этот двухступенчатый прием повторяется до тех
пор, пока не останется небольшое число неизвестных, для которых полученная система линейных алгебраических уравнений решается прямым
методом. Оценка числа арифметических операций для решения сеточного
уравнений Пуассона в полном многосеточном варианте дала почти опти

Предисловие
7

мальное количество операций O(N log2 N), где N – число неизвестных
исходной системы сеточных уравнений.
В параграфе 2.4 дано обоснование другой модификации многосеточного алгоритма, названной здесь квадратно-гнездовой. Геометрическая суть
редукции неизвестных в этой модификации состоит сначала в приближенном исключении центрального неизвестного в ячейке 3 × 3, а затем
в точном исключении средних неизвестных вдоль сторон получившегося квадрата. Это двухступенчатое прореживание уменьшает количество
неизвестных примерно в четыре раза и повторяется до остатка небольшого числа неизвестных и уравнений, которые решаются прямым методом.
На основе рекуррентного использования двухсеточных оценок параграфа 2.2 получена оценка сходимости полного многосеточного алгоритма и
подсчитано число необходимых арифметических операций. На этот раз
его верхняя оценка оказалась оптимальной: O(N).
В главе 3 обоснована сходимость, вычислены оценки скорости сходимости и вычислительной сложности многосеточного алгоритма для нескольких сеточных задач.
В параграфе 3.1 проведено обоснование сходимости для уравнения
Пуассона на прямоугольной составной сетке, полученной адаптацией части ячеек путем их деления на четыре равных части. В качестве основного
алгоритма в этой главе выбрана квадратно-гнездовая редукция. В итоге
получена оптимальная оценка скорости сходимости, как и в случае равномерной сетки.
В параграфе 3.2 исследуется многосеточный алгоритм для уравнения
теплопроводности на каждом шаге по времени. Полученная оценка скорости сходимости алгоритма делает неявные схемы конкурентоспособными
с методами расщепления и явными аппроксимациями по времени с точки
зрения экономичности решения. Но неявные схемы обладают некоторыми
преимуществами по аппроксимационным свойствам: отсутствие ограничений на шаг по времени, достижение меньшей схемной вязкости и более
высокого порядка аппроксимации по времени (как в схеме Кранка – Николсон).
В параграфе 3.3 алгоритм квадратно-гнездовой редукции применяется к сеточной задаче, аппроксимирующей уравнение Пуассона в области

Предисловие

с криволинейной границей. Проведен детальный анализ возможных ситуаций пересечения границы с равномерной прямоугольной сеткой и соответствующих свойств элементарных матриц жесткости. В итоге обоснована
сходимость многосеточного алгоритма, причем найденная оценка сходимости не уступает случаю прямоугольной области.

В параграфе 3.4 рассмотрена задача с переменными коэффициентами.
Здесь сначала применяется "внешний" итерационный процесс, основанный на спектральной эквивалентности сеточных операторов с переменными и постоянными коэффициентами. А внутри него на каждой итерации
применяется многосеточный алгоритм уже для задачи с постоянными коэффициентами. В современной терминологии это означает использование
многосеточного алгоритма как предобусловливателя во внешнем итерационном процессе. Доказано, что такая комбинация остается оптимальной
по числу арифметических операций.

В параграфе 3.5 изложены результаты вычислительных экспериментов. Здесь численно сопоставлены модификации многосеточного алгоритма, их сходимость, применимость для задач с переменными коэффициентами.

В главе 4 рассмотрен довольно распространенный алгебраический многосеточный метод с "шахматным" исключением неизвестных. В параграфе 4.1 описана общая схема алгоритма, при которой приближенная операция прореживания неизвестных вдвое проводится на каждом шаге редукции для неизвестных, чередующихся в шахматном порядке на равномерной прямоугольной сетке. В зарубежной литературе такое чередование
часто называется красно-черным. Отметим, что в этом алгоритме используется только проектирование на более редкую сетку без сглаживающих
итераций.

В параграфе 4.2 для исследования сходимости алгоритма и определения ее скорости применяется прямой анализ Фурье операторов перехода
в отличие от двух предыдущих глав, где использовалась спектральная
эквивалентность операторов задачи. Получен явный вид функции подавления гармоник Фурье в погрешности приближенного решения во всей
спектральной области для стандартного оператора проектирования Mh.

Предисловие
9

Предложена другая модификация оператора проектирования ˜
Mh, обеспечивающая вчетверо меньший коэффициент уменьшения погрешности.
В параграфе 4.3 с помощью вычислительных экспериментов продемонстрированы свойства операторов Mh и
˜
Mh, подтвердившие эффективность предложенного алгоритма с модифицированным оператором проектирования.
В параграфе 4.4 разбираются возможные тактики использования построенных в этой главе многосеточных алгоритмов. Для стандартного
оператора Mh представлена новая тактика, позволяющая достигнуть оптимального показателя в числе арифметических операций O(N) в отличие от стандартной рекомендуемой тактики с почти оптимальным числом
операций O(N log2 N). А для нового оператора ˜
Mh предложена тактика
с оптимальным показателем оценки cN числа арифметических операций
и константой c, в несколько раз меньшей, чем для оператора Mh.
Мы приносим благодарность Российскому научному фонду за поддержку части выполненной многолетней работы в рамках проекта № 1411-00147.

Введение

Многосеточный метод стал одним из эффективных итерационных методов решения систем уравнений, получаемых при дискретизации дифференциальных задач. Его алгоритмическую основу составляют два приема,
которые поясним на примере сеточной задачи, соответствующей довольно
мелкой сетке и имеющей большое число неизвестных. Первый из них (называемый далее восходящим) заключается в последовательном решении
сеточных задач, аналогичных исходной, но соответствующих более крупным сеткам. Решение начинается с самой грубой сетки, где оно может быть
осуществлено довольно экономично, например прямым методом. Затем
полученное решение интерполируется на более мелкую сетку и используется в качестве начального приближения в каком-либо итерационном
процессе на этом более высоком уровне. При одинаковом порядке размера шагов двух сеток такое начальное приближение уже имеет точность,
близкую к удовлетворительной. В итоге в итерационном процессе требуется существенно меньше итераций для доведения точности до требуемого
уровня. Этот прием был довольно распространен во времена ручного и
механизированного счета, когда при вынужденном переходе на более мелкую сетку пытались использовать всю дорого доставшуюся информацию
с грубой сетки.
Второй прием, который будем называть нисходящим, предложен в работе Р.П. Федоренко [25]. Он основан на быстрой сходимости некоторых
итерационных процессов для высокочастотных гармоник, вклад которых
в погрешность существенно подавляется за несколько итераций. Низкочастотные гармоники в погрешности подавляются гораздо медленнее и поэтому будут составлять ее основную часть. В результате погрешность можно рассматривать как некоторую плавно меняющуюся сеточную функцию, для которой можно выписать систему сеточных уравнений с невяз

Введение
11

кой в правой части. В принципе, эту систему можно решить и найти погрешность, но такой путь довольно трудоемок. Вместо этого полученной
системе поставим в соответствие дифференциальную задачу с достаточно плавным гладким решением, для которой снова построим сеточную
систему на сетке с более крупным шагом (например, в два или три раза
крупнее), т.е. на более низком уровне. Такую систему можно решить более
экономично, поскольку число неизвестных, например в двумерном случае,
уменьшится в четыре или девять раз. Ввиду гладкости ее решение будет
хорошо приближать искомую погрешность на мелкой сетке.
Несмотря на относительную алгоритмическую сложность многосеточного метода, в 1964 году Р.П. Федоренко [26] обосновал сходимость многосеточного метода для конечно-разностного аналога уравнения Пуассона
в квадрате, а в 1966 году Н.С. Бахвалов [4] доказал оптимальность метода по числу арифметических операций для достижения точности, согласованной с порядком сходимости. Оптимальность для дискретной задачи
с N неизвестными характеризуется числом O(N) арифметических действий для нахождения приближенного решения.
В итоге по асимптотическим оценкам эффективности многосеточный
метод опередил известные алгоритмы, но логическая сложность и громоздкое математическое обоснование на некоторое время завуалировали достоинства этого метода. На определенном этапе развития метода конечных элементов привлечение нового математического аппарата
и программных реализаций существенно снизило трудоемкость алгоритма и упростило его обоснование. Поэтому в конце 1970-х гг. начался рост
числа публикаций по многосеточным методам (в основном за рубежом),
продолжающийся в настоящее время.
По способу обоснования сходимости многосеточные методы можно разделить на три непересекающихся класса. Один из них объединяет так
называемые каскадные алгоритмы [74], в которых используется только
восходящий прием. Основу их обоснования составляет существенная близость решений сеточных задач на соседних уровнях, обычно вытекающая
из значительной гладкости решения дифференциальной задачи.
Второй класс, часто называемый алгебраическим, по-существу использует только нисходящий прием. Большая группа таких методов относится

Введение

к алгоритмам типа "черный ящик". Этот термин объединяет алгоритмы,
использующие минимальную информацию о свойствах и происхождении
дискретной задачи [52, 53]. Несмотря на вычислительные успехи таких методов, теоретическое обоснование их сходимости достигнуто лишь в редких случаях.

Другая группа методов использует структурную геометрическую или
топологическую информацию о способе формирования матрицы системы уравнений или расположении ее ненулевых элементов [32, 34, 35, 38,
59, 62]. В настоящей монографии мы будем активно привлекать поэлементную форму сборки глобальной матрицы жесткости в методе конечных элементов и другие общие принципы формирования матрицы в методе конечных разностей. Обоснование излагаемых алгоритмов опирается
на энергетическую эквивалентность (схожесть спектров) оператора сеточной задачи и некоторого оператора более простой структуры, который затем сводится к оператору более низкого уровня путем исключения части
неизвестных [62, 75]. Обратная процедура восстановления исключенных
неизвестных хотя и напоминает формулы интерполяции сеточных значений с нижнего уровня на более высокий, но таковой не является. Отметим, что при доказательстве сходимости обычно не налагаются условия
на гладкость дифференциального решения, что определяет несомненное
достоинство такого обоснования. Вместе с тем формирование эквивалентного оператора специальной, хотя и более простой, структуры – довольно
непростая задача. Поэтому круг дифференциальных задач, решаемых алгебраическими многосеточными методами, на текущий момент не очень
широк.

Третий класс будем называть классическим, поскольку он появился
первым и проработан наиболее детально [4, 25, 26, 29, 58]. Он использует и восходящий, и нисходящий прием. В сравнении с каскадным алгоритмом он менее требователен к гладкости дифференциального решения
и обладает большей областью применения. Но алгоритмическая структура
классического метода сложнее, чем у каскадного. А в сравнении с алгебраическими алгоритмами они более требовательны к гладкости дифференциального решения. Однако область их применения значительно ши