Использование полинейного рекуррентного метода с переменным параметром компенсации для решения разностных эллиптических уравнений Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 14, № 4, 2009

Использование полинейного рекуррентного метода

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

Л. Н. ФОМИНА Кемеровский государственный университет, Россия e-mail:lubafomina@mail .ru

Приводится алгоритм итерационного полинейного рекуррентного метода решения пятидиагональной системы линейных алгебраических уравнений (СЛАУ), возникающей при разностной аппроксимации двумерных уравнений динамики вязкой жидкости и тепломассопереноса. Тестовая задача решается модифицированным полинейным и полинейным рекуррентным методами как с постоянным, так и c переменным параметрами компенсации. Анализируется эффективность использования переменного параметра компенсации в сравнении с постоянным.

Ключевые слова: система линейных алгебраических уравнений, итерационный метод решения, переменный параметр компенсации.

Введение

Современные методы решения краевых задач, связанных с исследованием явлений гидродинамики и тепломассопереноса, сводятся, как правило, к разностной аппроксимации многомерных дифференциальных уравнений, что в свою очередь приводит к построению систем линейных алгебраических уравнений, матрицы которых имеют большую размерность и разреженно-упорядоченную структуру [1, 2]. Для одномерных по пространству задач матрицы СЛАУ имеют трехдиагональную структуру, решение таких задач легко находится прямым экономичным методом скалярной прогонки, в то время как для двумерных задач число диагоналей возрастает, как правило, до пяти, что существенно осложняет процедуру решения системы. В первую очередь эти осложнения связаны с необходимостью использования итерационных подходов при построении методов решения подобных СЛАУ [1-4].

Повышение эффективности итерационных методов решения СЛАУ, согласно идее, высказанной еще в монографии [3], определяется скоростью передачи информации от границ расчетной области ко всем ее внутренним узлам. В свою очередь скорость передачи информации и, как следствие, скорость сходимости напрямую зависят от степени неявности алгоритма метода [1-5]. По этой причине поточечные методы, такие как, например, метод последовательной верхней релаксации, хотя и позволяют достаточно надежно строить приближенное решение СЛАУ, но делают это медленно, поскольку за одну итерацию информация передается всего на один узел сеточного пространства.

В полинейных методах скорость сходимости удается повысить путем комбинации неявных прогонок вдоль сеточных линий выбранного направления и поточечного — по

© ИВТ СО РАН, 2009.

остальным направлениям [3]. Однако из-за того, что лишь по одной пространственной координате информация передается за одну итерацию, а по остальным поперечным направлениям действует все тот же механизм поточечной передачи, скорость сходимости этих методов продолжает оставаться относительно невысокой.

Другой подход, также позволяющий заметно повысить эффективность итерационных методов, представляет собой так называемый принцип компенсации Булеева [2]. Согласно этому принципу компоненты вектора решения, входящие в правую часть итерационной записи уравнений и, следовательно, берущиеся с предыдущей итерации, приближенно выражаются через компоненты, определяемые на текущей итерации. Тем самым увеличивается степень неявности разрешающего оператора и соответственно повышается скорость сходимости итераций.

В работе [5] предложена модификация полинейного метода, суть которой заключается в том, что при построении решения вдоль выбранной сеточной линии с поперечного направления используются не сами значения искомой функции с предыдущей итерации, а коэффициенты двухточечной связи, характеризующие взаимоотношения значений функции. Поэтому, хотя и приближенно, информация о решении со всего поперечного направления используется при построении решения вдоль всего продольного направления в пределах одной итерации. Такой подход в сочетании с использованием принципа компенсации позволяет на порядок сократить необходимое число итераций по сравнению с классическим вариантом полинейного метода.

В работе [6] также решена проблема одномоментной передачи информации от всех границ ко всем внутренним узлам области за одну итерацию путем построения практически полностью неявного алгоритма решения СЛАУ. Несмотря на свою относительную сложность, при решении тестовой задачи предлагаемый метод позволяет на различных сетках понижать начальную норму невязки на три-четыре порядка за одну-две итерации.

