A NEW VERSION OF THE LINE-BY-LINE RECURSIVE METHOD FOR SOLVING DIFFERENCE ELLIPTICAL EQUATIONS.pdf Современные задачи вычислительной гидромеханики обычно сводятся к решению систем линейных алгебраических уравнений. Матрицы таких систем имеют высокий порядок и разреженную структуру специального вида [1]. Многие современные методы решения подобных СЛАУ обладают высокими скоростными характеристиками [2-6], которые, прежде всего, проявляются при решении систем с «хорошими» матрицами. Одним из таких эффективных методов решения систем уравнений с матрицами положительного типа [7] является полинейный рекуррентный метод, различные варианты которого позволяют понижать первоначальную норму невязки систем с количеством уравнений ~105-^106 на 1(Н14 порядков за 1(Н102 итераций [8, 9]. Однако на практике нередко встречаются задачи с различного рода осложнениями в своей исходной постановке: например, со скачкообразно меняющимися коэффициентами при вторых производных, - которые при своих дискретизациях сводятся к СЛАУ с плохо обусловленными матрицами. В подобных ситуациях даже самые лучшие методы начинают испытывать трудности при построении решений для таких систем и тогда вопрос ставится не о быстром (эффективном) решении СЛАУ, а о решении системы вообще [10]. Исходя из этого, в настоящей работе рассматривается новый вариант полинейного рекуррентного метода, который обладает улучшенными «разрешающими» способностями по сравнению со своими предыдущими вариантами [8 - 10].Постановка задачиПусть имеет место двумерная по пространству краевая задача на прямоугольной области Q = {(х,у): 0 < х < Lx, 0 < у < L } . Внутри области Q поведение искомой функции Ф(х,у) описывается дифференциальным уравнениемun+v™=±rx™)+±fya*}+Stдх ду дх\ дх) ду\ ду)Об одном варианте полинейного ренуррентного метода решения21где U, V- компоненты скорости, vx ,vy - коэффициенты переноса, S - источник. На границе области Г в общем случае имеют место условия третьего родадф и ^аг\-ог Ф = сг,дп где aY, bY, сг - известные величины, а и - нормаль к границе.Разностная аппроксимация подобной задачи приводит к системе линейных уравнений с пятидиагональной матрицей положительного типа. Основная суть полинейного рекуррентного метода состоит в последовательности преобразований исходной системы уравнений с целью приведения ее к виду, в котором матрица системы имеет четырехдиагональную структуру (рис. 1).XX XXX XXX XXX XXX X X X XглX X X X XXX XXX XXX XXX XXX X X X XX)X X X X XXX XXX XXX XXX XXX ^-X X X X()X X X X XXX ххх XXX ххх XXX X X X XL/X X X X XXX XXX XXX XXX XXОXX ххх ххх XXX XXX X X X XглXX ххх XXX XXX XXX X X X XX)XX XXX XXX XXX XXX ^-X X X X()XX XXX XXX XXX XXX X X X XL/X X XXX XXX XXX XXРис. 1Причем практически на всех этапах этой последовательности проводимые преобразования являются эквивалентными, кроме одного этапа, на котором значение поправки решения в так называемом «внешаблонном» узле приближенно выражается через поправки решения в узлах основного разностного шаблона. Легко показать, что данный прием по своей сути аналогичен применению принципа компенсации Н.И. Булеева [7], который, как известно, заключается в добавлении к уравнению таких итерируемых членов, которые приводят к взаимному исключению ошибок итерационных приближений, если они являются гладкими функциями координат.Модификация методаНетрудно видеть, что преобразованную таким образом систему можно разрешить последовательными трехточечными скалярными прогонками по клеткам матрицы [9]. В вариантах полинейного рекуррентного метода с линейной (LR1) и квадратичной (LR2) экстраполяцией поправки решения искомое значение Фу во «внешаблонном» узле выражается соответственно следующим образом:-Jc+\ф:к+\:Фу+(к+\Z/+1[2(Ф"ij+2Jф+1)-(ф;2-ф+2)];(1)ф:к+\Фк +(У[з (ф£1 ■-к+\,-y+3J>[3(Ф*.-ф£+2)+ф£+з].(2)Здесь 8 - итерационный параметр компенсации (0 < 8 < 1), к - номер итерации, а экстраполяция проводится вдоль линии /' = const сеточного разбиения области ре-22А.А.Фомин,Л.Н. Фоминашения. Формулы (1), (2), являясь, по сути, разложением искомой функции в ряд Тейлора по рядом лежащим узлам, позволяют замкнуть итерационный алгоритм построения решения, не добавляя никакой новой информации и, тем самым, не переопределяя исходную систему уравнений. Легко видеть, что в случае сходимости решения lim Фг = Фг они вырождаются в тождества типа Фг = Фг . Основ-ным недостатком этих формул является их приближенный характер на протяжении всего итерационного процесса. Следовательно, при неудачном стечении обстоятельств, например при решении СЛАУ с плохо обусловленной матрицей, ошибка, порождаемая приближенным соотношением (1) или (2) может на некотором этапе сходимости решения привести к неулучшаемой ситуации. Например, в [10] показано как при решении тестовой задачи с плохо обусловленной матрицей отношение норм текущей невязки к первоначальной быстро достигало величины е0~2,0х 10~5 и при дальнейших итерациях уже никак не менялось.Как уже отмечалось, в цепочке последовательных преобразований системы полинейного рекуррентного метода только один этап является приближенным. В данном случае это использование формулы (1) или формулы (2). Если бы вместо них было бы использовано некоторое точное соотношение, связывающее компоненты искомого вектора решения, то метод утратил бы свой итерационный характер и стал бы прямым. В силу того, что построить такое соотношение на данный момент не удается, можно принять компромиссный вариант, а именно - использовать замыкающее соотношение с переменными коэффициентами, которые по мере исполнения итераций стремились бы к некоторым предельным величинам, точно связывающим точные значения компонентов вектора решения. В качестве таких замыкающих соотношений удобно использовать двухточечные связи из метода В.Г. Зверева [3], имеющие следующий вид:Ф^^Ф^+Л*-,(3)здесь |г , г|г - коэффициенты, принимающие асимптотически точные значения в случае сходимости решения. То есть, если имеют место пределы lim Й =£,*,-, lim г|* = г|*, lim Ф* = Ф* , то точно выполняются соотношения%=^%+1+Щ-(4)Следуя методике [3], несложно получить расчетные формулы вычисления |г ,г|г . Пусть разностная схема исходной СЛАУ, построенная методом контрольногообъема [11] для узла (/',/) имеет вид%%' =%ф.--1у+%Ф,-+1у+%Ф(,--1+%.Ф«+1+^-(5)Используя механизм принципа компенсации, можно добавить в обе части (5) выражение -&(aw.. +aE.. )Фг . Тогда с учетом (3) уравнение (5) приобретет вид-6(Ч +%)Ф*.+%(^._1Ф,+л|-1) + 6,,(6)или окончательноОб одном варианте полинейного ренуррентного метода решения23[an. -6{ащ. + аЩ]) -^._1%у. ]Фу = аЩ] Ф1]+1 ++ [bv +4-1% +%Ф*-1У +%■■ уточняется экстра-поляционные соотношения, уточняется экстраполяционные соотношения ->■ уточняется решение и т.д. В то время как в вариантах LR\ и LR2 подобный механизм «одноходовой»: уточняется только решение, а коэффициенты экстраполяци-24А.А.Фомин,Л.Н. Фоминаонных соотношений остаются неизменно приближенными (см. (1) - (2)). Отсюда следует, что, в принципе, нельзя исключить прекращение процесса сходимости решения в районе некоторого критического е0, что и произошло при решении тестовой задачи в работе [10].Вычислительный эксперимент. Для сравнительного исследования свойств варианта LRz была решена следующая краевая задача Дирихле в единичной квадратной области:(yxOx)x+(yyOy)y=S(x,y),(9)Ф|г = 0. Коэффициенты при старших производных и точное решение задачи известны:vx=l + 2[(x-0.5)2+(>>-0.5)2], vy =l + 2[0.5-(x-0.5)2-(>>-0.5)2],Ф(х,у) = 256[ху(1-х)(1-у)]2.Правая часть S(x,y) вычислялась путем подстановки vx, vy и Ф(х,у) в уравнение (9). Сеточное разбиение области составляло 401x401=160 801 узлов. Сетка равномерная. Решение считалось найденным при выполнении условия сходимости Ш h /\\R II2 < е , где R ,R° - соответственно текущая и начальная невязки, а е - заданная точность решения, которая в данном и последующих расчетах принималась равной 10~10. На рис. 2 представлены зависимости отношений норм невязок от номера итерации при решении данной задачи вариантами полинейного рекуррентного метода: LR1 (6=0,99960), LR2 (6 = 0,9999970) и LRz (6 = 0,9999250). В скобках приведены значения итерационного параметра компенсации для каждого из вариантов, позволивших достичь максимальной скорости сходимости. Тем самым была проведена своеобразная «оценка сверху» возможностей каждого из вариантов полинейного рекуррентного метода. Видно, что кривые вариантов LR\ и LR2 ведут себя почти одинаково и для достижения необходимой точности потребовалось по 25 итераций, в то время как вариант LRz построил решение за 15 итераций.20 30 кРис.2Рис.3Об одном варианте полинейного ренуррентного метода решения25Как известно, важной характеристикой процесса сходимости является величина средней скорости сходимости, вычисляемая по формуле [3]0*=-!in(|z*| /||z°|| ),(10)^£ VII \\2 II 112/V 'где Z , Z0 - соответственно текущая и начальная погрешности решения. На рис. 3 сравниваются кривые средних скоростей сходимости вариантов LR1, LR2 и LRz. Видно, что графики вариантов LR1 и LR2, так же как и на рис. 2, ведут себя почти одинаково, что вполне естественно. А график LRz расположен заметно выше, что тоже понятно, поскольку этот расчет по количеству итераций сошелся почти в два раза быстрее. Обращает на себя внимание тот факт, что по порядку величины в рассматриваемых расчетах средняя скорость сходимости Q может быть оцененакак 0(1). Более подробное исследование вопроса зависимости средней скорости сходимости от количества решаемых уравнений показало, что при увеличении числа уравнений в 100 раз: сетка 41x41 ->■ сетка 401x401, - средняя скорость сходимости уменьшается в 3 раза: Qk « 4,0 ->■ Qk « 1,4.Следует заметить, что минимальное число итераций еще не означает минимальные затраты машинного времени. В проведенных расчетах на варианты с применением LR\ и LR2 потребовалось 6,0 с, а на вариант LRz - 8,5 с. Получается, что на проведение одной итерации варианту LRz требуется приблизительно в два раза больше времени, чем вариантам LR\ и LR2.Для демонстрации улучшенных «разрешающих» свойств варианта LRz были решены две задачи повышенной сложности из [4]. Первая из этих задач формулируется в единичном квадрате следующим образом:-(DOx)x-(DOy)y=l;(11)дФдхдФФ\у=0= -1, ™ =1.i,-1. j,-,х=09х х=1Причем функция D определяется как D = 1000 для 0,1 < х, у < 0,9 и D = 1 в оставшейся части области. Шаг равномерной сетки равен 1/150. Графики отношений норм невязок в зависимости от номера итерации представлены на рис. 4. Видно, что вариант LR\ (6 = 0,999988) сошелся приблизительно за 160 итераций, вариант LRz (6 = 0,999930) -за» 1900 итераций, а вариант LR2 (6 = 0,999990) за 2000 итераций достиг отношения \\R h/llR lb ~ 10~4. Продолжение счета вариантом LR2 позволили достичь значений по точности ~ 10~8 за 5000 итераций, после чего вычислительный процесс был прекращен. Следует заметить, что оригинальным методом бисопряженных градиентов [4] решить данную задачу не удалось по причине расходимости решения [10].Коэффициенты уравнения определяются следующим образом: В(х, у) = = 2 ехр[2 (х2 + у2)]; F = 100 в небольшой квадратной области в центре единичного квадрата, в основной его части F = 0; коэффициенте по области меняется кусочно-постоянным образом и его значение скачкообразно варьируется от 10~5 до 104 согласно схеме, приведенной в [4] или [10]. Шаг равномерной сетки составляет 1/128. Как и в предыдущей задаче, расчеты проводились вариантами полинейного рекуррентного метода LR1 (6 = 0,950), LR2 (6 = 0,550), LRz (6 = 0,999927). На26А.А.Фомин,Л.Н. Фоминарис. 5 изображены соответствующие кривые сходимости решения. Видно, что отношения 11л*112/11л°11 2 для вариантов LRI и LR2 достаточно быстро, через 30 -ь 100 итераций, достигают значения ~ 2,0х10~5, после чего перестают уменьшаться. В то время как вариант LRz свел решение к требуемой точности за 1680 итераций. Обращают на себя внимание относительно невысокие значения итерационного параметра компенсации для вариантов LR\ и LR2. Данные значения 8 являются наибольшими, при которых расчеты выполнялись устойчиво с максимальной скоростью сходимости. При дальнейших малых увеличениях параметра компенсации происходил аварийный останов вычислительного процесса. Также следует отметить, что в отличии от предыдущей, данная задача была успешно решена различными вариантами метода бисопряженных градиентов за ~ 70 -ь 120 итераций.3 \gk3 1gЈРис.4'■ ■ -хLR2N LRIГЛч LRz\s 10"10\ \12Рис.5Вторая краевая задача повышенной сложности формулируется в единичном квадрате следующим образом:-(АФх)х-(АФу)у+В(х,у)Фх=Е;(12)У'У0, Ф|1.Ф|1, ф|I ф|\у=1\у=0ВыводыНа основании изложенного материала можно сделать следующие выводы.Разработан вариант LRz полинейного рекуррентного метода с асимптотически точными замыкающими соотношениями.Для достижения заданной точности решения для варианта LRz требуется приблизительно столько же машинного времени, сколько и для вариантов LR1 и LR2 полинейного рекуррентного метода.Время, необходимое для реализации одной итерации вариантом LRz при прочих равных условиях, приблизительно в два раза превосходит аналогичное время вариантов LR1 и LR2.Вариант LRz обладает более высокими «разрешающими» способностями по отношению к вариантам LR1 и LR2 безотносительно к их скорости сходимости.

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