О повышении эффективности полинейного рекуррентного метода решения разностных эллиптических уравнений

скачать О ПОВЫШЕНИИ ЭФФЕКТИВНОСТИ ПОЛИНЕЙНОГО РЕКУРРЕНТНОГО МЕТОДА РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ Л.Н. Фомина Кемеровский государственный университет, Кемерово ABOUT IMPROVEMENT OF THE EFFECTIVENESS OF LINE-BY-LINE RECURSIVE METHOD FOR SOLVING A DIFFERENCE ELLIPTICAL EQUATIONS. L.N. Fomina Kemerovo state university, Kemerovo The article regards the technology to improve the effectiveness of line-by-line recursive method for solving a SLAE with the five-diagonal positive type matrixes by a combination of line-by-line recursive method and the biconjugate gradient stabilized method. It is shown that the traditional way by building a specialized preconditioner based on line-by-line recursive method raises both resolving possibilities of the method and time which is needed to find a SLAE solution. A new method allowing not only to solve much more rigid SLAE, but also to spend much less time to solve them is successfully developed by so-called direct combining of algorithms of line-by-line recursive method and the biconjugate gradient stabilized method. There is no need to use the preconditioner. The possibilities of this method are demonstrated by means of computing experiment. The superiority of new method both over initial line-by-line recursive method and over the most effective at present time biconjugate gradient stabilized method with preconditioner based on explicit Buleev method is shown. Введение. Современные методы численного решения уравнений математической физики, как правило, требуют их разностной дискретизации, в результате чего возникают системы линейных алгебраических уравнений (СЛАУ) большой размерности, матрицы которых имеют разреженно-упорядоченную структуру [1-3]. Матрицы возникающих при этом систем могут быть плохо обусловлены не только в силу особенностей физических постановок задач, но и за счет роста количества переменных при использовании более подробных сеточных покрытий областей определения задач. Поэтому возможности решения более сложных задач математической физики опираются не только на дальнейший рост мощностей ЭВМ, но и на методы, обладающими повышенной эффективностью решения СЛАУ. В последнее время широкое распространение получили градиентные методы, использующие процедуру ортогонализации или биортогонализации векторов невязок [1,2,4], характеризующиеся достаточно высокими разрешающими возможностями при решении жестких систем. Однако, что касается эффективности этих методов с точки зрения быстродействия, то положение здесь оставляет желать много лучшего. Как известно [1,2], одним из выходов из создавшейся ситуации является использование предобуславливателей, построенных на базе того или иного релаксационного метода решения СЛАУ. Причем эффективность такого комбинированного метода напрямую зависит от эффективности используемого в нем релаксационного метода [5,6]. ^ Постановка задачи. Пусть имеет место двумерная по пространству краевая задача в единичной области . Внутри области поведение искомой функции Ф(x,y) описывается дифференциальным уравнением (1) где ^ U, V – компоненты скорости, – коэффициенты диффузии, S – источник. На линиях границы области Г в общем случае имеют место условия третьего рода , (2) где aГ, bГ, cГ – известные величины, а n – нормаль к границе вдоль координатной линии. Здесь U, V, – непрерывные дифференцируемые функции внутри и на границе , причем . Выбирая произвольные выражения для функций Ф, U, V, , а также произвольные коэффициенты аГ, bГ, можно всегда, путем вычисления источника S из (1) и коэффициента сГ из (2), сформулировать необходимую для решения тестовую задачу. В настоящей работе системы линейных алгебраических уравнений вида AФ = b получаются путем разностной аппроксимации задачи (1) - (2) методом контрольного объема (вариант экспоненциальной схемы) [3] на равномерной сетке, покрывающей расчетную область . Итерационный процесс построения решения полученных СЛАУ прекращается, когда отношение норм невязок достигает заданной точности , где Rk – невязка на k итерации, R0 – начальная невязка,  - точность решения. Все расчеты проводились на вычислительной системе Intel Core i5 750 2.66GHz, RAM 4Gb, Win32, Fortran. ^ Предварительный анализ. В статье [4] приводится алгоритм метода бисопряженных градиентов со стабилизацией и на примере решения нескольких модельных задач повышенной сложности, дискретизация которых приводит к системам с плохо обусловленными матрицами, демонстрируются его возможности. Одна из этих задач формулируется следующим образом: «Система линейных уравнений с несимметричной матрицей получена путем дискретизации краевой задачи на единичном квадрате со следующими граничными условиями Дирихле , на остальных границах u = 1. Коэффициенты уравнения определяются следующим образом: В(x,y) = 2 exp[2 (x2+y2)]; F = 100 в небольшой квадратной области в центре единичного квадрата, в основной его части F = 0; определение А следует из схемы области. Расчетная область покрывается прямоугольной сеткой с шагом 1/128». При этом оценка числа обусловленности для сетки 128×128 составляет МА  3,26×108. На рис. 1 показан процесс сходимости при решении данной модельной задачи различными методами: явным методом Н.И. Булеева, бисопряженным градиентным методом со стабилизацией без предобуславливателя (Bi-CGStab P) [4], полинейным рекуррентным методом первого порядка (LR1) [6-8], бисопряженным градиентным методом со стабилизацией с предобуславливателем, построенным на основе явного метода Н.И. Булеева (Bi-CGStab P B) [5]. Рис. 1 – Поведение отношений норм невязок Rk/R 0 в зависимости от номера итерации k. 1 – явный метод Н.И. Булеева; 2 – Bi-CGStab P; 3 – LR1; 4 – Bi-CGStab P B. Видно, что метод Bi-CGStab P, показывая свое характерное не монотонное поведение, за тысячу итераций не смог достичь точности даже 10-3. Релаксационные же методы (явный метод Н.И. Булеева и LR1) «остановились» в районе величины  2×10-5, причем эта величина является локально-критической (с 18 по 40 итерацию) и для метода Bi-CGStab P B, но он смог достичь заданной точности e = 10-8 за 57 итераций. Из поведения графиков на рис. 1 можно сделать следующий вывод: по отдельности ни явный метод Н.И. Булеева, ни Bi-CGStab P не в состоянии решить данную задачу с требуемой точностью, однако их комбинация успешно позволяет сделать это. В этой ситуации естественно предположить, что комбинация бисопряженного градиентного и полинейного рекуррентного методов также позволит создать более эффективный метод. ^ Градиентный метод решения СЛАУ. Алгоритм метода Bi-CGStab P для случая использования предобуславливателя имеет следующий вид [4]. 1. Ф0 –вектор начального приближения решения; 2. r0 = b – A Ф0; 3. – произвольный вектор, такой, что , например, ; 4. ; 5. ; 6. для k = 1,2,3, … 7. ; 8. ; 9. определение вектора y из решения системы ; 10. ; 11. ; 12. ; 13. определение вектора z из решения системы ; 14. ; 15. ; 16. ; 17. если Фk достигло требуемой точности – выход из цикла; 18. ; 19. конец цикла по k Здесь K – предобуславливатель, причем предполагается разложение K = K1K2, где K1 – нижняя, а K2 – верхняя треугольные матрицы. При наличии такого разложения строки 9 и 13 алгоритма Bi-CGStab P выполняются с линейной трудоемкостью относительно числа неизвестных. Следуя общей идеологии комбинации градиентного и релаксационного методов был создан предобуславливатель K = K1K2 путем неполной факторизации матрицы системы А на сомножители K1 и K2 таким образом, чтобы в основе этого разложения лежал алгоритм LR1. В целом такой подход не привел к желаемому результату, поскольку скорость сходимости нового метода заметно понизилась по отношению к исходному методу LR1. ^ Прямое сочетание алгоритмов LR1 и Bi-CGStab P. При детальном изучении особенностей алгоритма Bi-CGStab P можно заметить, что в нем используется только правое предобуславливание. При этом не применяется в явном виде разложение K на сомножители K1 и K2, но в неявном виде это разложение, естественно, предполагается для того, чтобы иметь возможность решить системы (строка 9 алгоритма Bi-CGStab P) и (строка 13 алгоритма Bi-CGStab P) с линейной трудоемкостью относительно числа неизвестных. Действительно, путем последовательного решения двух соответствующих пар систем и это легко сделать в силу треугольных структур матриц K1 и K2. По определению предобуславливатель – это легко обратимая матрица K в некотором смысле близкая к исходной, то есть . При этом считается, что качество предобуславливателя тем выше, чем ближе он к исходной матрице системы. Вполне естественно эту близость рассматривать как близость решений x и xA двух систем и . С этой точки зрения вектор y, который является точным решением системы (строка 9 алгоритма Bi-CGStab P), можно трактовать как приближенное решение системы . Соответственно вектор z (строка 13 алгоритма Bi-CGStab P) – как приближенное решение системы . Получается, что вместо того, чтобы точно решать вспомогательную систему с матрицей предобуславливателя K, можно просто найти приближенное решение той же вспомогательной системы, но уже с основной матрицей A, применяя метод LR1. В этом случае строки 9 и 13 алгоритма Bi-CGStab P должны быть заменены на следующие: … 9'. определение с помощью LR1 вектора приближенного решения y системы ; … 13'. определение с помощью LR1 вектора приближенного решения z системы ; … Иными словами идеологически (но не математически!) нет принципиальной разницы между тем, чтобы точно решить приближенную систему или тем, чтобы приближенно решить точную. Поскольку эффективность прямого сочетания двух методов зависит от точности приближения векторов y и z при минимальных затратах машинного времени для их нахождения, то отсюда следует, что лучше всего использовать такой метод, который за первую итерацию способен резко понизить начальную погрешность решения (соответственно и невязку). В этом смысле применение LR1 представляется наиболее оправданным. Дело в том, что данный метод характеризуется «фирменной» особенностью: он резко (на 2-4 порядка) понижает первоначальную невязку за первый итерационный проход [6-8]. В итоге на основании изложенных рассуждений можно определить новый метод, представляющий собой прямое сочетание алгоритмов LR1 и Bi-CGStab P, в дальнейшем для удобства называемый LR1sK. ^ Анализ эффективности. Для оценки эффективности метода LR1sK был проведен сравнительный анализ этого метода с методом Bi-CGStab P B путем решения СЛАУ, полученных разностной дискретизацией некоторых модельных задач. В качестве первой такой задачи была выбрана задача, сформулированная выше, но с сеточным разбиением области 1001×1001, что привело к увеличению числа обусловленности матрицы системы до МА  2,0×1010. Результаты расчета методами Bi-CGStab P B и прямого сочетания алгоритмов LR1 и Bi-CGStab P – LR1sK приведены на рис. 2. Алгоритму LR1sK потребовалось всего 23 итерации и времени 13,4 с для достижения точности e = 10-8, а алгоритму Bi-CGStab P B – соответственно 180 итераций и 51,2 с. Таким образом, предположение об ускорении алгоритма LR1 за счет его прямого сочетания с алгоритмом Bi-CGStab P оправдалось. В статье [4] вторая модельная задача формулируется следующим образом: «Система линейных уравнений следует из пятиточечной конечно-разностной дискретизации уравнения в частных производных на единичном квадрате с условием Дирихле на границе y = 0 и условиями Неймана на остальных границах. Функция D определяется как D = 1000 для 0,1  x, y  0,9 и D = 1 в оставшейся части области». Для проведения расчетов необходимо было доопределить постановку конкретными граничными условиями: . Оценка числа обусловленности матрицы системы в поставленной задаче при сеточном разбиении области 1001×1001 составила МА  6×108. На рис. 3 представлено поведение отношения норм невязок в зависимости от номера итерации. Заметны преимущества метода LR1sK (k = 25 итераций; t = 15,0 с) перед методом Bi-CGStab P B (k = 180; t = 50,8 с). Рис. 2 – Поведение отношения норм невязок Рис. 3 – Поведение отношения норм невязок 1 – LR1sK; 2 – Bi-CGStab P B 1 – LR1sK; 2 – Bi-CGStab P B В статье [4] приводится еще одна модельная задача, которая формулируется следующим образом: «Система линейных уравнений с несимметричной матрицей получена путем дискретизации краевой задачи на единичном квадрате с граничными условиями Дирихле на всех границах. При этом а = 20 exp[3,5 (x2+y2)]». Как и предыдущая, эта задача была доопределена конкретным условием u = 1 на всех границах области. Оценка числа обусловленности при сеточном разбиении 1001×1001 дает величину МА  6,8×105. Результаты расчетов данной задачи приведены на рис. 4. Хорошо видно, что алгоритму LR1sK потребовалось всего k = 9 итераций для достижения точности e = 10-8, при этом время, затраченное на эти итерации, составило t = 5,08 с против количества итераций, затраченных методом Bi-CGStab P B, k = 484 и времени t=132,5 с. Рис. 4 – Поведение отношения норм невязок Рис. 5 – Поведение отношения норм невязок 1 – LR1sK; 2 – Bi-CGStab P B 1 – LR1sK; 2 – Bi-CGStab P B В заключение следует отметить еще одну особенность метода LR1sK, выгодно отличающую его от исходного LR1. По определению алгоритм полинейного рекуррентного метода LR1 требует, чтобы коэффициенты исходной матрица системы уравнений A удовлетворяли условию строчного диагонального преобладания. Прямое же сочетание этого алгоритма LR1 с алгоритмом бисопряженных градиентов, представляющий метод LR1sK, позволяет ослабить это требование, так как градиентные методы способны находить решение систем, не удовлетворяющих данному условию. Иными словами LR1sK в состоянии устойчиво строить решения систем уравнений, коэффициенты которых не удовлетворяют условию строчного диагонального преобладания. Для иллюстрации этого утверждения решена следующая задача: в единичном квадрате поведение искомой функции определяется уравнением Гельмгольца ; (3) на границах области задается граничное условие первого рода . Здесь SC = (k2–2 2) u, = 1; u (x, y) = sin( x) sin( y) – точное решение задачи. Понятно, что поскольку k2 > 0, разностная дискретизация задачи приводит к СЛАУ с матрицей, у которой отсутствует диагональное преобладание коэффициентов. Результаты решения задачи (3) методами LR1sK и Bi-CGStab P B при сеточном разбиении области равномерной сеткой 1001×1001 (МА 107), заданной точности e = 10-8 и k2 =103 представлены на рис. 5. Необходимо отметить, что для достижения данной точности методу LR1sK потребовалось 193 итерации и времени 108,5 с, а методу Bi-CGStab P B – 2479 итераций и 170,7 с соответственно. Выводы. На основе представленных результатов расчетов можно сделать вывод, что новый метод LR1sK позволяет не только решать более жесткие СЛАУ (по отношению к методу LR1), но и тратить на их решение заметно меньше времени. Показано превосходство этого метода по быстродействию как над исходным полинейным рекуррентным методом, так и над одним из наиболее эффективных на данный момент методов бисопряженных градиентов со стабилизацией с предобуславливателем на базе явного метода Н.И. Булеева. Список литературы Yousef Saad. Iterative Methods for Sparse Linear Systems / Y. Saad. – N.Y.: PWS Publ. – 1996. – 460 p. Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений / В.П. Ильин. – Новосибирск: Изд-во ин-та математики. – 2000. – 345 c. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости / С. Патанкар. – М.: Энергоатомиздат. – 1984. – 152 c. Van der Vorst H.A. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems / H.A. van der Vorst // SIAM J. Sci. Stat. Comput. – 1992. – vol. 13. – № 2. – pp. 631–644. Старченко А.В. Сравнительный анализ некоторых итерационных методов для численного решения пространственной краевой задачи для уравнений эллиптического типа // Вестник ТГУ, Бюллетень оперативной научной информации.– Томск: ТГУ.–2003.– № 10.– C. 70-80. Фомин А.А. Сравнение эффективности высокоскоростных методов решения разностных эллиптических СЛАУ / А.А. Фомин, Л.Н. Фомина // Вестник Томского государственного университета. Математика и механика. – 2009. – № 2. – C. 71–77. Фомина Л.Н. Использование полинейного рекуррентного метода с переменным параметром компенсации для решения разностных эллиптических уравнений / Л.Н. Фомина // Вычислительные технологии. – ИВТ СО РАН. – 2009. – Т. 14. – № 4. – C. 108–120. Фомин А.А. Об одном варианте полинейного рекуррентного метода решения разностных эллиптических уравнений / А.А. Фомин, Л.Н. Фомина // Вестник Томского государственного университета. Математика и механика. – 2010. – № 2. – С. 20–27. Скачать 103,96 Kb. оставить комментарий Дата 02.10.2011 Размер 103,96 Kb. Тип Документы, Образовательные материалы