В работе [7] в целях повышения эффективности модифицированного полинейного метода [5] используется технология подбора компенсационного значения итерационного параметра индивидуально для каждого разностного уравнения на каждом итерационном шаге. Суть идеи заключается в минимизации влияния очередного приближения решения с предыдущего итерационного шага на построение следующего приближения решения. Для этого в правой части каждого разностного уравнения выделяется сумма слагаемых, которая возникает вследствие итерационного характера записи уравнения и соответственно связана с приближением решения на предыдущей итерации. Эта сумма приравнивается к нулю, а из полученного уравнения вычисляется значение параметра компенсации. Результаты расчетов тестовой задачи показывают, что в зависимости от сеточного разбиения число итераций, необходимых для уменьшения начальной нормы невязки на четыре порядка, снижается в 3-10 раз по отношению к соответствующим расчетам с постоянным параметром.

В настоящей работе по аналогии с [7] рассматривается технология использования переменного параметра компенсации применительно к методу [6] и анализируются результаты ее применения.

1. Постановка задачи

Рассматривается единичная область П = {(ж, у): 0 < х < 1, 0 < у < 1}, внутри которой стационарное поведение искомой функции Ф описывается обобщенным диффе-

у+11

и —ср— 5и79-Г<^Г9зё

—<

Л NwX_L±iN±pE

на- •••

2

01"

7

'[I т

И /'

/+7

Рис. 1. Три типа шаблонов, возникающих при разностной аппроксимации уравнения и граничных условий

ренциальным уравнением:

5

6ФУ

д

д Фч

дх у х дх ) + ду у у ду ) а на границе области Г соответствует граничным условиям третьего рода:

д Ф , л

аг -—+ Ьг Ф = сг. 1 дп г г

(1)

(2)

Здесь аг' Ьг, сг — известные величины, а п — нормаль к границе.

Разностная аппроксимация (1) производится, в общем случае, по пятиточечному шаблону (рис. 1) методом контрольного объема [3] и имеет следующий вид:

аР„ Ф] = аЕц Фг+1] + Фг-1] + аМ„ Ф]+1 + Ф]-1 + Ь

г]-

(3)

Здесь 1 < г < п, 1 < ] < т, где п,т — количество узлов сеточного разбиения П по координатам х' у соответственно. Разностное уравнение на границе области получается путем аппроксимации уравнения (1) с учетом того, что поток Ф по нормали к границе выражается через (2).

В общем случае это будет уравнение, расписанное по четырехточечному шаблону (рис. 1). Например, на граничной линии х = 0 или (г =1) оно имеет вид

ари Ф] = ави Ф2] + ам^ Ф]+1 + а^. Ф1]-1 + Ь

>ц.

(4)

Аналогичным образом разностное уравнение в угловых узлах области П расписывается, в общем случае, по трехточечному шаблону (рис. 1). Например, в углу (0,0) или (г = 1'] = 1) разностное уравнение имеет вид

арп Ф11 = аЁ11 Ф21 + амп Ф12 + Ьп.

В итоге получается система линейных алгебраических уравнений с разреженной пятидиагональной матрицей.

2. Общая идея метода решения

Рассматриваемый метод решения СЛАУ сводится к реализации следующих пяти шагов.

1. Для простоты изложения изначально предполагается, что существует процедура устойчивых преобразований уравнений системы, расписанных по шаблонам с центральными узлами на линиях г = 1 и г = 2 (при этом 1 < ] < т), такая, что позволяет получить уравнения, соответствующие шаблонам с центральными узлами на линии г = 2. Причем структура этих новых шаблонов идентична структуре граничных шаблонов исходных уравнений на линии г = 1 (рис. 2). В этом случае преобразованная система может быть разделена на две подсистемы: первая подсистема — это уравнения на линии г = 1, 1< ] < т, а вторая подсистема состоит из уравнений, соответствующих шаблонам с центральными узлами (г,^): 2 < г < п, 1< ] < т, — и эта вторая подсистема будет замкнутой, поскольку в ней присутствуют вновь выведенные "граничные" трехи четырехточечные уравнения на линии г = 2. Следовательно, она может быть решена независимо от уравнений на граничной линии г = 1 .

