CC BY
20
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2010
Математика и механика
№ 2(10)
УДК 519.632.4
А.А. Фомин, Л.Н. Фомина
ОБ ОДНОМ ВАРИАНТЕ ПОЛИНЕЙНОГО РЕКУРРЕНТНОГО МЕТОДА РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ
Рассматривается новый вариант замыкающего соотношения полинейного рекуррентного метода, представляющий собой двухточечную линейную связь значений искомого решения в соседних узлах сеточного разбиения области. Особенностью рассматриваемого замыкания являются асимптотически точные значения коэффициентов соотношения в случае сходимости решения. На примерах решения модельных задач показана более высокая эффективность данного варианта полинейного рекуррентного метода решения систем линейных алгебраических уравнений (СЛАУ) с плохо обусловленными матрицами по сравнению с его предыдущими вариантами.
Ключевые слова: система линейных алгебраических уравнений, итерационный метод решения, замыкающие соотношения.
Современные задачи вычислительной гидромеханики обычно сводятся к решению систем линейных алгебраических уравнений. Матрицы таких систем имеют высокий порядок и разреженную структуру специального вида [1]. Многие современные методы решения подобных СЛАУ обладают высокими скоростными характеристиками [2 - 6], которые, прежде всего, проявляются при решении систем с «хорошими» матрицами. Одним из таких эффективных методов решения систем уравнений с матрицами положительного типа [7] является полинейный рекуррентный метод, различные варианты которого позволяют понижать первоначальную норму невязки систем с количеством уравнений ~105^106 на 10^14 порядков за 10^102 итераций [8, 9]. Однако на практике нередко встречаются задачи с различного рода осложнениями в своей исходной постановке: например, со скачкообразно меняющимися коэффициентами при вторых производных, - которые при своих дискретизациях сводятся к СЛАУ с плохо обусловленными матрицами. В подобных ситуациях даже самые лучшие методы начинают испытывать трудности при построении решений для таких систем и тогда вопрос ставится не о быстром (эффективном) решении СЛАУ, а о решении системы вообще [10]. Исходя из этого, в настоящей работе рассматривается новый вариант полинейного рекуррентного метода, который обладает улучшенными «разрешающими» способностями по сравнению со своими предыдущими вариантами [8 - 10].
Пусть имеет место двумерная по пространству краевая задача на прямоугольной области О = {(х,у): 0 < х < Ьх, 0 < у < Ьу }. Внутри области О поведение искомой функции Ф(х,у) описывается дифференциальным уравнением
Постановка задачи
где и, V - компоненты скорости, Vх, Vу - коэффициенты переноса, £ - источник. На границе области Г в общем случае имеют место условия третьего рода
дФ . ^
аг-----+ Ъг Ф = сг,
дп
где аГ, ЪГ, сГ - известные величины, а п - нормаль к границе.
Разностная аппроксимация подобной задачи приводит к системе линейных уравнений с пятидиагональной матрицей положительного типа. Основная суть полинейного рекуррентного метода состоит в последовательности преобразований исходной системы уравнений с целью приведения ее к виду, в котором матрица системы имеет четырехдиагональную структуру (рис. 1).
Рис. 1
Причем практически на всех этапах этой последовательности проводимые преобразования являются эквивалентными, кроме одного этапа, на котором значение поправки решения в так называемом «внешаблонном» узле приближенно выражается через поправки решения в узлах основного разностного шаблона. Легко показать, что данный прием по своей сути аналогичен применению принципа компенсации Н.И. Булеева [7], который, как известно, заключается в добавлении к уравнению таких итерируемых членов, которые приводят к взаимному исключению ошибок итерационных приближений, если они являются гладкими функциями координат .
Модификация метода
Нетрудно видеть, что преобразованную таким образом систему можно разрешить последовательными трехточечными скалярными прогонками по клеткам матрицы [9]. В вариантах полинейного рекуррентного метода с линейной (LR1) и квадратичной (LR2) экстраполяцией поправки решения искомое значение Ф. во «внешаблонном» узле выражается соответственно следующим образом:
ф|+1 =ф|-+0[2(Фк;:1 -ф*-+1) -(Фк+2-ф.+2)]; (1)
ф^+1 = ф. +0 [3 (фк+1 - фк+2 ) + фк+3 ] - 0 [3(ф.+1 -ф|.+ 2) + фк+3 ] . (2)
Здесь 0 - итерационный параметр компенсации (0 < 0 < 1), k - номер итерации, а экстраполяция проводится вдоль линии i = const сеточного разбиения области ре-
шения. Формулы (1), (2), являясь, по сути, разложением искомой функции в ряд Тейлора по рядом лежащим узлам, позволяют замкнуть итерационный алгоритм построения решения, не добавляя никакой новой информации и, тем самым, не переопределяя исходную систему уравнений. Легко видеть, что в случае сходимо-
к * * * сти решения пт Ф,7 = Ф,7 они вырождаются в тождества типа Ф.7 = Ф„-. Основах 1 1 1 1
ным недостатком этих формул является их приближенный характер на протяжении всего итерационного процесса. Следовательно, при неудачном стечении обстоятельств, например при решении СЛАУ с плохо обусловленной матрицей, ошибка, порождаемая приближенным соотношением (1) или (2) может на некотором этапе сходимости решения привести к неулучшаемой ситуации. Например, в [10] показано как при решении тестовой задачи с плохо обусловленной матрицей отношение норм текущей невязки к первоначальной быстро достигало величины е0~2,0х 10-5 и при дальнейших итерациях уже никак не менялось.
Как уже отмечалось, в цепочке последовательных преобразований системы полинейного рекуррентного метода только один этап является приближенным. В данном случае это использование формулы (1) или формулы (2). Если бы вместо них было бы использовано некоторое точное соотношение, связывающее компоненты искомого вектора решения, то метод утратил бы свой итерационный характер и стал бы прямым. В силу того, что построить такое соотношение на данный момент не удается, можно принять компромиссный вариант, а именно - использовать замыкающее соотношение с переменными коэффициентами, которые по мере исполнения итераций стремились бы к некоторым предельным величинам, точно связывающим точные значения компонентов вектора решения. В качестве таких замыкающих соотношений удобно использовать двухточечные связи из метода В.Г. Зверева [3], имеющие следующий вид:
Ф1+1 = ^ф1+1 +пк., (3)
здесь , п. - коэффициенты, принимающие асимптотически точные значения в
случае сходимости решения. То есть, если имеют место пределы
Пт Й , Пт п. = п*, 1™ Ф. = Ф* , то точно выполняются соотношения
к^х 1 к^<» к^х
А А А А
Ф1 =^ Ф1+1 +п . (4)
Следуя методике [3], несложно получить расчетные формулы вычисления ^к , пк . Пусть разностная схема исходной СЛАУ, построенная методом контрольного объема [11] для узла (.,.) имеет вид
% Ф1 = %■ Ф.-11 + %. Ф.+11 + % Ф-1 + %■ Ф1+1 + Ъ.1 . (5)
Используя механизм принципа компенсации, можно добавить в обе части (5)
выражение -9 (ат, + аЕ,, )Фг1. Тогда с учетом (3) уравнение (5) приобретет вид
1 1
% 9 (аЕ у + аШу )] Ф1 = аЫу Ф1+1 + аШу Фк-1 ] + а.
Е. 1 Ф к+11
9 (%■ + %■ )Фк + % ^-1 Фу + п1 -1 ) + Ъу , (6)
или окончательно
[ ар.1 - 9 (аЕ.1 + аШц ) - ^-1%- ] Фу = %■ Ф1+1 +
+ [ЪУ + п1-1% + % Фк-11 + %■ Фк+11 -9 (аЩ + %■ )Ф(/ ]. (7)
Сравнение (7) с (3) позволяет записать рекуррентные двухточечные связи:
= %■ 7 [%■ -9 (% + аЕу )-4(1-1а5у ] ;
к = Ъ1 + п1 -1% + Фк-11 + аЩ] Фк+11 - 9 (%- + аЩ] )Ф(8)
У % -9 (% +% ) Ч*-1 % '
Если итерационный параметр компенсации 9 зависит от номера итерации, например, при использовании технологии его адаптации под очередное приближение вектора решения {ф* } [9,12], то итерируются оба коэффициента в соотношении (3). В противном случае вектор {^} вычисляется один раз, а итерируется только вектор свободного члена Щ} . Нетрудно видеть, что при сходимости решения существуют и пределы коэффициентов двухточечных связей. Иными словами, если бы заранее были известны векторы {4*} и {п*}, то подстановка их в
общий алгоритм полинейного рекуррентного метода позволила бы сразу без итераций построить вектор решения. Отсюда следует, что совершенно не обязательно, чтобы 4*1 и п* рассчитывались по формулам (8). Подойдет любой алгоритм, в
том числе и не зависящий явно от очередного приближения фк, лишь бы он позволял получать предельные значения коэффициентов, корректные с точки зрения выполнения соотношения (4). В этом смысле рассматриваемый метод (в дальнейшем обозначаемый как ЬЯг) перекликается с основной идеей метода Б.Н. Четве-рушкина [13], в котором, как известно, итерируется не сам вектор решения, а коэффициенты некоторого вычислительного представления его компонент в виде их линейной взаимозависимости. И только после сходимости этих коэффициентов с заданной точностью производится окончательное вычисление компонент искомого вектора решения за один проход по расчетной области.
В методе ЬЯг итерируются как сам вектор решения, так и коэффициенты двухточечной связи замыкающего соотношения (3). Именно эта особенность обеспечивает данному варианту полинейного рекуррентного метода повышенные «разрешающие» свойства. Действительно, в варианте ЬЯг замыкающие соотношения по мере сходимости решения становятся практически точными, в то время как в рассмотренных ранее вариантах (ЬК\, ЬК2) они продолжают оставаться приближенными, по сути - разложениями в ряд Тейлора с первым (ЬК\) или вторым (ЬК2) порядком точности [8, 9]. Отсюда следует, что по мере сходимости решения экстраполяционная формула (3) в ЬЯг все более и более уточняется, что и приводит, в конечном счете, к достижению заданной точности решения безотносительно к количеству затраченных на это итераций. Иными словами, реализуется «двухходовый» механизм сходимости: уточняется решение ^ уточняется экстраполяционные соотношения, уточняется экстраполяционные соотношения ^ уточняется решение и т.д. В то время как в вариантах ЬК\ и ЬК2 подобный механизм «одноходовой»: уточняется только решение, а коэффициенты экстраполяци-
онных соотношений остаются неизменно приближенными (см. (1) - (2)). Отсюда следует, что, в принципе, нельзя исключить прекращение процесса сходимости решения в районе некоторого критического е0, что и произошло при решении тестовой задачи в работе [10].
Вычислительный эксперимент. Для сравнительного исследования свойств варианта ЬИх была решена следующая краевая задача Дирихле в единичной квадратной области:
(V * Ф х) х + (V у Фу)у = 5 (х, у), (9)
Ф|г= 0.
Коэффициенты при старших производных и точное решение задачи известны: Vх = 1 + 2[(х-0.5)2 + (у-0.5)2] , Vу = 1 + 2[0.5-(х-0.5)2-(у-0.5)2] ,
Ф(х, У ) = 256 [ху (1 - х )(1 - у )]2.
Правая часть 5(х,у) вычислялась путем подстановки Vх, Vу и Ф(х,у) в уравнение (9). Сеточное разбиение области составляло 401x401=160 801 узлов. Сетка равномерная. Решение считалось найденным при выполнении условия сходимости 1И12/11 И0 ||2 <6 , где Ик, И0 - соответственно текущая и начальная невязки, а 6 - заданная точность решения, которая в данном и последующих расчетах принималась равной 10-10. На рис. 2 представлены зависимости отношений норм невязок от номера итерации при решении данной задачи вариантами полинейного рекуррентного метода: ЬИ1 (0=0,99960), ЬИ2 (0 = 0,9999970) и ЬИх (0 = 0,9999250). В скобках приведены значения итерационного параметра компенсации для каждого из вариантов, позволивших достичь максимальной скорости сходимости. Тем самым была проведена своеобразная «оценка сверху» возможностей каждого из вариантов полинейного рекуррентного метода. Видно, что кривые вариантов ЬИ1 и ЬЯ2 ведут себя почти одинаково и для достижения необходимой точности потребовалось по 25 итераций, в то время как вариант ЬЯг построил решение за 15 итераций.
Рис. 2
Рис. 3
Как известно, важной характеристикой процесса сходимости является величина средней скорости сходимости, вычисляемая по формуле [3]
сравниваются кривые средних скоростей сходимости вариантов ЬR1, ЬR2 и ЬRz. Видно, что графики вариантов LR1 и LR2, так же как и на рис. 2, ведут себя почти одинаково, что вполне естественно. А график ЬRz расположен заметно выше, что тоже понятно, поскольку этот расчет по количеству итераций сошелся почти в два раза быстрее. Обращает на себя внимание тот факт, что по порядку величины в рассматриваемых расчетах средняя скорость сходимости Qk может быть оценена
как 0(1). Более подробное исследование вопроса зависимости средней скорости сходимости от количества решаемых уравнений показало, что при увеличении числа уравнений в 100 раз: сетка 41*41 ^ сетка 401*401, - средняя скорость сходимости уменьшается в 3 раза: Qk и 4,0 ^ Qk и 1,4.
Следует заметить, что минимальное число итераций еще не означает минимальные затраты машинного времени. В проведенных расчетах на варианты с применением LR1 и ЬR2 потребовалось 6,0 с, а на вариант ЬRz - 8,5 с. Получается, что на проведение одной итерации варианту ЬRz требуется приблизительно в два раза больше времени, чем вариантам LR1 и LR2.
Для демонстрации улучшенных «разрешающих» свойств варианта ЬRz были решены две задачи повышенной сложности из [4]. Первая из этих задач формулируется в единичном квадрате следующим образом:
Причем функция D определяется как D = 1000 для 0,1 < х, у < 0,9 и D = 1 в оставшейся части области. Шаг равномерной сетки равен 1/150. Графики отношений норм невязок в зависимости от номера итерации представлены на рис. 4. Видно, что вариант LR1 (0 = 0,999988) сошелся приблизительно за 160 итераций, вариант ЬRz (0 = 0,999930) - за и 1900 итераций, а вариант ЬR2 (0 = 0,999990) за 2000 итераций достиг отношения 1^112 /IIR0IІ2 ~ 10-4. Продолжение счета вариантом ЬR2 позволили достичь значений по точности ~ 10-8 за 5000 итераций, после чего вычислительный процесс был прекращен. Следует заметить, что оригинальным методом бисопряженных градиентов [4] решить данную задачу не удалось по причине расходимости решения [10].
Коэффициенты уравнения определяются следующим образом: В(х, у) = = 2 ехр[2 (х2 + у2)]; F = 100 в небольшой квадратной области в центре единичного квадрата, в основной его части F = 0; коэффициент А по области меняется кусочно-постоянным образом и его значение скачкообразно варьируется от 10-5 до 104 согласно схеме, приведенной в [4] или [10]. Шаг равномерной сетки составляет 1/128. Как и в предыдущей задаче, расчеты проводились вариантами полинейного рекуррентного метода LR1 (0 = 0,950), ЬR2 (0 = 0,550), ЬRz (0 = 0,999927). На
(10)
где Ък, Z - соответственно текущая и начальная погрешности решения. На рис. 3
(DФх)х - (DФу)у = 1;
(П)
, дФ , дФ , дФ
Ф _п = 1, ---------- = 1, ----- =-1, --------- = 1.
у ду у=1 дх х=0 дх х=1
рис. 5 изображены соответствующие кривые сходимости решения. Видно, что отношения 1И |2/| И|2 для вариантов LR1 и LR2 достаточно быстро, через 30 -г- 100 итераций, достигают значения ~ 2,0х10-5, после чего перестают уменьшаться. В то время как вариант LRz свел решение к требуемой точности за 1680 итераций. Обращают на себя внимание относительно невысокие значения итерационного параметра компенсации для вариантов LR1 и LR2. Данные значения 9 являются наибольшими, при которых расчеты выполнялись устойчиво с максимальной скоростью сходимости. При дальнейших малых увеличениях параметра компенсации происходил аварийный останов вычислительного процесса. Также следует отметить, что в отличии от предыдущей, данная задача была успешно решена различными вариантами метода бисопряженных градиентов за ~ 70 -г- 120 итераций.
Рис. 4 Рис. 5
Вторая краевая задача повышенной сложности формулируется в единичном квадрате следующим образом:
-(АФ*)x - (АФу )у + B(x, у) Фx = F ; (12)
Ф у=1 = 0 Ф у=0 = 1 Ф *=1 = 1 Ф x=o =1.
Выводы
На основании изложенного материала можно сделать следующие выводы.
Разработан вариант LRz полинейного рекуррентного метода с асимптотически точными замыкающими соотношениями.
Для достижения заданной точности решения для варианта LRz требуется приблизительно столько же машинного времени, сколько и для вариантов LR1 и LR2 полинейного рекуррентного метода.
Время, необходимое для реализации одной итерации вариантом LRz при прочих равных условиях, приблизительно в два раза превосходит аналогичное время вариантов LR1 и LR2.
Вариант LRz обладает более высокими «разрешающими» способностями по отношению к вариантам LR1 и LR2 безотносительно к их скорости сходимости.
ЛИТЕРАТУРА
1. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. М.: Эдиториал УРСС, 1999. 248 с.
2. Yousef Saad. Iterative Methods for Sparse Linear Systems. N.Y.: PWS Publ., 1996. 460 p.
3. Зверев В.Г. Модифицированный полинейный метод решения разностных эллиптических уравнений // ЖВМ и МФ. 1998. Т. 38. № 9. С. 1553 - 1562.
4. 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. № 2. P. 631 - 644.
5. Старченко А.В. Сравнительный анализ некоторых итерационных методов для численного решения пространственной краевой задачи для уравнений эллиптического типа // Вестник ТГУ. Бюллетень оперативной научной информации. Томск: ТГУ, 2003. № 10. C. 70 - 80.
6. Ильин В.П. Методы бисопряженных направлений в подпространствах Крылова // Сибирский журнал индустриальной математики. 2008. Т. 11. № 4. С. 47 - 60.
7. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995. 288 c.
8. Фомин А.А., Фомина Л.Н. Полинейный рекуррентный метод решения СЛАУ с пятидиагональной матрицей // Четвертая Сибирская школа-семинар по параллельным и высокопроизводительным вычислениям (Томск, 9 - 11 октября 2007 г.). Томск: Дельтаплан, 2008. C. 192 - 201.
9. Фомина Л.Н. Использование полинейного рекуррентного метода с переменным параметром компенсации для решения разностных эллиптических уравнений // Вычислительные технологии. ИВТ Со РАН. 2009. Т. 14. № 4. C. 108 - 120.
10. Фомин А.А., Фомина Л.Н. Сравнение эффективности высокоскоростных методов решения разностных эллиптических СЛАУ // Вестник ТГУ. Математика и механика. 2009. № 2. C. 71 - 77.
11. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.
12. Zverev V.G. About the iteration method for solving difference equations // Lecture notes in computer science. Berlin; Heidelberg: Springer-Verlag, 2005. V. 3401. P. 621 - 628.
13. Четверушкин Б.Н. Об одном итерационном алгоритме решения разностных уравнений // ЖВМ и МФ. 1976. Т. 16. № 2. C. 519 - 524.
СВЕДЕНИЯ ОБ АВТОРАХ:
ФОМИН Александр Аркадьевич - кандидат физико-математических наук, руководитель
отдела информационных технологий, ОАО «Издательско-полиграфическое предприятие
«Кузбасс». E-mail: [email protected]
ФОМИНА Любовь Николаевна - старший преподаватель кафедры вычислительной математики Кемеровского государственного университета. E-mail: [email protected]
Статья принята в печать 10.05.2010 г.
