The comparison of high-speed methods efficiency for solvinga difference elliptical SLAE.pdf Разностная аппроксимация многомерных дифференциальных эллиптическихуравнений, описывающих задачи динамики вязкой жидкости и тепломассопере-носа приводит, как правило, к системам линейных алгебраических уравнений.При этом корректные способы аппроксимации дифференциальных операторовсохраняют свойство эллиптичности исходной краевой задачи на уровне разност-ных схем. Иными словами, необходимое условие однозначной разрешимости рас-сматриваемых задач: удовлетворение принципу максимума, - выполняется нетолько для дифференциальной постановки, но и для соответствующей ей разност-ной [1 - 3]. В подобных случаях полученные системы линейных алгебраическихуравнений, для подчеркивания преемственности свойств эллиптичности исходныхзадач, удобно называть разностными эллиптическими.Особенностью матриц таких СЛАУ является их высокий порядок и специфи-ческая разреженная структура. К тому же они могут быть плохо обусловленными.Для решения подобных систем существует большое количество хорошо зареко-мендовавших себя численных методов, и при этом постоянно разрабатываютсяновые, еще более эффективные алгоритмы [1, 3]. Понятно, что в связи с этим су-ществует проблема выбора наиболее подходящего метода для решения конкрет-ной задачи, который бы характеризовался высокой скоростью сходимости, вы-числительной устойчивостью и минимальными потребностями в компьютерныхресурсах. Подобный выбор, как правило, производится путем сравнения эффек-тивности различных методов при их использовании для решения одной и той жетестовой задачи. При этом необходимо заметить, что если эта задача будет отно-сительно простой [1, 4], то выбор нужного метода окажется затруднен, посколькувсе современные методы в подобной ситуации демонстрируют практически оди-наковую высокую эффективность: обычно для сходимости решения с разумнойстепенью точности требуется несколько десятков итераций. Следовательно, длясравнительного анализа характеристик современных методов желательно исполь-зовать тестовые задачи повышенной сложности. Такие задачи, в частности, сфор-мулированы в работе [5], дискретизация математических постановок которыхприводит к системам с плохо обусловленными матрицами.72 А.А. Фомин, Л.Н. ФоминаВ настоящей работе путем решения тестовых задач [5] сравниваются следую-щие высокоскоростные методы решения СЛАУ: стабилизированный метод би-сопряженных градиентов (Bi-CGStab) [5], он же с предобуславливателем на базенеполной LU факторизацией матрицы системы (Bi-CGStab P LU) [5], он же с пре-добуславливателем на базе явного метода Булеева (Bi-CGStab P B) [6]; модифици-рованный полинейный метод В. Г. Зверева (mLL) [4]; полинейный рекуррентныйметод с первым порядком аппроксимации (LR1) и со вторым порядком аппрокси-мации (LR2) [7, 8]. Здесь для удобства изложения материала в круглых скобкахприведены аббревиатуры названных методов.При решении всех задач разбиение расчетной области проводилось с равномер-ным шагом, а дискретизация уравнений выполнялась по методике Патанкара [9].Решение считалось найденным в случае достижения следующего условия сходимо-сти 02 2/ kR R < ε , где Rk, R0 - соответственно текущая и начальная невязки, а ƒ -заданная точность решения, которая во всех расчетах принималась равной 10-10.Для методов, в которых присутствует итерационный параметр, при решении каж-дой задачи он подбирался индивидуально путем численного эксперимента, цельюкоторого была минимизация количества итераций, необходимых для достижениязаданной точности. Таким образом, производилась своеобразная оценка сверхувозможностей каждого из методов.Задача 1. В качестве отправной точки исследования проведено сравнение вы-бранных методов на решении относительно простой задачи для эллиптическогоуравнения с переменными коэффициентами без смешанных производных( ) ( ) (, ) x yx x y yν u + ν u = S x y в единичном квадрате с краевыми условиями первогорода. Коэффициенты при старших производных и точное решение задачи известны:( ) ( ) 2 21 2 0,5 0,5 xx y ƒ = + ⎡ − + − ⎤ ⎣ ⎦ ,( ) ( ) 2 2y 1 2 0,5 0,5 0,5x y ƒ = + ⎡ − − − − ⎤ ⎣ ⎦ ,( ) [ ( ) ( )]2u x, y = 256 xy 1− x 1− y .Правая часть S(x,y) вычисляется путем подстановки xν , yν и u(x,y) в уравнение.Сеточное разбиение области составляет 101101=10 201 узлов.На рис. 1 представлено поведение отношения норм невязок в зависимости отномера итерации. Хорошо видно, что кривые 4,5,6 (стабилизированного методабисопряженных градиентов) ведут себя немонотонным образом, что, впрочем,было показано еще в работе [5]. При этом использование предобуславливателейповышает эффективность метода, причем предобуславливатель, построенный наоснове явного метода Булеева, позволил достигнуть наибольшего ускорения схо-димости. Данный результат полностью согласуется с выводами работы [6]. По-скольку разные методы для расчета одной итерации потребляют разное время, тоэффективность того или иного метода следует также оценивать с точки зрениясуммарного времени, необходимого для сходимости решения с заданной точно-стью. В данном случае эти времена составили, при прочих равных условиях, сле-дующие значения: LR2 - 0,342 с, LR1 - 0,373 с, mLL - 0,519 с, Bi-CGStab P B -0,715 с, Bi-CGStab P LU - 1,499 с, Bi-CGStab - 1,810 с. В целом, методы не вариа-ционного типа (кривые 1, 2, 3) показали себя как более скоростные по отношениюк стабилизированному методу бисопряженных градиентов.Сравнение эффективности высокоскоростных методов решения 731 10 100 k0-1-3-5-9-11-7ε = 10−10lg RkR0123456Рис. 1. 1 - LR2; 2 - LR1; 3 - mLL; 4 - Bi-CGStab P B;5 -Bi-CGStab P LU; 6 - Bi-CGStabЗадача 2. В [5] вторая тестовая задача формулируется следующим образом:Система линейных уравнений следует из пятиточечной конечно-разностной дис-кретизации уравнения в частных производных ( ) ( ) 1 x x y yDu Du − − = на единич-ном квадрате с условием Дирихле на границе y = 0 и условиями Неймана на ос-тальных границах. Сеточный шаг выбирался из условия, чтобы система линейныхуравнений имела 150150=22 500 неизвестных. Функция D определяется какD = 1000 для 0,1 ≤ x, y ≤ 0,9 и D = 1 в оставшейся части области.Для проведения расчетов необходимо было доопределить постановку кон-кретными граничными условиями. В данном случае принято:01 0 11, 1, 1, 1yy x xu u uuy x x== ==  = = = − =  .На рис. 2 показана динамика отношений норм невязок в зависимости от номе-ра итерации k для различных методов решения задачи. Видно, что за разумное ко-личество итераций (161) сходимость решения удалось достичь только методомLR1. Использование метода mLL достаточно быстро, за первую сотню итераций,выводит отношение норм невязок на асимптоту ≈ 6,310-3. Обе реализации методабисопряженных градиентов (Bi-CGStab P LU и Bi-CGStab P B) дают неустойчивоеповедение решения, в итоге приводящее к его расходимости. Любопытно поведе-ние отношения норм невязок при применении метода LR2: характеризуясь отно-сительно медленной сходимостью, оно, тем не менее, практически монотонно по-нижается и на приведенном графике за 500 итераций достигает величины ≈ 5,010-3.Дальнейший расчет, ограниченный 5000 итерациями, позволил достигнуть вели-чины отношения ≈ 3,610-8. Здесь следует отметить, что для достижения макси-мальной скорости сходимости в алгоритмах LR1 и LR2 в значении итерационногопараметра пришлось учес 5 и даже 6 значащих цифр. Отсюда можно предполо-жить, что использование в вычислениях арифметики чисел с двойной точностью74 А.А. Фомин, Л.Н. Фоминаоказалось недостаточным для эффективной реализацией алгоритма LR2 при ре-шении данной задачи, поскольку теоретически алгоритм LR2 должен быть болеебыстрым, чем LR1.1 10 100 klg RkR00-1-3-5-9-11-7ε = 10−101313452Рис. 2. 1 - LR2; 2 - LR1; 3 - mLL; 4 - Bi-CGStab P B;5 - Bi-CGStab P LUЗадача 3. В [5] третья тестовая задача формулируется следующим образом:Система линейных уравнений с несимметричной матрицей получена путемдискретизации краевой задачи [( ) ] / 2 1 xx yy x xu u au au − − + + = на единичномквадрате с граничными условиями Дирихле на всех границах. При этома = 20 exp[3,5 (x2+y2)]. Разбиение области произведено прямоугольной сеткой сшагом 1/201.Как и в предыдущей задаче, потребовалось доопределить искомую функциюна границах области, а именно: u 1Γ= . Необходимо заметить, что в данной по-становке матрица исходной СЛАУ не является симметричным оператором. Приэтом решение было построено всеми пятью методами. Динамика соответствую-щих норм невязок в зависимости от номера итерации k представлена на рис. 3.Методу LR1 (LR2) для достижения заданной точности потребовалось 10 итераций(0,911 с и 0,903 с соответственно), mLL - 19 итераций (1,164 с), Bi-CGStab P B -49 итераций (3,377 с), Bi-CGStab P LU - 116 итераций (7,843 с). Как и в первой за-даче, процесс сходимости решения при использовании метода Bi-CGStab носитявно немонотонный характер. Взаимное расположение кривых говорит о пре-имущественных скоростных характеристиках методов LR1 (LR2) и mLL передBi-CGStab с учетом, разумеется, оптимального подбора итерационных параметровдля каждого из методов.Сравнение эффективности высокоскоростных методов решения 751 10 100 k0-1-3-5-11ε = 10−10lg RkR012345Рис. 3. 1 - LR2; 2 - LR1; 3 - mLL; 4 - Bi-CGStab P B;5 - Bi-CGStab P LUЗадача 4. В [5] четвертая тестовая задача формулируется следующим образом:Система линейных уравнений с несимметричной матрицей получена путем дис-кретизации краевой задачи ( ) ( ) ( , ) x x y y xAu Au B x y u F − − + = на единичном квад-рате со следующими граничными условиями Дирихле 10yu== , на остальныхграницах u = 1. Коэффициенты уравнения определяются следующим образом:В(x,y) = 2 exp[2 (x2+y2)]; F = 100 в небольшой квадратной области в центре еди-ничного квадрата, в основной его части F = 0; определение А следует из схемыобласти. Расчетная область покрывается прямоугольной сеткой с шагом 1/128.На рис. 4 представлена схема области с точностью до ее графического пред-ставления в тексте статьи [5]. Поскольку численные размеры элементов схемыобласти в статье не указаны, их пришлось доопределять путем непосредственныхизмерений изображения.Динамика отношений норм невязок для различных методов приведена нарис. 5. Видно, что требуемой точности ¦ = 10-10 удалось достичь только методамиBi-CGStab P B и Bi-CGStab P LU. При этом количество итераций составило дляBi-CGStab P B - 63, для Bi-CGStab P LU - 119. В это же время методы LR1 (LR2) иmLL остановились в районе величины ƒ* 2,010-5, которая, что характерно,оказалась локально-критической и для метода Bi-CGStab P B и Bi-CGStab P LU:достигнув этого значения, отношение норм невязки несколько десятков итерацийоставалось практически постоянной (с ∼20 по ∼ 40 итерацию для Bi-CGStab P B ис ∼ 40 по ∼ 70 итерацию для Bi-CGStab P LU) и только потом, после небольшогоувеличения, понизилась до требуемой величины. Напрашивается естественныйвывод, что значения ƒ* есть функция постановки задачи (СЛАУ), а не применяе-мого для ее решения метода. Что же касается стабилизированного метода бисоп-ряженных градиентов без предобуславливателя Bi-CGStab, то он, очень медленноснижая отношение норм невязок, стремится также к величине ƒ* 2,010-5.76 А.А. Фомин, Л.Н. Фомина00,20,30,450,550,70,810,2 0,3 0,45 0,55 0,7 0,8 10,49 0,51F=102A = 104A = 10−5A = 102xyРис. 41 10 100 klg RkR00-1-3-5-9- 1 1-7ε = 10−101123454 56Рис. 5. 1 - LR2; 2 - LR1; 3 - mLL; 4 - Bi-CGStab P B;5 - Bi-CGStab P LU; 6 - Bi-CGStabСравнение эффективности высокоскоростных методов решения 77Анализ полученных результатов позволяет сделать следующие выводы.1. Для выявления предельных возможностей рассматриваемых методов в плантестовых расчетов необходимо включать задачи с плохо обусловленными матри-цами.2. Для задач с хорошими матрицами, которые характеризуются относитель-но небольшими числами обусловленности, методы [4,7] показали более высокуюэффективность по отношению к стабилизированному методу бисопряженныхградиентов [5].3. Наиболее эффективная реализация стабилизированного метода бисопря-женных градиентов Bi-CGStab P B использует итерационный параметр, что, всвою очередь, приводит к потере одного из основных преимуществ методов ва-риационного типа, которые в общем случае не требуют оптимизации итерацион-ного параметра.4. При решении задач с плохо обусловленными матрицами желательно иметь варсенале несколько принципиально различных методов, что, безусловно, повыша-ет вероятность успешного решения таких задач.

Скачать электронную версию публикации