2. Последовательное применение п.1 алгоритма вдоль оси х (увеличение индекса г, прямой ход) путем комбинации "граничных" уравнений на линии г с "приграничными" на линии г+1 позволяет постепенно уменьшать число уравнений в последующих замкнутых подсистемах, состоящих из уравнений, соответствующих линиям г ... п.

3. В итоге при г = п — 1 в подсистеме останется 2т уравнений, соответствующих шаблонам, представленным на рис. 3. Последнее применение предполагаемой процедуры преобразований к уравнениям с шаблонами на линиях п — 1 и п позволяет получить замкнутую систему трехточечных по г] уравнений вдоль линии г = п. Трехточечная структура этой системы обусловлена естественным усечением исходных разностных уравнений на границе области х = 1 (г = п). Такая подсистема может быть решена трехточечной скалярной прогонкой независимо от остальных уравнений преобразованной СЛАУ.

1 :.:.. У1

Iсн—Ёй—¿- М1СИ—[¿и

]+111

1 ||0......0|.....Ос

и

-Й

1У1

I

м

и

о!

-й IV-'

1сП о 51 о > !<5 $!

О

..... х

тО

з *

I

.1.

I

п-1

& А;

1 —1ф-

I _:

п-1

Рис. 2. Схема полинейного рекуррентного преобразования на разностные шаблоны, подобные граничным

Рис. 3. Схема замыкания полинейного рекуррентного преобразования на правой границе расчетной области

4. Найденные компоненты подвектора решения }, ] = 1,т подставляются в уравнения, соответствующие, в общем случае, четырехточечным шаблонам на линии г = п — 1 преобразованной СЛАУ, что заново приводит к замкнутой системе трехточечных по ] уравнений, которая снова решается прогонкой.

5. Последовательное применение п. 4 алгоритма вдоль оси х (уменьшение индекса г, обратный ход) позволяет, в итоге, определить все компоненты вектора решения {Ф^}.

Нетрудно видеть, что изложенная пятишаговая идея метода кратко может быть сформулирована следующим образом: предполагается существование устойчивого преобразования исходной СЛАУ, в результате воздействия которого структура матрицы системы примет верхний блочнотреугольный вид. Размерность блоков в данном случае будет равна т. При этом диагональные блоки представляют собой подматрицы с ненулевыми элементами на главной и прилегающих к ней диагоналях. Блоки, прилегающие сверху к диагональным блокам, — это диагональные подматрицы с ненулевыми элементами на главных диагоналях. Все остальные блоки матрицы преобразованной СЛАУ — нулевые.

3. Процедура преобразования

Для вывода предполагаемого в п. 1 разд. 2 преобразования сначала рассматриваются первые три уравнения исходной СЛАУ:

«Рп Ф 11 = «ЕП Ф21 + «N11 Ф12 + Ьп; (5)

«Р12 Ф12 = «Е12 Ф22 + «N12 Ф13 + ^12 Ф11 + Ь12; (6)

«Р13 Ф13 = «Е13 Ф23 + «N13 Ф14 + «^13 Ф12 + Ь13. (7)

Комбинация выражений (5) и (6) с целью исключения Ф11 позволяет записать уравнение

ар12 Ф12 = «N12 Ф13 + «Ё12 Ф22 + «5^12 Ф21 + в12, (8)

где аР12 = «Р12 — «N11 «512 М\1 , а5Ё12 «Ец «512 /«Р11, в12 = («512/«Р11 )Ь11 + Ь12.

