CC BY
115
УДК 519.612.2
Анализ эффективности итерационных методов решения систем линейных алгебраических уравнений
© И.К. Марчевский, В.В. Пузикова МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Для выбора оптимального в смысле вычислительной эффективности итерационного метода решения систем линейных алгебраических уравнений, возникающих при дискретизации дифференциальных уравнений в частных производных, помимо скорости сходимости следует учитывать такие характеристики системы и метода, как число обусловленности, коэффициент сглаживания, показатель «затратности». Последние две характеристики вычисляют по коэффициентам усиления гармоник, которые позволяют судить о сглаживающих свойствах итерационного метода и его «затратности», т. е. о том, насколько хуже метод подавляет низкочастотные компоненты ошибки по сравнению с высокочастотными. Предложен способ определения коэффициентов усиления гармоник, основанный на использовании дискретного преобразования Фурье. В качестве примера приведён анализ эффективности метода В1СО&аЪ с 1Ьи и многосеточным предобусловли-ванием при решении разностных аналогов уравнений Гельмгольца и Пуассона.
Ключевые слова: разреженные линейные системы, предобусловливание, сглажива-тели, дискретное преобразование Фурье, многосеточные методы.
Введение. Основной объём вычислительной работы при численном моделировании технических систем, физических явлений и технологических процессов, как правило, приходится на решение систем линейных алгебраических уравнений (СЛАУ), возникающих при дискретизации соответствующих дифференциальных или интегродиффе-ренциальных уравнений [1].
С практической точки зрения наиболее важной характеристикой итерационного метода решения СЛАУ является время, затрачиваемое на решение с заданной точностью конкретной задачи на конкретной ЭВМ, которое определяет вычислительную эффективность метода применительно к данной задаче. В то же время при теоретическом исследовании итерационных методов основное внимание уделяется скорости сходимости, т. е. скорости уменьшения ошибки, а также кластеризации собственных значений матрицы перехода от итерации к итерации. Однако опыт применения различных итерационных методов при решении задач механики сплошной среды показывает, что прямой связи между скоростью сходимости метода и его вычислительной эффективностью для конкретных задач может не быть. Это связано с тем, что теоретические оценки строятся для наиболее общего случая и не учитывают специфику решаемой задачи, а также в них не учитывается вычислительная трудоемкость каждой итерации. Поэтому для выбора оптимального в смысле вычислительной эффектив-
ности итерационного метода помимо скорости сходимости следует также учитывать и другие характеристики: число обусловленности, коэффициент сглаживания, показатель «затратности».
Постановка задачи. Рассмотрим систему линейных алгебраических уравнений
Ax = b, x, b e RNd, A e M(R)Nd xNd, det A * 0, (1)
возникающую при дискретизации уравнения Lu = f, где L — дифференциальный либо интегродифференциальный оператор, возникающий при решении задач математической физики; u — искомая функция; f — известная функция. Будем считать, что невырожденная матрица A системы является разреженной и не обладает специальными свойствами (симметрией, положительной определенностью).
Целью настоящей работы является построение методики анализа вычислительной эффективности различных итерационных методов решения задачи (1) и априорного выбора метода решения данной задачи, обладающего высокой вычислительной эффективностью. Подчеркнём, что при таком анализе будут учитываться лишь свойства матрицы системы и матрицы перехода от итерации к итерации, в то время как технические детали реализации алгоритма (распараллеливание операций, эффективность использования кэш-памяти и т.д.) в рамках данной работы не рассматриваются. Поэтому вполне возможно, что решение систем, аналогичных рассматриваемым, итерационными либо прямыми методами, реализованными в известных библиотеках (PETSc [2], HYPRE [3], INTEL MKL [4] и др.), будет давать результаты, отличающиеся от представленных в работе.
Анализ скорости изменения итерационной ошибки. Для решения системы (1) применим итерационный метод, записанный в форме метода, основанного на расщеплении:
Mx"+1 = Nx" + b, M - N = A. (2)
Здесь x" — n -е итерационное приближение к искомому решению x .
Пусть r" = b - Ax" — вектор невязки, z" = x - x" — итерационная ошибка, pn = xn+1 - xn — вектор коррекции. В дальнейших рассуждениях будем пользоваться тем, что
z n+1= M-1Nz ", (3)
и если собственные векторы yk матрицы M N перехода от итера-
Nd
ции к итерации образуют базис в RNd , то z0 = yk, где ah — не-
k=1
которые константы, соответственно
а _1 _1
х1 = М_1Nz0 = ^Х,М N, •.., жп = М_1Nz= ^ак(Х,М к)"^. (4)
к=1 к=1
Обычно [5, 6] для анализа свойств итерационного метода исследуются коэффициенты усиления отдельных гармоник, возникающих из решения разностного аналога одномерной спектральной задачи с однородными граничными условиями на равномерной сетке. Собственные числа возникающей при этом матрицы являются коэффициентами усиления соответствующих гармоник, поскольку из (4) видно, что за одну итерацию вклад к -й гармоники в ошибку изменяется в
N раз, а собственные векторы естественным образом делятся на низкочастотные гармоники, соответствующие номерам к = 1, [0,5Ыа ], и высокочастотные — при к = [0,5Л"а ] +1, ; здесь [а] — целая
часть числа ае Я .
Это обстоятельство является чрезвычайно важным для анализа эффективности итерационного метода: низкочастотные гармоники являются длинноволновыми (длина волны Х > 4И, И — величина шага по пространству), а высокочастотные — коротковолновыми (2И <Х < 4И ). Наличие длинно- и коротковолновых компонент ошибки обеспечивает жесткость задачи (3), заключающуюся в различной скорости их убывания.
Такие итерационные методы, как метод Якоби, Гаусса - Зейделя и другие редко используются для решения жёстких задач, поскольку за первые несколько итераций они практически полностью устраняют высокочастотные компоненты ошибки, тогда как низкочастотные подавляются очень медленно. Таким образом, указанные методы действуют как фильтры низких частот, поэтому их называют сглажива-телями, а максимальный коэффициент усиления высокочастотных гармоник — коэффициентом сглаживания итерационного метода р5т .
На этом эффекте основаны многосеточные методы [5-9]: на подробной сетке с характерным линейным размером ячейки И в результате выполнения нескольких итераций сглаживателя подавляются высокочастотные компоненты ошибки, а затем осуществляется переход на более грубую сетку с размером ячейки 2И . Высокочастотные (коротковолновые) компоненты ошибки исходной задачи на этой сетке станут неразрешимыми, а половина низкочастотных (на исходной сетке) компонент — высокочастотными на грубой сетке. На новом сеточном уровне снова выполняется несколько итераций сглаживателя для подавления появившихся высокочастотных компонент ошибки, после чего снова осуществляется переход на следующий сеточный уровень. На самой грубой сетке задача решается либо прямым методом, либо итерационным с приемлемой точностью. Такой подход по-
зволяет понизить «затратность» метода: на каждом сеточном уровне подавляются только те компоненты ошибки, которые являются на нем высокочастотными (коротковолновыми).
Вычисление коэффициентов усиления гармоник. Поскольку скорость убывания компоненты ошибки определяется коэффициентом усиления соответствующей гармоники, для анализа эффективности численного метода при решении конкретной задачи необходимо уметь вычислять эти коэффициенты. В общем случае собственные векторы матрицы М перехода от итерации к итерации не являются гармониками в указанном выше смысле, поэтому для определения коэффициентов усиления гармоник в данной работе предлагается пользоваться следующей схемой.
Для того чтобы перейти в (4) от собственных векторов к гармоникам используем дискретное преобразование Фурье (ДПФ) [10]. Поскольку для х е КМя его ДПФ X е См* является сопряжённо-симметричным, а для аппроксимации исходной дифференциальной задачи (линеаризованной на шаге) используется разностная схема, точная
. ^
для решения, равного константе (т. е. = £хк = 0), обратное преоб-
к=1
разование Фурье (ОПФ) будет иметь следующий вид:
Хк = £ 2|Хп+1^ Ыкм 1)П + Arg(Xи+l)], к = 1, N. (5)
V )
N п=1
Для анализа изменения итерационной ошибки (4) необходимо исследовать коэффициенты усиления первых гармоник (гармоники
с более высокими номерами на сетке с ячейками неразрешимы). Из (5) следует, что для получения коэффициентов усиления для такого числа гармоник должно выполняться соотношение М, = 2МЛ +1,
т. е. для ук е необходимо получить ДПФ хРк е С2м<1 +1.
Заметим, что любая сеточная функция х е может быть интерпретирована как дискретный сигнал, а её ДПФ X е См* — как дис-
м, -1
кретные отсчёты спектральной функции 8(ш) = £ хк+1в-гшк этого
к=0
сигнала (т. е. спектральные отсчёты), соответствующие частотам 2кп
шп =
м
я
Хп+1= (Шп), п = 0, М, -1. (6)
Из определения спектральной функции следует, что если добавить к набору отсчётов дискретного сигнала х е некоторое коли-
чество нулей, Я (ш) не изменится, а ДПФ полученного при этом набора отсчётов х = {.% ..., , 0,.., 0} даст большее число спектральных отсчётов, соответствующих частотам, более тесно расположенным в интервале [0; 2к), т. е. при дополнении сигнала нулями спектральное разрешение ДПФ повышается (рис. 1).
Рис. 1. Повышение спектрального разрешения ДПФ при дополнении
сигнала нулями:
а — исходный сигнал и модуль его ДПФ; б — сигнал, дополненный нулями,
и модуль его ДПФ
Именно этим приёмом воспользуемся для получения ХРк е С2+1.
Пусть ар =008^22Р(-/"11) + Фр] , 3 = 1^ , р = 17^, фр = V 2 + 1
= Лг§(Хр) е Я. Тогда согласно (5), собственные векторы матрицы
перехода от итерации к итерации можно представить в виде следующей линейной комбинации гармоник:
=-1-2 I *Р+11 ар, к = 1, Ыа. (7)
у +0,5р+и ' ' 1
N<1
В качестве модельного примера рассмотрим ситуацию, когда
N< к
в (4) Х0 = 2 у , т. е. итерационная ошибка есть сумма всех собствен-
к=1
ных векторов матрицы перехода от итерации к итерации. Поскольку
N — коэффициент усиления к -го собственного вектора, суммарное усиление с учётом (7) есть
Nd _1 Nd _1 1 N1 :хмN=:хмN —Ц-:\*р+1\пр= к=1 к=1 ^+и,5 р=1
—-—:
+ 0,5 р=1
М N | ш к х к \ ^ р+1
пр, (8)
л к=1 ^
т. е. коэффициентом усиления р -й гармоники является сумма
Я = V х м lN \ Ч>к \ Яр к \ ^ р+1 \.
к=1
Наиболее информативными с точки зрения анализа вычислительной эффективности итерационного метода являются максимальный коэффициент усиления высокочастотных гармоник (коэффициент сглаживания) рзт и максимальный коэффициент усиления низкочастотных гармоник ятах, а также их отношение Р = тах , которое нар
зовём показателем «затратности» итерационного метода для рассматриваемой задачи. Это отношение показывает, во сколько раз медленнее затухают низкочастотные гармоники по сравнению с высокочастотными, т. е. насколько плохо метод преодолевает жёсткость задачи (3). Чем ближе Р к единице и чем ближе к нулю спектральный
радиус матрицы итераций р(М, тем лучше метод подходит для решения данной задачи. Как показывают вычислительные эксперименты, два различных итерационных метода, обладающие близкими значениями спектрального радиуса, но имеющие на решаемой задаче различные значения показателя «затратности», могут существенно различаться по вычислительной эффективности.
Вычислительные эксперименты. В качестве примеров систем
рассмотрим разностные аналоги уравнений Гельмгольца ((А_к )и =
к е Я\{0}) и Пуассона (Ар = /2), возникающих при решении хорошо
известной модельной задачи о расчёте течения вязкой несжимаемой жидкости в каверне. Данное течение описывается уравнениями Навье -Стокса; переход к следующему временному слою осуществляется в несколько этапов: сначала решается система линейных алгебраических уравнений — разностный аналог уравнения Гельмгольца для прогноза поля скоростей, а затем — уравнения Пуассона для поправки давления.
Одновременное рассмотрение двух линейных систем связано с тем, что их матрицы обладают различными характеристиками с
точки зрения численного решения итерационными методами. На примере этих задач можно увидеть, как различия в свойствах дискрети-зируемого оператора отражаются на величине введённого показателя «затратности». Кроме того, их сравнение показывает, что упомянутые ранее многосеточные методы, несмотря на свою высокую эффективность, не являются универсальными, и в некоторых случаях (например, для разностного аналога уравнения Гельмгольца) лучше работают более «простые» методы.
Для решения систем использовался метод BiCGStab (метод би-сопряжённых градиентов со стабилизацией), предложенный Х. ван дер Ворстом в 1992 г. [9, 11, 12]. Преимуществом данного метода является высокая скорость сходимости и возможность его использования для решения систем достаточно общего вида: симметрия и положительная определенность матрицы не имеют принципиального значения. Поскольку итерационные коэффициенты ап, Рп и юп вычисляются из условия минимума нормы невязки на итерации
(
N
гп )2 ), матрица перехода, а значит, и коэффициенты уси-
]=1
ления гармоник будут меняться от итерации к итерации. С одной стороны, это делает метод весьма эффективным, поскольку гармоники, имеющие на п -й итерации довольно высокие коэффициенты усиления, на (п +1) -й могут иметь, наоборот, низкие коэффициенты усиления. Таким образом исключается ситуация, когда вклад высокочастотных гармоник уменьшается быстро, а вклад низкочастотных — медленно. С другой стороны, на отдельных итерациях вклад некоторых гармоник может не уменьшаться, а увеличиваться (рис. 2), поэтому при решении систем методом BiCGStab наблюдается немонотонная зависимость нормы невязки от номера итерации (рис. 3).
Рис. 2. Коэффициенты усиления гармоник для метода BiCGStab на разных итерациях при решении разностного аналога уравнения Пуассона (=16)
Повысить скорость сходимости метода можно путем применения различных предобуславливателей. Их влияние будет выражаться в умножении матрицы перехода от итерации к итерации метода BiCGStab
на матрицу перехода предобуславливателя М, N = М - А. Таким
образом, анализируя характеристики матрицы М ^, можно судить о повышении скорости сходимости итерационного процесса.
140 п
Рис. 3. Зависимость нормы невязки от номера итерации при решении тестовой СЛАУ методом ЫС081аЬ ( = 654481)
Одним из простейших и часто используемых методов предобу-словливания является неполное Ьи-разложение матрицы, тогда М = Ьи « А . Проанализируем, насколько эффективным является использование такого подхода для решения рассматриваемых систем.
Представляющие интерес характеристики систем и итерационного метода для различных сеток с ячейками приведены в табл. 1 и 2
(МА = тах
X к
/ тт
к
X к
— оценка числа обусловленности).
Как видно из табл. 1 и 2, показатель «затратности» Р практически не зависит от , в то время как остальные характеристики зависят от размерности задачи весьма существенно. Это очень важно, поскольку
г-1*
при больших вычисление ХМ* ^ , ^к е К" и | ф |е К "" ~ превращается в очень громоздкую задачу, поэтому для рассматриваемой системы (1) и метода (2) важно иметь характеристику, которая определяется видом дискретизируемого оператора Ь и мало зависит от , т. е. от величины шага по пространству.
Для разностного аналога уравнения Гельмгольца рассматриваемый метод весьма эффективен: для достижения необходимой точности достаточно выполнения 1-2 итераций. Это согласуется с вышеизложенными предположениями о влиянии спектрального радиуса матрицы перехода и показателя «затратности» на скорость сходимости:
для рассмотренных случаев р(М) не превышает 2 -10-6, а Р не превышает 1,15.
фк |е к 2 Ыа~
Таблица 1
Числовые характеристики матрицы перехода 1Ьи-предобуславливателя при решении разностного аналога уравнения Гельмгольца
12 56 240
мА 1,00 1,00 1,00
р(М _:]Ч) 1,5 -10 _9 2,8 -10_8 1,7 -10_6
Р зт 1,7-10 9 7,3 -10_8 1,5 -10 5
Р 1,01 1,14 1,09
Таблица 2
Числовые характеристики матрицы перехода 1Ьи-предобуславливателя при решении разностного аналога уравнения Пуассона
16 64 256
мА 85,88 589,89 3292,74
р(М 0,82 0,97 0,99
р.зт 0,13 0,24 0,30
Р 3,34 3,25 3,32
Численное решение разностного аналога уравнения Пуассона требует 30-50 итераций (при порядка нескольких сотен), при этом спектральный радиус матрицы перехода приближается к единице при измельчении сетки (т. е. с увеличением ), а показатель «затратности» Р«3,3. Следует отметить, что число обусловленности матрицы системы растёт с увеличением . Для повышения скорости решения разностного аналога уравнения Пуассона следует использовать предобусловливатель, позволяющий более эффективно подавлять низкочастотные компоненты ошибки по сравнению с ГЬЦ-предобуславливателем. Такому требованию отвечает многосеточный (МО) предобуславливатель. Многосеточное предобусловливание означает, что вместо построения матрицы М и решения на каждой итерации метода Б1СОБ1аЬ системы вида МуП = рП находится прибли-
к П П *-» *-» А
женное решение системы Ау = р с исходной матрицей А многосеточным методом в модификации, описанной в [13].
В указанной модификации есть три параметра, значения которых можно варьировать. Это параметр релаксации ю в используемом
в качестве сглаживателя методе ADLJ (Alternating Damped Line Jacobi), число предсглаживаний npre и число постсглаживаний npost.
В [14] указано, что для подобного подхода оптимальным выбором будет ш = 6 - 2л/7 « 0,7085, npre =0 и npost = 1 или npost =2 . Однако рассчитав спектральный радиус p(M-1N) матрицы перехода [6] и показатель «затратности» Р (рис. 4), можно сделать вывод, что для систем, полученных при дискретизации уравнения Пуассона, в двумерном случае наилучшим выбором является ш = 0,8559, npre = 1 и
= 1. При этих параметрах наблюдается наименьший спектраль-
n
post
ный радиус при близком к единице показателе «затратности»: р(М"^) * 0,16, р^ * 0,11, р* 1,25 при Nd = 16 и р(М"^) * 0,25, р^ * 0,15, р* 1,16 при Nd = 64.
p(M_1N)
o(M_1N)
Рис. 4. Зависимость характеристик многосеточного предобусловливания от параметра релаксации ю в АБЫ-сглаживателе при решении
уравнения Пуассона для Nd =16 (слева) и Nd =64 (справа):
1 — n =0,n
pre ' post
4 — n = 0 ,n t
pre ' post
- 1 ; 2 — n =0 , n
pre post
3 ; 5 — n =1, n
pre post
2 ; 3 — n
1, n
post
1
2 ; 6 — и _„ = 2 , n
post
7 — n =0, n t =5; 8 — n =0, n t =4
pre post pre post
В табл. 3 приведено число итераций, необходимых для решения уравнения Пуассона в тестовой задаче при Nd = 17 760 и ограниче-
нии на норму невязки s = 10 6 на первом шаге по времени. Значения
в двух нижних строках приведены для сравнения. Необходимость рассмотрения при анализе вычислительной эффективности итерационного метода помимо спектрального радиуса матрицы перехода от итерации к итерации показателя «затратности» видна из второй и четвёртой строк таблицы: при этих параметрах р(М у многосеточных методов довольно близки, а Р различаются почти в 2 раза; в итоге методу с большим значением показателя «затратности» требуется почти в 3 раза больше итераций для получения решения задачи с заданной точностью.
Таблица 3
Число итераций при решении разностного аналога уравнения Пуассона методом Б1С081аЬ с различными предобуславливателями
Предобуславливатель ю п рге п » р(М Р Число итераций
1Ш - - - - - 158
Мв 1,0000 1 2 0,21 5,04 29
Мв 0,8559 1 1 0,16 1,25 9
Мв [14] 0,7085 0 2 0,22 2,53 11
Мв 0,1000 0 2 0,77 1,01 23
На рис. 5 в логарифмическом масштабе показаны значения собственных чисел матрицы перехода от итерации к итерации для методов, перечисленных в табл. 3. Видно, что кластеризация собственных чисел имеет место только для метода, соответствующего ю = 0,1, что вместе с близостью показателя «затратности» Р к единице и объясняет достаточно малое число итераций, несмотря на сравнительно большой спектральный радиус. Для остальных методов явно выраженной кластеризации собственных чисел не наблюдается; это подтверждает необходимость анализа величины показателя «затратности», что наглядно проявляется при сравнении методов, соответствующих ю = 1,0 и ю = 0,7085.
Следует отметить, что трудоемкость одной итерации метода Б1СОБ1аЬ с многосеточным предобусловливанием в несколько раз выше, чем при использовании ГЬЦ-предобусловливания, тем не менее при использовании многосеточного предобусловливания с оптимальными параметрами общее время решения линейной системы удается сократить примерно в 4 раза.
Приведённые выше примеры решения задач относительно малой размерности (до 100 000), которые возникают при решении задач математической физики на достаточно грубых сетках, следует рассматривать как источник информации о коэффициентах усиления гармоник и показателе «затратности».
Рис. 5. Собственные числа матриц перехода для многосеточных методов
Рис. 6. Структура временных затрат при решении тестовой задачи расчёта течения в каверне методом БЮв81аЬ с различными предобусловливателями (ОрепРОЛМ, Ыа = 10000 ):
а — время счёта 243 с; б — время счёта 76 с; в — время счёта 68 с
Аналогичные результаты получаются и при решении данной задачи средствами пакета ОрепРОЛМ [15] на сетке 100 х100: наименьшее время счёта также получается при использовании ГЬЦ-предобу-словливания для решения разностного аналога уравнения Гельмголь-ца и многосеточного предобусловливания для решения разностного аналога уравнения Пуассона. На рис. 6 приведена структура временных затрат при решении тестовой задачи расчёта течения в каверне: серым обозначена доля затрат на решение уравнения Гельмгольца, белым — доля затрат на решение уравнения Пуассона, чёрным — прочие затраты. В результате выбора итерационного метода и параметров многосеточного предобусловливателя для данной задачи время счёта удалось сократить почти в 4 раза по сравнению с использованием методов решения, приведенных в ОрепРОЛМ в соответст-
вующем учебном примере. При этом важно, что параметры, подобранные в результате анализа итерационной ошибки на грубой сетке, остаются весьма эффективными и при решении той же задачи на более подробной сетке.
Таким образом, можно предложить эффективную методику выбора итерационного метода решения линейных систем, возникающих при дискретизации некоторой задачи:
1) необходимо построить разностный аналог задачи для достаточно грубой сетки;
2) для этой задачи строится матрица перехода от итерации к итерации;
3) вычисляются собственные числа и собственные векторы матрицы перехода;
4) с помощью вышеописанной схемы вычисляются коэффициенты усиления гармоник;
5) оценивается близость спектрального радиуса матрицы перехода к нулю и показателя «затратности» к единице;
6) если р(Ми р оказываются достаточно велики, следует выбрать другой метод (предобусловливатель) и вернуться к п. 2.
Заключение. Определение коэффициентов усиления гармоник позволяет получить некоторые характеристики итерационного метода, по которым можно судить о его эффективности при решении той или иной задачи. Для определения этих коэффициентов удобно использовать дискретное преобразование Фурье.
Расчеты показывают, что для итерационного метода решения СЛАУ отношение наибольшего коэффициента усиления низкочастотной гармонии к наибольшему коэффициенту усиления высокочастотной гармоники (показатель «затратности») определяется видом дискретизируемого дифференциального оператора Ь и мало зависит от Nd, т. е. от густоты сетки. Это позволяет предложить эффективную методику выбора параметров используемого итерационного метода и предобусловливателя на основе анализа величины спектрального радиуса матрицы и показателя «затратности», который может быть выполнен на достаточно грубых сетках.
На примере расчета течения вязкой несжимаемой среды в квадратной каверне показано, что для решения разностного аналога уравнения Гельмгольца для прогноза скорости целесообразно использовать итерационный метод ЫС081аЬ с ГЬЦ-предобусловливанием, в то время как для разностного аналога уравнения Пуассона более эффективным оказывается многосеточное предобусловливание с параметрами, выбранными в соответствии с предложенной методикой.
ЛИТЕРАТУРА
[1] Зарубин В.С., Кувыркин Г.Н. Особенности математического моделирования технических устройств. Математическое моделирование и численные методы, 2014, № 1, с. 5-17.
[2] PETSc. URL: http://www.mcs.anl.gov/petsc
[3] Scalable linear solvers: HYPRE. URL: http://computation.llnl.gov/casc/linear_solvers/sls_hypre.html
[4] Intel Math Kernel Library 11.0. URL: http://software.intel.com/en-us/intel-mkl
[5] Wesseling P. An introduction to multigrid methods. Chichester, John Willey & Sons Ltd., 1991, 284 p.
[6] Ольшанский М.А. Лекции и упражнения по многосеточным методам. Москва, Изд-во ЦПИ при механико-математическом факультете МГУ им. М.В. Ломоносова, 2003, 163 с.
[7] Федоренко Р.П. Введение в вычислительную физику. Москва, Изд-во Моск. физ.-техн. ин-та, 1994, 528 с.
[8] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 590 с.
[9] Saad Y. Iterative Methods for Sparse Linear Systems. New York, PWS Publ., 1996, 547 p.
[10] Сергиенко А.Б. Цифровая обработка сигналов. Санкт-Петербург, Питер, 2002, 608 с.
[11] Van der Vorst H.A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for solution of non-symmetric linear systems. SIAM J. Sci. Stat. Comp., 1992, no. 2, pp. 631-644.
[12] Ильин В.П. Методы бисопряжённых направлений в подпространствах Крылова. Сибирский журнал индустриальной математики, 2008, т. IX, № 4, с. 47-60.
[13] Пузикова В. В. Решение систем линейных алгебраических уравнений методом BiCGStab с предобуслoвливанием. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2011, спец. вып. «Прикладная математика», с. 124-133.
[14] Van Kan J., Vuik C., Wesseling P. Fast pressure calculation for 2D and 3D time dependent incompressible flow. Numer. Linear Algebra Appl., 2000, no. 7, pp. 429-447.
[15] OpenFOAM. URL: http://www.openfoam.com
Статья поступила в редакцию 01.12.2014
Ссылку на эту статью просим оформлять следующим образом:
Марчевский И.К., Пузикова В.В. Анализ эффективности итерационных методов решения систем линейных алгебраических уравнений. Математическое моделирование и численные методы, 2014, № 4, с. 37-52.
Марчевский Илья Константинович родился в 1983 г., окончил МГТУ им. Н.Э. Баумана в 2005 г. Канд. физ.-мат. наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. Автор более 30 научных публикаций в области исследования движения и устойчивости конструкций в потоке среды, вычислительной гидродинамики, высокопроизводительных вычислений, элементарной математики. e-mail: [email protected]
Пузикова Валерия Валентиновна родилась в 1992 г., окончила МГТУ им. Н.Э. Баумана в 2014 г. Ассистент, аспирантка кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. Автор 5 научных публикаций в области вычислительной гидродинамики и численных методов линейной алгебры. e-mail: [email protected]
Performance analysis of iterative methods of combined linear algebraic equations solution
© I.K. Marchevsky, V.V. Puzikova Bauman Moscow State Technical University, Moscow, 105005, Russia
When sampling partial differential equations one has to solve a system of linear algebraic equations. To select the optimal in the sense of the computational efficiency of iterative method for solving such equations, in addition to the rate of convergence we should take into account such characteristics of the system and method, as the condition number, the smoothing factor, the indicator "costs on." The last two characteristics are calculated by the coefficients of harmonics amplification that give evidence of the smoothing properties of the iterative method and its "costs on", i. e. how worse the method suppresses frequency components of the error as compared with the high-frequency ones. The suggested method of determining harmonic gain factors is based on of the discrete Fourier transform. As an example, an analysis of the effectiveness of the BiCGStab method with ILU and multigrid preconditioning when solving difference analogues of the Helmholtz and Poisson equations is described.
Keywords: sparse linear systems, preconditioning, smoothers, discrete Fourier Transform, multigrid methods.
REFERENCE
[1] Zarubin V.S., Kuvyrkin G.N. Matematicheskoe modelirovanie i chislennye metody - Mathematical Modeling and Computational Methods, 2014, no. 1, pp. 5-17.
[2] PETSc - Portable, Extensible Toolkit for Scientific Computation. Available at: http://www.mcs.anl.gov/petsc.
[3] Scalable linear solvers: HYPRE. Available at: http://computation.llnl.gov/casc/linear_solvers/sls_hypre.html.
[4] Intel Math Kernel Library 11.0. Available at: http://software.intel.com/en-us/intel-mkl.
[5] Wesseling P. An introduction to multigrid methods. Chichester, John Willey & Sons Ltd., 1991, 284 p.
[6] Ol'shansky M.A. Lektsii i uprazhneniya po mnogosetochnym metodam [Lectures and Exercises on Multigrid Methods]. Moscow, Lomonosov MSU Publ., 2003, 163 p.
[7] Fedorenko R.P. Vvedenie v vychislitel'nuyu fiziku [Introduction to Computational Physics]. Moscow, Moscow Institute of Physics and Technology Publ., 1994, 528 p.
[8] Galanin M.P., Savenkov E.B. Metody chislennogo analiza matematicheskikh modeley [Methods of Numerical Analysis of Mathematical Models]. Moscow, BMSTU Publ., 2010, 590 p.
[9] Saad Y. Iterative Methods for Sparse Linear Systems. New-York, PWS Publ., 1996, 547 p.
M.K. MapneecKuü, B.B. nysuKoea
[10] Sergienko А.B. Tsifrovaya obrabotka signalov [Digital Signal Processing]. St. Petersburg, Piter Publ., 2002, 608 p.
[11] Van der Vorst H.A. SIAMJ. Sci. Stat. Comp, 1992, no. 2, pp. 631-644.
[12] Il'in V.P. Sibirskiy Zhurnal Industrial'noy Matematiki - Siberian Journal of Industrial Mathematics, 2008, vol. 9, no. 4, pp. 47-60.
[13] Puzikova V.V. VestnikMGTU im. N.E. Baumana. Seriya: Estestvennye nauki -Herald of the Bauman Moscow State Technical University, Series: Natural Science, 2011, special issue "Applied Mathematics", pp. 124-133.
[14] Van Kan J., Vuik C., Wesseling P. Numer. Linear Algebra Appl, 2000, no. 7, pp. 429-447.
[15] OpenFOAM - The Open Source Computational Fluid Dynamics (CFD) Toolbox. Available at: http://www.openfoam.com.
Marchevsky I.K. (b. 1983) graduated from Bauman Moscow State Technical University in 2005. Ph.D. (Phys. & Math.), assoc. professor of the Applied Mathematics Department at Bauman Moscow State Technical University. Author of over 30 publications in the field of the structures in the flow motion and stability investigation, computational fluid dynamics, high performance computations and elementary mathematics. e-mail: [email protected]
Puzikova V.V. (b. 1992) graduated from Bauman Moscow State Technical University in 2013. A postgraduate at the Applied Mathematics Department of Bauman Moscow State Technical University. Author of 5 publications in the field of computational fluid dynamics and numerical methods in linear algebra. e-mail: [email protected]
