CC BY
35
Вычислительные технологии
Том 13, Специальный выпуск 4, 2008
Алгоритм ускорения сходимости в методе коллокации
и наименьших квадратов*
В. И. Исаев
Институт теоретической и прикладной механики им. С. А. Христиаповича СО РАН, Новосибирск, Россия e-mail: [email protected]
A new version of the multistep least squares algorithm for acceleration of the convergence rate for an iterative process is proposed and implemented. The method was applied for the collocation and least squares method for solution of the Navier—Stokes equations. Advantages of the proposed version over the previous ones are shown.
Введение
Построение приближенного решения краевых задач для уравнений математической физики, имеющего необходимую точность, часто связано со значительными затратами ресурсов ЭВМ. Во многих численных методах для нахождения решения используются итерационные процессы. Низкая скорость сходимости итераций существенно сужает возможности применения метода на практике, так как из-за большого времени расчета вычисления на достаточно подробных сетках становятся недоступными. Существуют различные подходы, позволяющие добиться ускорения сходимости итераций [1-4]. Так, например, в работах М.А. Ольшанского для этого используется геометрический многосеточный метод, включающий аддитивный многосеточный метод и классические V- и \¥-циклы [4]. В статье II. Саада [1] изложены алгоритмы ускорения, в которых используются подпространства Крылова. В работах А.Г. Слепцова [2, 3] описан многошаговый метод наименьших квадратов ускорения сходимости итерационного процесса х-7, ] = 0,1В этом методе через каждые к1 шагов к текущему приближению хп добавляется поправка, являющаяся линейной комбинацией к1 векторов невязок рассматриваемой системы уравнений г7 = х7+1 — х-7, ] = (п — к1),... , (п — 1), к1 < п. Поправка ищется из условия минимума функционала, составленного из квадратов невязок уравнений системы, эквивалентной исходной. С ростом номера итерации нормы невязок в случае сходящегося итерационного процесса стремятся к нулю. При этом в матрице переопределенной системы линейных алгебраических уравнений (СЛАУ), к решению которой сводится реализация рассматриваемого метода ускорения, могут появиться столбцы, близкие к линейно зависимым. В результате возникает необходимость решения на ЭВМ плохо обусловленной задачи линейной алгебры, что существенно снижает эффективность алгоритма. В [2, 3] затронута эта проблема, однако приемлемого
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00080-а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
способа ее решения для произвольного числа невязок к1 не предложено. А,Г, Слепцовым варианты метода реализованы па ЭВМ только для к1 < 2. Более того, в [2, 3] для решения этой плохо обусловленной СЛАУ используется метод наименьших квадратов (\ 111К). который приводит к дополнительному ухудшению обусловленности. Однако в [2, 3] получены хорошие результаты ускорения сходимости итераций в само- и в несамосопряженном случаях. Все это говорит о том, что метод обладает хорошими возможностями, но его можно улучшить,
В данной работе для решения указанной переопределенной СЛАУ предложен новый вариант ортогонального метода линейной алгебры, описанного в работах [5, 6], В новом варианте учитывается возможность наличия в матрице системы линейно зависимых или близких к линейно зависимым столбцов. Получаемое решение в случае отсутствия ошибок округления и невырожденной СЛАУ совпадает с тем, что получается в методе наименьших квадратов. При этом не происходит ухудшения обусловленности задачи линейной алгебры, которое имеет место в МНК, Применение нового варианта ортогонального метода позволило решить проблему обобщения алгоритма ускорения сходимости на произвольное число невязок,
1. Описание метода ускорения
Рассмотрим систему линейных алгебраических уравнений
Ах = Ь, (1)
решение которой ищется при помощи итерационного процесса. Здесь векторы х, Ь € Ыр; А — матрица размера р х р, р € N. Пусть задано некоторое начальное приближение хо-Итерационная последовательность строится следующим образом:
хга+1 = Тхп + ¥, п = 0,1,... (2)
Здесь Т — матрица размера р х р. Предположим, что итерационный процесс сходится. Тогда векторы х^ стремятся при п — то к решению х системы
х = Тх + ¥, (3)
п
ции будем называть вектор Гп = Тхп + ¥ — хп = х„+1 — х^, погрешностью — вектор
Уга = х — хп = х — хп+1 + гп = Уп+1 + Гп- (4)
Невязки и погрешности удовлетворяют соотношениям
Гп+1 = Т Гп, Уп+1 = т Уп. (5)
Из равенств (4) и (5) можно получить систему
(I — Т) уп = Гп, (6)
связывающую компоненты векторов уп и гп, где I — это единичная матрица размера р х р. Умножив уравнение (6) на Т-1 слева и заменив вектор Т-1гп на гп-1, получим
(Т 1 — I) Уп = Гп-1.
(7)
В многошаговом методе наименьших квадратов ускорения сходимости, предложенном в [2, 3], через каждые fc итераций к текущему приближению хп добавляется поправка
k—i
Уп* = Y^ ai rn-k+i- (8)
i = i
Уточненный таким образом вектор хп используется далее для выполнения следующих k итераций. Коэффициенты ai} i = 1, ..., (k — 1), ищутся из условия минимума квадрата нормы невязки уравнения (7), возникающей при подстановке в него приближенного вектора погрешности уп* вместо yn, т. е. из условия
(T-i — I)Уп* — Гп-1 2 ^ min . (9)
2 «I,...,afc-1
Вторая норма вектора v = {vi, i = 1, ..., m} определяется здесь равенством
m
l|v||2 = Y v?.
i=1
Задача о нахождении минимума (9) эквивалентна решению переопределенной СЛАУ:
(T-1 — 1)Уп* = Гп-1, (10)
в которой неизвестными являются коэффициенты ai, i = 1, ..., (k — 1), Система (10) может быть записана в виде
ai rп-k+i — rп-k )+ ••• + afc I rп—i — rп-2 )= —rп-i- (11)
Для решения СЛАУ (11) в данной работе предложено использовать ортогональный метод, эквивалентный МНК [6], При больших n в случае сходящегося итерационного процесса элементы матрицы системы (11) малы, так как нормы ||rn||^ —> 0, Здесь
n—
||v||^ = max |vj|. Чтобы избежать арифметических операций с числами, близкими
i=1,..., m
к машинному нулю, проводилась нормировка столбцов матрицы при помощи замены переменных:
ai = 11 rn-k+i — rra-fc+(i-1)i = 1) — 1)-
Однако таким способом не решаются проблемы, возникающие при наличии среди векторов
Zi = (rn-k+i - rn-k+(i-1)), i = 1) • • •) (k - 1),
близких к линейно зависимым, В этом случае задача линейной алгебры (11) будет плохо обусловлена. Причина плохой обусловленности заключается в использовании слишком большого количества невязок в (8), Здесь близкими к линейно зависимым называются такие векторы vb • • •, vg ( ||vj= 1 j = 1, • • • , (k — 1)), для которых
31 (1 < l < q) V i = l (1 < i < q) 3 ki max =0: v - E
l<i<q,i=l
l<i<q,i=l
<
где параметр зависимости e1 — фиксированная малая величина, В данной работе ее значение полагалось равным 10-15, Для представления вещественных чисел использовался тип long double (С++), В качестве критерия отсутствия в матрице СЛАУ (11)
близких к линейно зависимым столбцов использовалось условие на модули ведущих элементов ортогонального исключения [6]. Все они должны быть больше величины
Рассмотрим случай, когда указанный критерий не выполняется. Пусть на некотором этапе решения переопределенной СЛАУ (11) при проведении ортогонального исключения в .7-м столбце ее матрицы, где 1 < . < (к — 1), модуль ведущего элемента стал меньше величины £1, Из этого следует, что невязки г = 1, ..., близки к линейно зависимым. Несложно показать, что тогда вектор ъг1, где . < г1 < (к — 1), и векторы |ъг, г = 1, ..., (. — 1)} близки к линейно зависимым с параметром зависимости, равным С£1; где конетанта С зависит от нормы матрицы Т и от к. В этом случае не имеет смысла использовать невязки ъг, г = (к — 1) для вычисления поправки. Поло-
жим коэффициенты аг, г = (к — 1), равными нулю и перепишем с учетом этого
систему (11). Из полученной таким способом переопределенной СЛАУ найдем константы аг, г = 1, ..., (. — 1), Заметим, при этом можно воспользоваться результатами уже проведенного до (. — 1)-го этапа ортогонального исключения в системе (11). В итоге к вектору х добавляется поправка
т.е. вместо (к — 1) векторов для определения поправки (8) используются только первые (j — 1) невязок. Поправки yn* и yn** мало отличаются в случае отсутствия ошибок округления, а обусловленность задачи линейной алгебры, к решению которой сводится реализация метода ускорения сходимости, после сокращения числа невязок станет существенно лучше.
В методе ускорения вместо (2) можно использовать последовательность Хп = x^, n = 1, 2,..., которая составлена го вычисленных по формуле (2) векторов xn, взятых через каждые s шагов. Матрица оператора перехода от приближения xn к Хп+1 здесь равна Ts, Как показали численные эксперименты, этот прием в некоторых случаях позволяет добиться дополнительного ускорения сходимости,
2. Численные эксперименты
Предложенный вариант алгоритма ускорения сходимости использован в методе колло-кации и наименьших квадратов решения краевой задачи для стационарных уравнений Навье—Стокса в двумерной прямоугольной области [6, 7], Эксперименты проводились с несколькими решениями при различных числах Рейпольдса, На рис, 1 и 2 изображены графики, полученные для одного тестового решения при Re = 100, Результаты, аналогичные приведенным здесь, имели место во всех остальных экспериментах,
n
псевдопогрешности еп, равной ||xn+1 — xn где xn — n-e приближение для вектора неизвестных. Кривая, имеющая на рис, 1 большую толщину, соответствует случаю, когда используется метод ускорения сходимости, другая кривая — случаю, когда он не применяется. Видно, что использование ускорения позволило сократить в данном эксперименте число итераций, необходимых для получения приближенного решения, более чем в четыре раза. Итерационный процесс в методе КИК решения краевых задач для уравнений Навье—Стокса начиная с некоторого шага становится близким к
i =1
-1 -3 -5 -7 -9
0 1500 3000 4500 6000
Рис. 1. Зависимость псевдопогрешности от номера итерации (в = 5, к = 50)
3000 250020001500 1000
250023002100 1900 1700-
20 40 60 80 а
Рис. 2. Зависимость числа итераций от параметра к при в = 10 (а) и параметра в при к = 50 (б)
линейному, о чем свидетельствует то, что график для случая, когда ускорение не производится, в логарифмическом масштабе мало отличается от прямой. Поэтому здесь оправданно применение многошагового метода наименьших квадратов, в котором существенно используется независимость матрицы оператора перехода Т (2) от текущего приближения. Каждая "ступенька" на графике в случае использования метода ускорения соответствует однократному прибавлению поправки (8). Видно, что при этом псевдопогрешность существенно уменьшается.
На рис. 2, а приведен график зависимости от параметра к числа итераций, потребовавшихся для получения приближенного решения. Видно, что с увеличением числа невязок до 20 растет эффективность метода ускорения. Если значение параметра к > 20, то она практически не ухудшается по сравнению со случаем к = 20, оставаясь примерно на одном и том же уровне. Таким образом, поскольку в более ранних вариантах метода ускорения сходимости [2, 3] при вычислении поправки (8) использовалось не более двух невязок, предложенный здесь вариант позволяет добиться большего по сравнению с ними ускорения. На рис. 2, б изображен график зависимости числа итераций от параметра 5. Здесь пунктир — это усредненная кривая, полученная методом наименьших квадратов. Видно, что увеличение параметра 5 также может привести к дополнительному уменьшению числа итераций.
Таким образом, применение в методе КНК изложенного здесь алгоритма ускорения сходимости позволило существенно уменьшить время, необходимое для проведения расчетов. В некоторых случаях (при больших числах Рейнольдса) итерационный процесс из расходящегося становился сходящимся после использования предложенного здесь метода ускорения, т.е. происходило, как сказано в работах [2, 3], его исправление. Все это позволило использовать в численных экспериментах более подробные сетки, а так-
же расширить диапазон чисел Re, при которых можно успешно применять метод KHK для уравнений Навье—Стокса, В результате при расчете течения в каверне с движущейся верхней крышкой на сетке 640 х 640 наряду с известными угловыми вихрями обнаружены еще вихри второго порядка малости в углах каверны при Re = 1000, Характерные значения функции тока в угловых вихрях второго порядка в 103 раз меньше, чем в вихрях первого порядка, и в 105 раз меньше, чем в основном центральном вихре. При использовании сетки 640 х 640 в варианте метода KHK с аппроксимацией давления полиномами второго порядка общее количество неизвестных задачи около 6 • 106, Решение систем такой большой размерности стало возможным в данной работе во многом благодаря использованию рассмотренного здесь метода ускорения сходимости.
Заключение
Предложен и реализован новый вариант многошагового метода наименьших квадратов для ускорения сходимости итерационного процесса. Решена проблема его обобщения на произвольное количество невязок. Показано, что новый вариант метода позволяет добиться большего ускорения, чем те, что предложены ранее [2, 3], Использование алгоритма ускорения в совокупности с улучшениями метода KHK, предложенными в работе [6], позволило провести расчет течения в каверне на такой подробной сетке, как 640 х 640 Re = 1000
второго порядка малости.
Список литературы
fl] Saad Y. Krvlov subspace methods for Solving large unsymmetric linear systems // Mathematics of Computation. 1981. Vol. 37, N 155. P. 105-126.'
[2] Слепцов А.Г. Об ускорении сходимости линейных итераций. 1 // Моделирование в механике: Сб. науч. тр. /РАН Сиб. отд-ние. ВЦ; ИТПМ. 1989. Т. 3(20), № 3. С. 132-147.
[3] Слепцов А.Г. Об ускорении сходимости линейных итераций. 2 // Моделирование в механике: Сб. науч. тр. /РАН Сиб. отд-ние. ВЦ; ИТПМ. 1989. Т. 3(20), № 5. С. 132-147.
[4] Ольшанский М.А. Равномерные по параметру многосеточные и итерационные методы: Дис. ... доктора физ.-мат. наук. \!.. 2006. 285 с.
[5] Шапеев A.B. О двух алгоритмах численных методов линейной алгебры // Матер. XXXVIII Междунар. науч. студен, конф. "Студент и научно-технический прогресс": Тез. докл. Новосибирск, 2000. С. 16.
[6] Исаев В.И., Шапеев В.П., Еремин С.А. Исследование свойств метода коллокации и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье^Стокса. // Вычисл. технологии. 2007. Т. 12, № 3. С. 53-70.
[7] Сёмин Л.Г., Шапеев В.П. Метод коллокаций и наименьших квадратов для уравнений Навье^Стокса // Вычисл. технологии. 1998. Т. 3, № 3. С. 72-84.
Поступим в редакцию 20 февраля 2008 г.