Такое преобразование устойчиво, поскольку деление производится на элемент главной диагонали матрицы СЛАУ, который всегда: 1) отличен от нуля; 2) по модулю больше суммы модулей остальных элементов строки матрицы или равен этой сумме [1-3]. Комбинация выражений (7) и (8) с целью исключения Ф12 позволяет записать уравнение

«Р13 Ф13 = «N13 Ф14 + «Е13 Ф23 + «5Е13 Ф22 + «55Е^ Ф21 + /?13, (9)

где аР13 = «Р13 — «N12 «513 /аР12 , а5Е13 = «Е12 «513 /аР12 , а55Е13 = а5Е12 «513 /аР12 ,

013 = («513/«Р12 )в12 + Ь13-

Поскольку затруднительно предложить прямое устойчивое преобразование, позволяющее исключить из (9) слагаемое, содержащее Ф21, не привлекая новых неизвестных, следует применить компенсационный подход. Пусть АФк31 = Ф^31 — Фк — поправка компоненты вектора решения в узле (г,^) между к и к+1 итерациями. Линейная экстраполяция АФ^1 через АФ^1 и АФ^1 для равномерной сетки имеет вид АФк+ = 0 ( 2АФк+1 — АФ^1) или

Фк^1 = Фкк1 + 0 [2 Ф2+1 + ФкЗ1 — 2 Фкк2 — ФУ . (10)

Здесь в — параметр компенсации. В силу того, что компенсация не может быть более чем полной [2], параметр в имеет естественное ограничение: 0 < в < 1. Подстановка (10) в (9) позволяет записать уравнение

архз Ф^1 = Фм+1 + «Е13 Ф2+1 + а^з Ф2+1 + Аз, (11)

где

аР13 = аР13 — а^12 а^13/аР12 , аЁ13 = аЕ13 — ва^Ё12 а^13/аР12 ,

а^е13 = (аЕ12 + 2 вasE12 )(а^13 /аР12))

Аз = &13 + [в12 + аsEl2 (Ф^ - 2в + в Ф|з)] ^13/а^)•

Как и при выводе (8), данное преобразование также является устойчивым. Вместо линейной формулы для поправки решения можно воспользоваться квадратичной, которая, в общем случае, будет точнее и соответственно может повысить скорость сходимости:

дФкз1 = в [ з (дФкз - дФ^зз1)+дФк4+1

или

Фкз1 = ф21+в [ з (Фк2+1 - Ф2+1) + Фкз1 - з (Фк2 - ф£з) - Фк4] • (12)

Подстановка (12) в (9) позволяет вместо (11) записать уравнение

аР13 Ф?21 = ^N13 Ф?4+1 + амЕ13 Ф2+1 + «Е^ Ф2+1 + аsElз Ф2+1 + Аз, (13)

где

аР13 = аР13 - а^2 /аР12 ) аМЕ13 = ваSEl2 /аР12 > аЕ13 — аЕ13 3ваSE12 аS13 / аР12, аsElз = (аЕ12 + 3ваsEl2) (aslз / «Р12),

Аз = &13 + [в12 + аsEl2 (ф2х + 3в Ф2з - зв Ф22 - в Ф*^] (aslз /«Р12 )•

Для определенности дальнейшие преобразования проводятся для варианта метода с линейной экстраполяцией поправки решения. Преобразования варианта с квадратичной экстраполяцией хотя и несколько более громоздкие, но принципиального отличия от линейного случая не имеют.

Комбинация исходного уравнения с центральной точкой в узле (1,4)

аР14 Фм1 = аЕ14 Ф221 + аМ14 Ф1+1 + а^4 Ф1+1 + &14

с уравнением (11) с целью исключения Ф1+1 и последующим применением формулы линейной экстраполяции (10) вновь приводит к уравнению типа (11). В общем случае для произвольного ] будет иметь место уравнение

«Р1. Ф1+1 = а^. Ф1++11 + аЕ1. Ф2+1 + аsElj. Ф22-1 + ву, (14)

