CC BY
95
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2009 Математика и механика № 2(6)
УДК 519.632.4
А.А. Фомин, Л.Н. Фомина СРАВНЕНИЕ ЭФФЕКТИВНОСТИ ВЫСОКОСКОРОСТНЫХ МЕТОДОВ РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ СЛАУ
Проводится сравнительный анализ эффективности высокоскоростных методов решения пятидиагональных систем линейных алгебраических уравнений (СЛАУ), возникающих при разностной аппроксимации двумерных уравнений динамики вязкой жидкости и тепломассопереноса. Решаются четыре тестовые задачи. Анализируется эффективность использования рассматриваемых методов.
Ключевые слова: система линейных алгебраических уравнений, итерационный метод решения, тестовые задачи повышенной сложности.
Разностная аппроксимация многомерных дифференциальных эллиптических уравнений, описывающих задачи динамики вязкой жидкости и тепломассопереноса приводит, как правило, к системам линейных алгебраических уравнений. При этом корректные способы аппроксимации дифференциальных операторов сохраняют свойство эллиптичности исходной краевой задачи на уровне разностных схем. Иными словами, необходимое условие однозначной разрешимости рассматриваемых задач: удовлетворение принципу максимума, - выполняется не только для дифференциальной постановки, но и для соответствующей ей разностной [1 - 3]. В подобных случаях полученные системы линейных алгебраических уравнений, для подчеркивания преемственности свойств эллиптичности исходных задач, удобно называть разностными эллиптическими.
Особенностью матриц таких СЛАУ является их высокий порядок и специфическая разреженная структура. К тому же они могут быть плохо обусловленными. Для решения подобных систем существует большое количество хорошо зарекомендовавших себя численных методов, и при этом постоянно разрабатываются новые, еще более эффективные алгоритмы [1, 3]. Понятно, что в связи с этим существует проблема выбора наиболее подходящего метода для решения конкретной задачи, который бы характеризовался высокой скоростью сходимости, вычислительной устойчивостью и минимальными потребностями в компьютерных ресурсах. Подобный выбор, как правило, производится путем сравнения эффективности различных методов при их использовании для решения одной и той же тестовой задачи. При этом необходимо заметить, что если эта задача будет относительно простой [1, 4], то выбор нужного метода окажется затруднен, поскольку все современные методы в подобной ситуации демонстрируют практически одинаковую высокую эффективность: обычно для сходимости решения с разумной степенью точности требуется несколько десятков итераций. Следовательно, для сравнительного анализа характеристик современных методов желательно использовать тестовые задачи повышенной сложности. Такие задачи, в частности, сформулированы в работе [5], дискретизация математических постановок которых приводит к системам с плохо обусловленными матрицами.
В настоящей работе путем решения тестовых задач [5] сравниваются следующие высокоскоростные методы решения СЛАУ: стабилизированный метод би-сопряженных градиентов (Б1-СОБ(аЪ) [5], он же с предобуславливателем на базе неполной ЬЦ факторизацией матрицы системы (Б1-СОБ(аЬ Р ЬЦ) [5], он же с предобуславливателем на базе явного метода Булеева (Б1-СОБ(аЬ Р Б) [6]; модифицированный полинейный метод В. Г. Зверева (тЬЬ) [4]; полинейный рекуррентный метод с первым порядком аппроксимации (ЬЯ1) и со вторым порядком аппроксимации (ЬЯ2) [7, 8]. Здесь для удобства изложения материала в круглых скобках приведены аббревиатуры названных методов.
При решении всех задач разбиение расчетной области проводилось с равномерным шагом, а дискретизация уравнений выполнялась по методике Патанкара [9]. Решение считалось найденным в случае достижения следующего условия сходимости ||лк||2 /||Л0||2 < 6 , где Як, Я0 - соответственно текущая и начальная невязки, а е -
заданная точность решения, которая во всех расчетах принималась равной 10-10. Для методов, в которых присутствует итерационный параметр, при решении каждой задачи он подбирался индивидуально путем численного эксперимента, целью которого была минимизация количества итераций, необходимых для достижения заданной точности. Таким образом, производилась своеобразная «оценка сверху» возможностей каждого из методов.
Задача 1. В качестве отправной точки исследования проведено сравнение выбранных методов на решении относительно простой задачи для эллиптического уравнения с переменными коэффициентами без смешанных производных (Vх и х) х + (Vу и у) у = Б (х, у) в единичном квадрате с краевыми условиями первого рода. Коэффициенты при старших производных и точное решение задачи известны: Vх = 1 + 2[(х - 0,5)2 + (у - 0,5)2 ] ,
^ = 1 + 2 [0,5 - (х - 0,5 )2 - (у - 0,5 )2 ] ,
и (х, у) = 256 [ху (1 - х)(1 - у )]2 .
Правая часть Б(х,у) вычисляется путем подстановки Vх, vy и и(х,у) в уравнение. Сеточное разбиение области составляет 101x101=10 201 узлов.
На рис. 1 представлено поведение отношения норм невязок в зависимости от номера итерации. Хорошо видно, что кривые 4,5,6 (стабилизированного метода бисопряженных градиентов) ведут себя немонотонным образом, что, впрочем, было показано еще в работе [5]. При этом использование предобуславливателей повышает эффективность метода, причем предобуславливатель, построенный на основе явного метода Булеева, позволил достигнуть наибольшего ускорения сходимости. Данный результат полностью согласуется с выводами работы [6]. Поскольку разные методы для расчета одной итерации потребляют разное время, то эффективность того или иного метода следует также оценивать с точки зрения суммарного времени, необходимого для сходимости решения с заданной точностью. В данном случае эти времена составили, при прочих равных условиях, следующие значения: ЬЯ2 - 0,342 с, ЬЯ1 - 0,373 с, тЬЬ - 0,519 с, Б1-СОБ1аЬ Р Б -
0,715 с, Б1-СОБ1аЬ Р ЬЦ - 1,499 с, Б1-СОБ1аЬ - 1,810 с. В целом, методы не вариационного типа (кривые 1, 2, 3) показали себя как более скоростные по отношению к стабилизированному методу бисопряженных градиентов.
Рис. 1. 1 -ЬК2; 2 -ЬК1; 3 - шЬЬ; 4 - Ві-СОБґаЬ РВ; 5 -Ві-СОБіаЬ Р Ш; 6 - Ві-СОБіаЬ
Задача 2. В [5] вторая тестовая задача формулируется следующим образом: «Система линейных уравнений следует из пятиточечной конечно-разностной дискретизации уравнения в частных производных -(Вих)х - (Биу) = 1 на единичном квадрате с условием Дирихле на границе у = 0 и условиями Неймана на остальных границах. Сеточный шаг выбирался из условия, чтобы система линейных уравнений имела 150x150=22 500 неизвестных. Функция Б определяется как Б = 1000 для 0,1 < х, у < 0,9 и Б = 1 в оставшейся части области».
Для проведения расчетов необходимо было доопределить постановку конкретными граничными условиями. В данном случае принято:
и у=о
ди 1 ди = -1, ди
ду ^=1 ’ дх х=0 дх
= 1.
X=1
На рис. 2 показана динамика отношений норм невязок в зависимости от номера итерации к для различных методов решения задачи. Видно, что за разумное количество итераций (161) сходимость решения удалось достичь только методом ЬЯ1. Использование метода шЬЬ достаточно быстро, за первую сотню итераций, выводит отношение норм невязок на асимптоту ~ 6,3-10-3. Обе реализации метода бисопряженных градиентов (Бі-ЄОЗїаЬ Р Ьи и Бі-ЄОБїаЬ Р В) дают неустойчивое поведение решения, в итоге приводящее к его расходимости. Любопытно поведение отношения норм невязок при применении метода ЬЯ2: характеризуясь относительно медленной сходимостью, оно, тем не менее, практически монотонно понижается и на приведенном графике за 500 итераций достигает величины ~ 5,0-10-3. Дальнейший расчет, ограниченный 5000 итерациями, позволил достигнуть величины отношения ~ 3,6-10-8. Здесь следует отметить, что для достижения максимальной скорости сходимости в алгоритмах ЬЯ1 и ЬЯ2 в значении итерационного параметра пришлось учесть 5 и даже 6 значащих цифр. Отсюда можно предположить, что использование в вычислениях арифметики чисел с двойной точностью
оказалось недостаточным для эффективной реализацией алгоритма ЬЯ2 при решении данной задачи, поскольку теоретически алгоритм ЬЯ2 должен быть более быстрым, чем ЬЯ1.
Рис. 2. 1 - ЬК2; 2 - ЬК\; 3 - шЬЬ; 4 - В1-С08гаЪ Р В;
5 - Вг-СОЪгаЪ Р Ш
Задача 3. В [5] третья тестовая задача формулируется следующим образом: «Система линейных уравнений с несимметричной матрицей получена путем дискретизации краевой задачи -ыхх - и№ + [(аи)х + аих]/2 = 1 на единичном
квадрате с граничными условиями Дирихле на всех границах. При этом а = 20 ехр[3,5 (х2+у2)]. Разбиение области произведено прямоугольной сеткой с шагом 1/201».
Как и в предыдущей задаче, потребовалось доопределить искомую функцию на границах области, а именно: и|г = 1. Необходимо заметить, что в данной постановке матрица исходной СЛАУ не является симметричным оператором. При этом решение было построено всеми пятью методами. Динамика соответствующих норм невязок в зависимости от номера итерации к представлена на рис. 3. Методу ЬЯ1 (ЬЯ2) для достижения заданной точности потребовалось 10 итераций (0,911 с и 0,903 с соответственно), тЬЬ - 19 итераций (1,164 с), Б1-СОБ1аЬ Р Б -49 итераций (3,377 с), Б1-СОБ1аЬ РЬи- 116 итераций (7,843 с). Как и в первой задаче, процесс сходимости решения при использовании метода Б1-СОБ1аЬ носит явно немонотонный характер. Взаимное расположение кривых говорит о преимущественных «скоростных» характеристиках методов ЬЯ1 (ЬЯ2) и тЬЬ перед Б1-СОБ1аЬ с учетом, разумеется, оптимального подбора итерационных параметров для каждого из методов.
Рис. 3. 1 -ЬК2; 2 -ЬК\; 3 - шЬЬ; 4 - В1-С08гаЪ РВ; 5 - Вг-СОЪгаЪ Р Ш
Задача 4. В [5] четвертая тестовая задача формулируется следующим образом: «Система линейных уравнений с несимметричной матрицей получена путем дискретизации краевой задачи -(Аих)х - (Аиу) + В(х, у) их = ^ на единичном квадрате со следующими граничными условиями Дирихле и|у=1 = 0, на остальных
границах и = 1. Коэффициенты уравнения определяются следующим образом: В(х,у) = 2 ехр[2 (х2+у2)]; Г = 100 в небольшой квадратной области в центре единичного квадрата, в основной его части Г = 0; определение А следует из схемы области. Расчетная область покрывается прямоугольной сеткой с шагом 1/128».
На рис. 4 представлена схема области с точностью до ее графического представления в тексте статьи [5]. Поскольку численные размеры элементов схемы области в статье не указаны, их пришлось доопределять путем непосредственных измерений изображения.
Динамика отношений норм невязок для различных методов приведена на рис. 5. Видно, что требуемой точности е = 10-10 удалось достичь только методами В1-С081аЪ Р В и В1-С081аЪ Р Ьи. При этом количество итераций составило для В1-СОБ(аЪ Р В - 63, для В1-СОБ1аЪ Р Ьи - 119. В это же время методы ЬЯ1 (ЬЯ2) и тЬЬ «остановились» в районе величины 5* « 2,0-10-5, которая, что характерно, оказалась локально-критической и для метода В1-СОБ1аЪ Р В и В1-СОБ1аЪ Р Ьи: достигнув этого значения, отношение норм невязки несколько десятков итераций оставалось практически постоянной (с ~20 по ~ 40 итерацию для В1-С081аЪ Р В и с ~ 40 по ~ 70 итерацию для В1-С081аЪ Р Ьи) и только потом, после небольшого увеличения, понизилась до требуемой величины. Напрашивается естественный вывод, что значения 5* есть функция постановки задачи (СЛАУ), а не применяемого для ее решения метода. Что же касается стабилизированного метода бисоп-ряженных градиентов без предобуславливателя В1-СОБ1аЪ, то он, очень медленно снижая отношение норм невязок, стремится также к величине 5*«2,0-10-5.
Рис. 4
Рис. 5. 1 - ЬК2; 2 - ЬК\; 3 - шЬЬ; 4 - В1-С08гаЪ РВ; 5 - В1-С08гаЪ Р Ьи; 6 - В1-С08гаЪ
Анализ полученных результатов позволяет сделать следующие выводы.
1. Для выявления предельных возможностей рассматриваемых методов в план тестовых расчетов необходимо включать задачи с плохо обусловленными матрицами.
2. Для задач с «хорошими» матрицами, которые характеризуются относительно небольшими числами обусловленности, методы [4,7] показали более высокую эффективность по отношению к стабилизированному методу бисопряженных градиентов [5].
3. Наиболее эффективная реализация стабилизированного метода бисопряженных градиентов Bi-CGStab P B использует итерационный параметр, что, в свою очередь, приводит к потере одного из основных преимуществ методов вариационного типа, которые в общем случае не требуют оптимизации итерационного параметра.
4. При решении задач с плохо обусловленными матрицами желательно иметь в арсенале несколько принципиально различных методов, что, безусловно, повышает вероятность успешного решения таких задач.
ЛИТЕРАТУРА
1. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 590 с.
2. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. М.: Эдиториал УРСС, 1999. 248 с.
3. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995. 288 с.
4. Зверев В.Г. Модифицированный полинейный метод решения разностных эллиптических уравнений // ЖВМ и МФ. 1998. Т. 38. № 9. С. 1553 - 1562.
5. Van der Vorst H.A. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. V. 13. No. 2. P. 631 - 644.
6. Старченко А.В. Сравнительный анализ некоторых итерационных методов для численного решения пространственной краевой задачи для уравнений эллиптического типа // Вестник ТГУ. Бюллетень оперативной научной информации. Томск: ТГУ, 2003. № 10. C. 70 - 80.
7. Фомин А.А., Фомина Л.Н. Полинейный рекуррентный метод решения разностных эллиптических уравнений. Кемерово: КемГУ, 2007. 78 с. Деп. в ВИНИТИ 06.04.2007, № 385-В2007.
8. Фомин А.А., Фомина Л.Н. Полинейный рекуррентный метод решения СЛАУ с пятидиагональной матрицей // Четвертая Сибирская школа-семинар по параллельным и высокопроизводительным вычислениям (Томск, 9 - 11 октября 2007 г.). Томск: Дельтаплан, 2008. C. 192 - 201.
9. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.
СВЕДЕНИЯ ОБ АВТОРАХ:
ФОМИН Александр Аркадьевич - кандидат физико-математических наук, руководитель отдела информационных технологий ОАО «Издательско-полиграфическое предприятие “Кузбасс”», Россия (г. Кемерово). E-mail: fomin_aa@mail.ru.
ФОМИНА Любовь Николаевна - старший преподаватель кафедры вычислительной математики Кемеровского государственного университета. E-mail: lubafomina@mail.ru
Статья принята в печать 26.03.2009 г.