где коэффициенты вычисляются аналогично (11).

Повторение подобной цепочки преобразований, но уже не по возрастанию индекса ] от нижней границы ] = 1, а по его убыванию от верхней границы ] = т позволяет выписать подобное (14) уравнение вида

1Р1] Ф1+1 = а^. ф1+-11 + 7Е^ Ф2+1 + 7*Ец Ф12+11 + ¿у • (15)

Рис. 4. Схемы комбинаций исходного уравнения с преобразованиями в целях получения уравнения, соответствующего "зеркальному" шаблону: а — линейная экстраполяция, б — квадратичная экстраполяция

В итоге, для одного и того же центрального узла (1, j) имеют место три уравнения: (4), (14), (15). Причем в этих уравнениях заложена информация, которая, во-первых, содержится в уравнениях с центральными узлами только на линии г =1, во-вторых, передана с границ у = 0 (уравнение (14)) и у = 1 (уравнение (15)).

Вычитание из (4) уравнений (14) и (15) приводит к уравнению вида

РРц ф17+1 = Ф+1 + РЯВц Ф + Рмвц Ф^1! + Ян, (16)

разностный шаблон которого "зеркален" шаблону (4) за исключением того, что его центральный узел (1, j), в отличие от (4), расположен "на краю" шаблона. На рис. 4 представлены шаблоны соответствующих преобразований для вариантов как с линейной экстраполяцией поправки решения (рис. 4, а), так и с квадратичной (рис. 4, б).

Комбинация (16) с уравнением (3), записанным для г = 2, с целью исключения Фу-1 приводит к уравнению вида

«Р2; Ф2+1 = ащ Ф2++1! + авь, Ф3+1 + Ф2+-1! + Ьъ, (17)

структура которого идентична структуре уравнения (4), что и требовалось выполнить согласно предположению п. 1 разд. 2 алгоритма преобразований.

Необходимо заметить, что общая идея метода решения изложена для случая, когда глобальное уменьшение числа уравнений подсистемы (прямой ход) и вычисления подвекторов решения (обратный ход) производятся вдоль оси х, а локальные преобразования уравнений — вдоль оси у. Поскольку алгоритм можно симметрично развернуть, сделав координату у глобальной, а х — локальной, под одной итерацией в дальнейшем будет пониматься последовательность двух глобальных проходов по обеим координатам. Причем приближение решения, полученное при прохождении вдоль первой координаты, используется в качестве известного для определения поправки решения при прохождении вдоль второй координаты.

Подробные расчетные формулы вычислительного алгоритма и обоснование корректности метода изложены в [6].

4. Оценка эффективности метода

Эффективность предлагаемого метода оценивается путем сравнения результатов решения задачи Дирихле разными численными методами. Задача формулируется в единичной области П уравнением (1) при условии

vx = 1 + 2

(x - 0.5)2 + (y - 0.5)2J , Vy = 1 + 2 [0.5 - (x - 0.5)* - (y - 0.5)'

В качестве решения выбирается пробная функция

u (x, y) = 256 [xy (1 - x) (1 - y)] 2, (18)

которая отличается от пробной функции из работ [1,5] w (x, y) = const x y (1 - x) (1 - y) тем, что имеет четвертую степень зависимости от координат, а не вторую. Дело в том, что для случая квадратичной зависимости решения от координат экстраполяционная формула (12) является точным выражением. В этом случае предлагаемый метод в варианте квадратичной экстраполяции поправки решения будет прямым, что и показали проверочные расчеты, когда независимо от начального приближения, сетки и требуемой точности решение "сходилось" за одну итерацию. По той же причине вариант метода с линейной экстраполяцией поправки решения (формула (10)), как и модифицированный полинейный метод [5], является прямым для случая линейного решения.

Как следствие из (18) краевое условие первого рода задачи имеет вид Фг = 0. Система линейных уравнений получается путем аппроксимации исходной задачи схемой вида (3) на равномерной сетке с одинаковым числом узлов по каждой координате. При этом вектор правой части системы рассчитывается как произведение матрицы коэффициентов системы на вектор точного решения {uj}, компоненты которого вычислены в узлах сеточного разбиения по формуле (18). Динамика сходимости решения

Rk

\2 / п г\2

отслеживается по значению отношения норм векторов невязок

Rk

/ IIR0 II, где

2 л 1/2

{е b- к+^еj+*k>+1+j+j } • (19)

Суммирование в (19) производится по всем узлам сеточного разбиения расчетной области. В качестве начального приближения берется единичный вектор Ф0 = 1. Реше-

Rk

/ ||R0|| < е, где е — заданная

ние считается наиденным, если выполняется условие точность.

Для сравнения были выбраны следующие методы: последовательной верхней релаксации (SOR) [1], полинейный (LL) [3], модифицированный полинейный метод (mLL) [5]. Поскольку эффективность используемых методов, включая предлагаемый, существенно зависит от величины параметра релаксации 0, каждый расчет каждым методом проводился при оптимальном для данного метода значении 0opt. Делалось это для того, чтобы сравнительные характеристики сходимости были получены при наилучших для каждого метода условиях (т. е. проводилась своеобразная "оценка сверху" возможностей каждого метода). Подбор 0opt производился серией пробных расчетов с целью минимизации числа итераций для достижения заданной точности.

По аналогии с аббревиатурой тестовых методов для предлагаемого метода введены сокращенные наименования: полинейный рекуррентный метод с линейной экстраполяцией — LR1, с квадратичной экстраполяцией — LR2.

Рис. 5. Поведение отношений норм невязок в зависимости от номера итерации: 1 — SOR, 2 — LL, 3 — mLL, 4 — LR1, 5 — LR2 lg

Расчеты проводились на сетке 101 х 101 при е = 5 • 10-14. Поведение отношений норм невязок в зависимости от номера итерации k представлено на рис. 5. Видно, что за первый итерационный шаг метод LR2 понижает начальную норму невязки более чем на четыре порядка, LR1 — примерно на 3 порядка, mLL — примерно на два порядка. Классические методы SOR и LL дают за первую итерацию понижение отношения нормы невязки около одного порядка. Общее минимальное число итераций, необходимых для достижения заданной точности, также различается на порядки: для SOR и LL точность тяготеет к величине 103, для mLL — к 102, а для LR1 и LR2 — к первому порядку. Такой результат однозначно демонстрирует преимущество методов, обеспечивающих передачу информации от одного расчетного узла к другому практически за одну итерацию, над методами с механизмом поточечно-полинейной передачи.

5. Анализ применения переменного итерационного параметра

В работе [7] формула вычисления переменного значения параметра компенсации при проходе вдоль x-координаты с учетом требования 0 < вj < 1 имеет вид

вгз = max j0,min I1, [aE..Фк+1+ aw..Ф^ j / [(ащ + aw..}} •

Несколько расширяя данную технологию вычисления итерационного параметра, можно принять eij = e0eVj, где во — const, параметр задачи, а в- рассчитывается по технологии [7]. При этом здесь и далее в0 £ [0, 1].

Аналогичным образом для полинейного рекуррентного метода можно построить алгоритм вычисления значения переменного параметра для каждого итерационного уравнения СЛАУ. Для этого необходимо в уравнениях типа (11) или (13), но выписанных в общем случае для узла (i,j), приравнять нулю многочлены в круглых скобках, используемых в выражениях для расчета величин ^j и 8ij. Например, для уравнения

Фгj aNij ^j+1 + aEtj Фг+1j + asEjj Фг+1-1 + Aj,

в котором ву = Ьу + ву-1 + «5Е,

^3-1

Фк+1у-2 —20 Фк+1у-1

ходимо положить Фк+1у-2 — 20 Фк+1у-1 + 0 Фк+13

+ 0 ф+и)] («5;; /«Р,;-1), необ-0, откуда следует выражение для

вычисления 0У.

Тогда, с учетом ограничения 0 < 0^ < 1, при использовании метода с линейной экстраполяцией поправки решения (ЬЯ1) в случае глобального направления вдоль оси х (индекс г) нетрудно получить выражение для вычисления переменного параметра компенсации:

— при возрастании ]

0у = 0о ■ ш1п{1,шах{0, Ф^+1 ¿-2/ (2Ф?+1 у-1 — Ф^+1

при убывании ]

0^ = 0о ■ шт {1, шах <¡0, ФгА+1 ¿+2/ (2Ф

к _Фк

¿+1¿+1 Фг+13

.

(20)

(21)

Подобным же образом для варианта метода с использованием квадратичной экстраполяции поправки решения (ЬЯ2) в общем случае получаются следующие расчетные выражения:

— при возрастании ]

0гз = 0о ■ шш{1,шах^Ф^ ¿.-2/ (а (Ф*+13-1 — Ф^+1+ Ф^+13+^} ,

при убывании ]

0у = 0о ■ ш1п{1,шах{0, Ф^+1¿+2/ (а (ф^+1 ¿+1 — Ф^+ Ф^ ¿-1)}} .

(22)

(23)

Поскольку в приграничных узлах в методе ЬЯ2 применяется линейная экстраполяция поправки решения, для определения параметра компенсации в этом случае необходимо использовать формулы типа (20) и (21).

Аналогичные формулы выписываются и для глобального направления вдоль оси у (индекс ]) с естественным учетом того, что в качестве известных значений приближения решения берутся величины, полученные при прохождении вдоль глобального направления оси х.

Для выяснения эффективности применения переменного параметра компенсации были проведены сравнительные расчеты задачи разд. 4 с точностью е = 5 ■ 10-14 на сетках 101 х 101, 201 х 201 и 401 х 401 методами шЬЬ, ЬЯ1 и ЬЯ2. При этом в половине расчетов параметр 00 полагался равным 1 (не оптимизировался), а в другой половине — подбирался исходя из требования минимизации числа итераций, необходимых для достижения требуемой точности.

В таблице представлены результаты расчетов по числу затраченных итераций. Для метода ЬЯ2 в силу того, что в некоторых случаях оптимальные результаты достигались при 00 = 1.0, эти результаты отнесены к колонкам с оптимизированным 00, а в соответствующих колонках с не оптимизированным 0о проставлены прочерки.

По результатам видно, что применение переменного параметра хорошо проявляет себя для методов шЬЬ и ЬЯ1 при отсутствии оптимизации, когда 0 либо 00 принимался равным 1 на всех расчетных сетках. Исключение составил расчет методом шЬЬ с переменным параметром на сетке 401 х 401, в котором отношение / ||Л°|| в лучшем случае достигало примерно 10-3. В то же время для метода ЬЯ2 применение переменного параметра не дает сколько-нибудь заметного эффекта как для не оптимизированного 0

Результаты расчетов

Метод Сетка Итерационный параметр

не оптимизированный оптимизированный

постоянный переменный постоянный переменный

е Число итераций 0о Число итераций е Число итераций 0о Число итераций

Модифицированный полинейный {тпЩ 101x101 1 229 1 83 0.997 36 0.9975 35

201x201 1 449 1 188 0.9992 58 0.9992 53

401x401 1 877 1 Нет сходимости 0.9998 104 0.9998 118

Полинейный рекуррентный с линейной экстраполяцией (LR1) 101x101 1 52 1 26 0.997 16 0.998 17

201x201 1 109 1 63 0.9989 24 0.9991 25

401x401 1 3363 1 387 0.9996 35 0.9996 36

Полинейный рекуррентный с квадратичной экстраполяцией (LR2) 101x101 1 21 1 18 0.99992 16 0.99999 17

201x201 1 26 - - 0.99999 20 1.0 23

401x401 - - - - 1.0 31 1.0 28

(либо в0), так и для оптимизированного. В целом для методов шЬЬ и ЬЯ1 нетрудно заметить, что для оптимизированных в (либо в0) число итераций сокращается в 2-3 раза и более и практически не зависит от того, использован переменный параметр или нет. Примечательно, что в этих случаях оптимизированные значения в и в0 практически совпадают.

а б

Рис. 6. Динамика отношений норм невязок при расчетах методом шЬЬ: а — сетка 101 х 101, б — сетка 401 х 401; параметр компенсации постоянный: 1 — 9 = 1.0, 3 — 9 = ; переменный: 2 — 0о = 1.0, 4 — 0о = 9ор

На рис. 6 представлена динамика отношения Кк / ||Д0|| в зависимости от номера итерации на различных сетках при использовании метода шЬЬ. Из поведения графиков хорошо видно, что во всем диапазоне требования по точности применение переменного параметра компенсации не исчерпывает всего потенциала ускорения расчетов, а в случае оптимизации 9 иным способом — практически не влияет на скорость сходимости. Об этом говорит заметное различие между кривыми 1 и 2, в то время как кривые 3 и 4 практически сливаются. К тому же здесь выявляется еще одна любопытная деталь применения переменного параметра: оно может приводить к понижению устойчивости вычислительного процесса на начальных итерациях либо вообще не позволяет достигнуть наперед заданной точности (рис. 6, б).

Подобные результаты на всем диапазоне точности наблюдаются и при использовании метода ЬЯ1 (рис. 7) за исключением того, что в данных расчетах не выявлено негативного влияния применения переменного параметра на устойчивость вычислений.

Рис. 7. Динамика отношений норм невязок при расчетах методом ЬЕ1: а — сетка 101 х 101, б — сетка 401 х 401; параметр компенсации постоянный: 1 — 9 = 1.0, 3 — 9 = ; переменный: 2 — 0о = 1.0, 4 — 0о = 9оРь

Рис. 8. Динамика отношений норм невязок при расчетах методом ЬЯ2: а — сетка 101 х 101, б — сетка 401 х 401; параметр компенсации постоянный: 1 — 9 = 1.0, 3 — 9 = 0,^; переменный: 2 — 9о = 1.0, 4 — 9о = 9ор

И, наконец, на рис. 8 приведены результаты расчетов методом LR2. Плотно примыкающие друг к другу кривые свидетельствуют о том, что значение параметра компенсации близко к 1, а использование переменного параметра практически не сказывается на повышении скорости сходимости.

Заключение

По результатам изложенного материала можно сделать следующие выводы.

1. Применение переменного параметра компенсации не исчерпывает всего потенциала оптимизации.

2. Не для каждого метода использование переменного параметра эффективно: метод LR 2 практически не чувствителен к нему.

3. Оптимизация параметра компенсации, как постоянного, так и переменного, дает фактически одинаковый результат по динамике сходимости методов.

4. Применение переменного параметра может заметно понижать устойчивость итерационного алгоритма, а при определенных условиях не позволяет достигнуть сходимости решения.

Список литературы

[1] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

[2] Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995.

[3] Патанкар С. Численные методы решения задач тепломассобмена и динамики жидкости. М.: Энергоатомиздат, 1984.

[4] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидродинамика и теплообмен. Т. 1. М.: Мир, 1990.

[5] Зверев В.Г. Модифицированный полинейный метод решения разностных эллиптических уравнений // ЖВМ и МФ. 1998. Т. 38, № 9. C. 1553-1562.

[6] Фомин А.А., Фомина Л.Н. Полинейный рекуррентный метод решения разностных эллиптических уравнений. Деп. в ВИНИТИ. 2007. № 385-В2007.

[7] Zverev V.G. About the iteration method for solving difference equations // Lecture notes in computer science. Berlin; Heidelberg: Springer-Verlag, 2005. Vol. 3401. P. 621-628.

Поступила в редакцию 11 ноября 2008 г.