CC BY
283
Наука к Образование
МГТУ им. Н.Э. Баумана
Сетевое научное издание
Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 08. С. 14-31.
ISSN 1994-0408
Б01: 10.7463/0815.0791351
Представлена в редакцию: Исправлена:
© МГТУ им. Н.Э. Баумана
03.08.2015 18.08.2015
УДК 519.6
Итерационный метод решения СЛАУ на основе механической аналогии
Берчун Ю. В.1'", Бурков П. В.1, Чиркова А. С.1, Прокопьева С. М.1, Рабкин Д. Л.1, Лукьянов А. А.1
1МГТУ им. Н.Э. Баумана, Москва, Россия
у b erchmi Sjnail-iu
В статье предложен новый принцип построения итерационного метода решения СЛАУ, основанный на численном интегрировании системы линейных ОДУ второго порядка. Указанная система представляет собой математическую модель динамики механической систе-мы. Рассматривается критерий останова интегрирования системы ОДУ в рамках итерации алгоритма решения СЛАУ. Предложен алгоритм аппроксимации решения системы ОДУ для сокращения вычислительной сложности метода. Показано, что полученную аппроксимацию можно эффективно использовать для однопараметрической оптимизации на каждой итера -ции метода.
Ключевые слова: СЛАУ, итерационные методы, методы установления, механическая аналогия, аппроксимация
Введение
На сегодняшний день большинство численных методов решения систем дифференциальных уравнений (как обыкновенных, так и в частных производных) сводится к решению систем линейных алгебраических уравнений (СЛАУ) вида АХ = В, где А - матрица коэффициентов, В - вектор свободных членов, а N - число неизвестных в векторе Х. Тем самым проблема эффективного (т.е. с требуемой точностью и за приемлемое время) решения СЛАУ является определяющей для возможности практического применения устойчивых схем интегрирования в программном обеспечении, предназначенном для моделирования технических систем самого различного назначения. Между тем, наличие мощного инструментария математического моделирования, опирающего на вычислительную мощность современных компьютеров, на сегодняшний день является необходимым фактором при разработке новых технологий и инновационных образцов техники.
Методы решения СЛАУ традиционно подразделяются на прямые и итерационные. На их работу оказывают влияние целый ряд факторов, однако в первую очередь говорят о
размерности системы. При программной реализации прямых методов, в теории дающих точный результат, сказываются ограничения средств вычислительной техники. Прямые методы не учитывают фактора высокой степени разреженности исходных матриц, в процессе решения могут быть получены практически полностью заполненные матрицы. Как следствие, для хранения структур данных требуется объём оперативной памяти, который, в общем случае, пропорционален квадрату размерности СЛАУ.
Для того, чтобы этого избежать, применяют перестановку строк и столбцов матрицы до запуска самого метода решения. Например, для симметричных матриц существует алгоритм Катхилла-Макки (СиШШ-МсКее), который позволяет привести исходную матрицу к ленточному виду с минимальной шириной этой ленты. Применения подобных перестановок (самих по себе ресурсоёмких, т.к. они, как правило, основаны на полном обходе графа с использованием поиска в ширину) сокращает не только требуемый объём памяти, но и процессорное время, необходимое для работы прямого метода решения СЛАУ. Таким образом, рассмотренные перестановки (перенумерация) в сочетании с увеличением доступных объёмов оперативной памяти на современных компьютерах расширяет возможности применения прямых методов с точки зрения размерности задач.
Однако это не исключает проблемы округлений при вычислениях в условиях ограниченной разрядной сетки процессора. Округления неизбежны при любых вычислениях на компьютере, но их влияние напрямую зависит не просто от общего числа арифметических операций (которое определяют при оценке вычислительной сложности алгоритма), а от числа арифметических действий, которые производятся с одним и тем же элементов данных. Особенно следует выделить операции сложении-вычитания, в ходе которых возможно получение большого количества неверных знаков мантиссы результата (например, при вычитании двух близких чисел). Для прямых методов решения СЛАУ характерным является вычисление очередных элементов матриц на основе ранее вычисленных. Например, если вычислительную сложность метода Гаусса оценивают как 0(М3), то это значит, что при вычислении правого нижнего элемента матрицы в процессе прямого хода приходится учитывать 0(М3) ошибок округления.
Для снижения влияния ошибок округления в ходе машинных вычислений прямые методы решения СЛАУ несколько модифицируются за счёт выбора главного элемента. Однако для плохо обусловленных (жёстких) систем это не является гарантией корректности решения. Кроме того, выбор главного элемента предполагает перестановку (перенумерацию) строк и столбцов матрицы, что может вступить в противоречие с логикой сокращения числа ненулевых элементов, которая обсуждалась выше.
Поэтому для решения СЛАУ большой размерности (особенно с высокой степенью разреженности) широко применяются итерационные методы. К числу их преимуществ следует отнести как возможность реализации и использования структур данных, учитывающих разреженность исходной задачи (итерационные методы как правило вообще не меняют матрицу), так и отсутствие длинных цепочек арифметических операций. Тем не менее, итерационные методы в общем случае требуют большего объёма процессорного времени. Кроме того, проблемой итерационных методов является вопрос их устойчивость.
Вместе с тем следует отметить, что вычислительная сложность итерационных методов заключается в многократном умножении матрицы на вектор. Эта операция хорошо подходит для параллельной реализации. В современных условиях она может быть ускорена с использованием векторных расширений набора команд процессора, реализована в многопоточном режиме на многоядерных центральных процессорах, вынесена для расчёта на графические процессоры. Т.е. имеются аппаратные средства, которые позволяют добиться достаточно высокой производительности даже на отдельно взятом рабочем месте.
На сегодняшний день итерационным методом, обладающим наиболее высокими показателями скорости сходимости и устойчивости, является стабилизированный метод бисопряжённых градиентов (Biconjugate gradient stabilized method, BiCGStab). Вместе с тем, на практике находят применение и более простые методы. Причём имеются перспективы их развития с учётом появления новых высокопроизводительных аппаратных средств, реализация на которых оказывает влияние на схему самих методов. Важным направлением исследований является построение предобуславливателей. Кроме того, не следует отказываться и от поисковой деятельности, в связи с чем целесообразно проанализировать те идеи, которые лежат в основе схем построения итерационных методов.
Нельзя провести чёткую классификацию итерационных методов решения СЛАУ. Например, метод наискорейшего спуска (МНС) очевидно заимствован из теории оптимизации, но его можно трактовать и как определённую модификацию для методов установления. Метод сопряжённых градиентов (МСГ, conjugate gradient, CG), с одной стороны, является развитием МНС, однако его (и другие методы, построенные на его основе) также можно отнести к методам проекционного типа (методы на основе подпространств Крылова). Вместе с тем можно классифицировать предпосылки, отталкиваясь от которых можно прийти к итерационной схеме.
Простейшие итерационные методы (метод простой итерации, метод Зейделя, метод Якоби, метод релаксации) основываются на идее выражения /-ой неизвестной из уравнения с соответствующим номером. Этот подход является чисто алгебраическим и основан на «расщеплении» исходной матрицы. В нём нет элементов математического анализа, поэтому ожидать от него качественных результатов маловероятно. Действительно, сходимость методов этой группы обеспечивается лишь при выполнении определённых условий.
Предпосылкой для создания целой группы методов является сведение задачи решения СЛАУ к решению задачи поиска минимума некоторого квадратичного функционала. В качестве такого функционала обычно выбирается квадрат нормы вектора невязок, что является характерной чертой метода наименьших квадратов (МНК). Также существует постановка оптимизационной задачи, в которой в качестве целевой функции выбирается квадрат вектора поправок (метод минимальных поправок).
Таким образом, в качестве итерационного метода решения СЛАУ в такой вариационной постановке может быть применён любой итерационный метод локальной
многопараметрической оптимизации. В одном из наиболее простых вариантов - это МНС, скорость его сходимости может быть улучшена, таким образом приходим к идее МСГ.
Следует отметить, что линейная структура задачи даёт возможность получить аналитические зависимости для рассчитываемых коэффициентов (вместо процедуры однопараметрической оптимизации в общем случае). Однако наиболее эффективные схемы методов удаётся построить лишь для задач, в которых матрица является симметричной и положительно определённой. Впрочем, такие задачи имеют большое прикладное значение. Например, при реализации метода конечных элементов всегда формируется матрица, отвечающая указанным требованиям.
Кроме того, методы матричной алгебры позволяют прийти к идее применения предобуславливателей, которые повышают скорость сходимости градиентных методов за счёт устранения «оврагов» на поверхности целевой функции, делая её линии уровня равноудалёнными (в идеале) от точки решения.
С формальной точки зрения линейной алгебры, рассмотренные выше методы могут быть описаны в терминах крыловских подпространств, образуя группу так называемых проекционных методов. В их основе лежит метод Петрова-Галёркина, который основан на выборе двух координатных систем элементов. При этом приближенное решение находится в виде линейной комбинации по одной базисной системе, а невязка ортогональна другой базисной системе.
Следует упомянуть о стохастических методах, хотя они и не имеют широкого распространения. Тем не менее, имеются примеры применения метода Монте-Карло и генетических алгоритмов для решения СЛАУ.
Ещё одну группу методов составляют методы установления, в основе которых лежит подбор такой системы однородных линейных дифференциальных уравнений первого порядка (где матрица коэффициентов очевидно связана с матрицей исходной СЛАУ), для которой вне зависимости от выбора начального приближения, решение задачи Коши в пределе приводит к устойчивой точке, которая и является решением СЛАУ. Для интегрирования используется численная конечно-разностная аппроксимация с применением явных методов.
Очевидно, что интегрирование с малым шагом хотя и приводит к результату (при выполнении определённых требований к матрице), но не способно обеспечить высокий уровень производительности. Выбор шага может быть осуществлён, исходя из минимальной невязки СЛАУ при движении в заданном направлении. Поэтому МНС можно рассматривать как частный случай метода установления с большим шагом.
Основной идеей предлагаемого метода является описание процесса перемещения точки приближённого решения в течении одной итерации с помощью системы обыкновенных линейных дифференциальных уравнений более высокого порядка (в данном случае - второго). Её интегрирование осуществляется на основе предположения о полиномиальном характере решения. При этом для выбора шага интегрирования используется принцип минимальной невязки. После окончания итерации задача решается заново с новыми начальными условиями. При этом точка постепенно приближается к точке оптимума с точки зрения задачи МНК.
1. Принципы построения предлагаемого итерационного метода решения
СЛАУ
Рассмотрим механическую аналогию, благодаря которой можно будет получить искомую систему дифференциальных уравнений. Сперва опишем её в двумерном варианте (рис. 1), далее она будет обобщена и для произвольной размерности N. Пусть имеются две бесконечные прямолинейные направляющие, расположенные на плоскости. Они пересекаются в некоторой точке. На направляющие надеты кольца, к которым присоединены пружины. Вторые концы пружин соединены с некоторым грузом. Груз без какого-либо сопротивления может двигаться в пространстве (в данном случае - по плоскости). Очевидно, что система находится в равновесии только в точке пересечения направляющих. Во всех остальных точках пространства возникают силы, которые приводят груз в движение.
Точка райнойеоия (минимум потенциальной энергии!
Рис. 1. Общая схема двумерной механической аналогии поставленной задачи
Если бы груз испытывал сопротивление среды, то в результате переходного процесса он оказался бы в точке равновесия. В условиях отсутствия сопротивления система войдёт в режим колебаний, при этом в общем случае она не пройдёт через точку равновесия (рис. 2). Важно, что равновесие в этой точке является устойчивым. Это принципиальное отличие от методов установления.
Рис. 2. Пример траектории движения груза
Для сокращения трудоёмкости последующих вычислений, нормируем исходную матрицу и вектор правых частей СЛАУ. Следует помнить, что коэффициенты в г-м уравнении системы представляют собой координаты вектора нормали для соответствующей прямой (плоскости или гиперплоскости). Поделив все коэффициенты уравнения на норму вектора нормали, получим в качестве коэффициентов направляющие косинусы этого вектора, а правая часть после нормализации получает смысл расстояния (с учётом направления) от начала координат до соответствующей прямой (плоскости, гиперплоскости). Тогда невязка уравнения приобретает смысл расстояния (с учётом направления) от анализируемой точки до соответствующей прямой (плоскости, гиперплоскости).
Согласно второму закону Ньютона, динамика системы описывается следующей системой обыкновенных линейных дифференциальных уравнений второго порядка:
2 Т1=1^ + 1¥$ = т-а = т-1х = \Хх12]■ (1)
Произведём замену переменных и получим систему обыкновенных линейных дифференциальных уравнений первого порядка (выражения приведены в матричной форме, для произвольной размерности):
й = Р, V = к ,
) т ■ (2) I х = V
Здесь:Ап - матрица СЛАУ после нормировки;
Вп - вектор свободных членов СЛАУ после нормировки;
АпТ - транспонированная матрица Ап.
Аппроксимируем дифференциальные уравнения (2) конечными разностями, после чего применим явный метод Эйлера:
У1+1 = СЛ1АпТ(Вп - Ап71) + V1
т 4 у
Х1+1 = х1+Т1- м
Для сокращения записей положим с = 1 и ш=1 и в дальнейшем записывать их не будем.
Продолжим проведение аналогий с механикой. Каждая деформированная пружина обладает потенциальной энергией, которая прямо пропорциональна квадрату расстояния от груза до направляющей. Таким образом, потенциальная энергия системы пропорциональна квадрату невязки (с точки зрения СЛАУ). Т.е. в естественном режиме получена задача МНК (линии уровня - рис. 3).
Рис. 3. Вид линий уровня потенциальной энергии
Ускорение груза определяется усилиями, создаваемыми пружинами. Они, в свою очередь, пропорциональны расстоянию от груза до направляющих. Таким образом вектор ускорения груза сонаправлен с антиградиентом в задаче МНК. X2
Начав движение, груз будет набирать скорость. Потенциальная энергия системы будет уменьшаться, но появится кинетическая энергия. Моделирование следует продолжать до тех пор, пока кинетическая энергия будет расти. Как только начнётся обратный процесс - это признак того, что система оказалась на «дне оврага» поверхности потенциальной энергии (достаточно близко к его оси). На этом итерация завершается -фактически, обнуляется скорость, и весь процесс начинается заново, но уже с точки, которая находится ближе к «оси оврага». Колебательная составляющая движения в этом случае сокращается, энергия системы в основном тратится на разгон груза в направлении точки равновесия. Итерации прекращаются, когда норма вектора невязки перестанет превышать некоторое пороговое значение (рис. 4).
Процесс является безусловно сходящимся, т.к. на каждой следующей итерации гарантированно снижается полная энергия механической системы - аналога.
Х2 у ^у у ' У \
у у У у У У \
У У у у у у _ .—. / \
У у у у у у Напр \ являющие прямые /
У У У У У У XI
4 6 8 10 12
Рис. 4. Пример итерационного приближения к точке равновесия
2. Аппроксимация системы ОДУ
Описанный в предыдущем разделе итерационный процесс с практической точки зрения имеет принципиальный существенный недостаток - это численное интегрирование с достаточно малым шагом системы ОДУ (2) для получения траектории груза. Необходимо предложить механизм, позволяющий при сравнительно малых вычислительных затратах получить оценку положения точки в конце очередной итерации.
Аппроксимируем решение системы ОДУ (2) полиномом четвёртой степени, коэффициенты которого являются векторами размерности N. Тогда выражение для радиус-вектора принимает следующий вид:
х(0 = + + с^2 + с^3 + с^4. (4)
Продифференцировав уравнения системы ОДУ (2) дважды, получим дополнительные уравнения, которые позволят в дальнейшем сформулировать начальные условия для задачи Коши:
АпТ(Вп - Апх)
¿(111)=_АПТАП£(1) . (5)
¿(Л0= -АпТАпх№
Запишем начальные условия, введя попутно вспомогательные обозначения х0, у0, а0
и ъо.
х I Г=0 — х0 хЩ =щ = 0
п=о и
¿(П) | с=о=АпТ Св п-АпГ0) = сТ0 (6)
¿ОЮI =-АпТАпхЩ =0 к=о к=о
1 х^) I = -АпТАп хI = -АпТАпа; = Ь0 < к=о п=о и и
С учётом начальных условий, находим коэффициенты полинома в выражении (4). Окончательно оно принимает вид:
X ( + (7)
Теперь, если мы сумеем оценить время £ *, через которое груз окажется в положении минимума потенциальной энергии, двигаясь по траектории, то координаты точки вычисляются путём подстановки этого времени в полученное выражение. Выражение для потенциальной энергии, очевидно, потребует выполнения дополнительных операций умножения матрицы на вектор, наиболее ресурсоёмких с вычислительной точки зрения. Однако применение механической аналогии позволяет перейти (в силу закона сохранения энергии) от задачи поиска минимума потенциальной энергии к задаче поиска максимума кинетической энергии. Для этого необходимо продифференцировать (7) и получить уравнение для скорости:
V ( г)=~а0 £ + 3. (8)
Кинетическая энергия пропорциональна квадрату нормы скорости, тогда:
I I V(0| | 2 = {у(о(о ) = ■ I2 + в(о?;ЬЭ ■ ^+±(Ьо]ъ1) ■ I6. (9)
Для поиска экстремума этой функции продифференцируем полученное выражение
по
^^ = 2 (а?;а?) ^ + 1(0^;^)-£ 3 +~6(К-Х) ■ £ 5. (10)
Это в конечном счёте приводит к уравнению:
2 (оМ + 4 (0а;-Х) ■ е2+-в(Г0Х) ■ I*4 = о . (11)
Аналогичный результат может быть получен и из уравнения ( V| г*; а| г*) = 0 , которое вытекает из условия касания траектории линии уровня потенциальной энергии в момент прохождения точки её локального минимума.
Полученное уравнение по сути является квадратным относительно £*. Рассчитаем его дискриминант и найдём ближайший к нулю корень:
Б =в6(ю;Х) 2-4| I I I Ч | К | | 2 (12)
МрйоЪ3^. (13)
л! (Ь0;Ь о) 1 1
Очевидно, выражение (13) справедливо лишь при Б > 0. С учётом допущений при формулировке (4), это условие на практике в общем случае не выполняется. При будем оценивать по упрощённой формуле:
Г = 2^ I -^ШЭ. (14)
(Ъ0-,Ъ0) 1 1
Следует отметить, что подкоренное выражение заведомо положительно, поскольку (аЪО; Ъо ) < 0. Это утверждение можно доказать от противного, опираясь на механическую аналогию. Если бы указанное скалярное произведение было неотрицательным, то вдоль направления вектора аЪО скорость продолжала бы бесконечно увеличиваться, уводя тем самым точку от положения равновесия.
Построим участок траектории до прохождения ближайшей точки, где будет достигнут максимум кинетической энергии (на рис. 5). Нижняя, более длинная кривая соответствует интегрированию системы ОДУ (2) явным методом Эйлера с заведомо малым шагом. Таким образом, эту траекторию можно считать эталонной. Верхняя (более короткая) кривая получена в результате аппроксимации решения с использованием выражения (7).
Х2
А м лпроксимания ногочленом 4 степени \ [ /
\ у/ Явньи малыг 1 метод Эйлера с I шагом XI
9.6 9.8 10
Рис. 5. Пример приближения траектории полиномиальной аппроксимацией
Траектории практически совпадают, однако оценка времени £ * по формуле (13) является более консервативной (в силу ограниченности степени аппроксимирующего полинома). В данном случае это не является недостатком, поскольку сокращает амплитуду колебаний на следующей итерации.
Предложенный метод итерационного решения СЛАУ был реализован в виде программы на языке С++ в среде разработки Visual Studio 2013. На данном этапе применялись наиболее простые структуры данных, не учитывающие особенности исходной матрицы (такие как, например, доля ненулевых элементов, симметрия, ширина профиля и др.).
Всего проведено два типа тестов. В тесте №1 проведено решение множества задач с матрицами размерности 2х2 в зависимости от начального приближения и угла между прямыми. Определено число итераций, необходимых для достижения точности в 15 верных знаков мантиссы.
В тесте №2 исследована скорость сходимости предложенного метода в сравнении с МНС и МСГ (по мере роста размерности) при решении СЛАУ с симметричной матрицей, полученной в результате математического моделирования некоторого процесса методом конечных элементов.
Для проведения теста №1 возьмём две прямые, проходящие через начало координат (рис. 6). Первая совпадает с осью абсцисс, вторая повернута на угол а, который варьируется от 1 до 90 градусов с шагом в 1 градус. Точки начального приближения выбираются на окружности единичного радиуса и параметризуются при помощи угла ß, варьируемого в интервале от 0 до 90 градусов с шагом в 1 градус.
3. Тестирование метода
R=1
Рис. 6. Параметризация исходных данных для выполнения теста №1
В результате тестирования была получена поверхность (рис. 7), характеризующая связь числа итераций, необходимых для сходимости метода с максимально возможной машинной точностью (для типа double).
Для проведения теста №2 решим методом конечных элементов следующее дифференциальное уравнение с приведенными краевыми условиями I рода:
Результаты решения данного уравнения приведены в таблице. Стоит отметить, что предложенный метод проиграл МНС и МСГ на малых размерностях, однако при увеличении сложности задачи и, соответственно, росте размерности предложенный метод становится конкурентоспособным по сравнению с методом сопряженных градиентов, оставляя МНС позади.
Количество итераций
Рис. 7. Зависимость числа итераций от параметров тестовой задачи
Число КЭ, функции формы Предложенный метод РБМ (МНС) СОМ (МСГ)
10 Лин. 30 23 9
Квадр. 185 123 19
100 Лин. 4379 2076 94
Квадр. 13708 10430 198
500 Лин. 58250 46951 474
Квадр. 19504 234852 1046
1000 Лин. 17727 179334 949
Квадр. 16985 894750 2107
2500 (лин.) 19976 1050883 2373
5000 (лин.) 19425 - 4746
10000 (лин.) 12981 - 9494
25000 (лин.) 14057 - 23735
(Прим.: последние 3 теста для МНС были прерваны до достижения заданной точности в силу ограничения по времени.)
Заключение
В работе предложен оригинальный метод итерационного решения СЛАУ. Представляет интерес его дальнейшее исследование, при этом можно выделить ряд приоритетных направлений:
• решение жёстких и сверхжёстких систем;
• возможности применения предобуславливателей;
• эффективная программная реализация с учётом специфики структуры матрицы;
• введение демпфирования в механической аналогии;
• усложнение моделей пружин в механической аналогии.
Разработчики планируют провести тестирование предложенного метода на матрицах, получаемых в процессе решения реальных прикладных задач, в первую очередь - задач механики деформируемого твердого тела. Также для дальнейшего тестирования метода будут использованы матрицы из коллекции Ма1;пхМагке1;, которые являются
эталонными для оценки качества (точности и производительности) методов решения СЛАУ.
Список литературы
1. Александров А.А., Димитриенко Ю.И. Математическое и компьютерное моделирование - основа современных инженерных наук // Математическое моделирование и численные методы. 2014. № 1. С. 3-4.
2. Маничев В.Б., Глазкова В.В., Кожевников Д.Ю., Кирьянов Д.А., Сахаров М.К. Решение систем линейных алгебраических уравнений с удвоенной точностью вычислений на языке Си // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 2011. № 4. С. 25-36.
3. Жук Д.М., Маничев В.Б., Сахаров М.К. Сравнение современных решателей жестких систем обыкновенных дифференциальных уравнений с решателями Си библиотеки SADEL // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 8. С. 283-300. Б01: 10.7463/0812.0445558
4. Трудоношин В.А. Эволюция средств моделирования динамических систем // Информационные технологии. 2012. № 10. С. 4-7.
5. Карпенко А.П. Современные алгоритмы поисковой оптимизации. Алгоритмы, вдохновленные природой: учеб. пособие. М.: Изд-во МГТУ им. Н.Э. Баумана, 2014. 446 с.
6. Карпенко А.П., Чернов С.К. Решение систем линейных алгебраических уравнений методом предобуславливания на графических процессорных устройствах // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. № 1. С. 185-214. Б01: 10.7463/0113.0525190
7. Желдаков А.В., Федорук В.Г. Параллельный алгоритм решения систем линейных алгебраических уравнений с многодиагональной матрицей коэффициентов // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. № 7. С. 257-272. Б01: 10.7463/0713.0590785
8. Матвеева Н.О., Горбаченко В.И. Решение систем линейных алгебраических уравнений на графических процессорах с использованием технологии СЦОА // Известия Пензенского государственного педагогического университета им. В.Г. Белинского. 2008. № 12. С. 115-120.
9. Саад Ю. Итерационные методы для разреженных линейных систем. В 2 т. Т. 1: пер. с англ. М.: Изд-во МГУ, 2013. 326 с.
10. Уоткинс Д. Основы матричных вычислений: пер. с англ. М.: Бином, 2009. 664 с.
11. Шевцов Г.С., Крюкова О.Г., Мызникова Б.И. Численные методы линейной алгебры. СПб.: Лань, 2011. 496 с.
12. Вержбицкий В.М. Основы численных методов. М.: Высшая школа, 2009. 848 с.
13. Тыртышников Е.Е. Матричный анализ и линейная алгебра. М.: Физматлит, 2007. 480 с.
14. Тыртышников Е.Е. Методы численного анализа. М.: Академия, 2007. 320 с.
15. Джордж А., Лю Дж. Численное решение больших разраженных систем уравнений: пер. с англ. М.: Мир, 1984. 333 с.
16. Эстребю О., Златев З. Прямые методы для разреженных матриц: пер. с англ. М.: Мир, 1987. 120 с.
Science and Education of the Bauman MSTU, 2015, no. 08, pp. 14-31.
DOI: 10.7463/0815.0791351
Received: Revised:
03.08.2015 18.08.2015
Science^Education
of the Bauman MSTU
ISS N 1994-0408 © Bauman Moscow State Technical Unversity
Mechanical Analogy-based Iterative Method for Solving a System of Linear Equations
Yu.V. Berchun1'*, P.V. Burkov1, A.S. Chirkova1, "t berchmigmail.m
S.M. Prokopieva1, D.L. Rabkin1, A.A. Lykyanov1
1Bauman Moscow State Technical University, Moscow, Russia
Keywords: SLE, iterative methods, the relaxation methods of establishing, mechanical analogy,
approximation
The paper reviews prerequisites to creating a variety of the iterative methods to solve a system of linear equations (SLE). It considers the splitting methods, variation-type methods, projection-type methods, and the methods of relaxation.
A new iterative method based on mechanical analogy (the movement without resistance of a material point, that is connected by ideal elastically-linear constraints with unending guides defined by equations of solved SLE). The mechanical system has the unique position of stable equilibrium, the coordinates of which correspond to the solution of linear algebraic equation. The model of the mechanical system is a system of ordinary differential equations of the second order, integration of which allows you to define the point trajectory. In contrast to the classical methods of relaxation the proposed method does not ensure a trajectory passage through the equilibrium position. Thus the convergence of the method is achieved through the iterative stop of a material point at the moment it passes through the next (from the beginning of the given iteration) minimum of potential energy. After that the next iteration (with changed initial coordinates) starts.
A resource-intensive process of numerical integration of differential equations in order to obtain a precise law of motion (at each iteration) is replaced by defining its approximation. The coefficients of the approximating polynomial of the fourth order are calculated from the initial conditions, including higher-order derivatives. The resulting approximation enables you to evaluate the kinetic energy of a material point to calculate approximately the moment of time to reach the maximum kinetic energy (and minimum of the potential one), i.e. the end of the iteration.
The software implementation is done. The problems with symmetric positive definite matrix, generated as a result of using finite element method, allowed us to examine a convergence rate of the proposed method. The proposed method has shown the comparable results with the most commonly used method of conjugate gradients. There are prospects for development of the proposed method, in particular through the preconditioned use.
References
1. Aleksandrov A.A., Dimitrienko Yu.I. Mathematical and computer modeling is the basis of modern engineering sciences. Matematicheskoe modelirovanie i chislennye metody, 2014, no. 1, pp. 3-4. (in Russian).
2. Manichev V.B., Glazkova V.N., Kozhevnikov D.Iu., Kir'ianov D.A., Sakharov M.K. Solving the Linear Algebraic Equation Systems with Double Precision of Computations in C Language. Vestnik MGTU im. N.E. Baumana. Ser. Priborostroenie = Herald of the Bauman Moscow State Technical University. Ser. Instrument Engineering, 2011, no. 4, pp. 25-37. (in Russian).
3. Zhuk D.M., Manichev V.B., Sakharov M.K. Comparison of modern solvers stiff systems of ordinary differential equations solvers C library SADEL. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2012, no. 8, pp. 283-300. DOI: 10.7463/0812.0445558 (in Russian).
4. Trudonoshin V. A. The Evolution of Software for Dynamic System Simulation. Informacionnye tehnologii = Information Technology, 2012, no. 10, pp. 4-7. (in Russian).
5. Karpenko A.P. Sovremennye algoritmy poiskovoi optimizatsii. Algoritmy, vdokhnovlennye prirodoi [Modern algorithms of search engine optimization. Nature-inspired optimization algorithms]. Moscow, Bauman MSTU Publ., 2014. 446 p. (in Russian).
6. Karpenko A.P., Chernov S.K. Solving systems of linear algebraic equations by preconditioning on graphics processing units. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2013, no. 1, pp. 185-214. DOI: 10.7463/0113.0525190 (in Russian).
7. Zheldakov A.V., Fedoruk V.G. Parallel algorithm for solving systems of linear algebraic equations with a multi-diagonal coefficient matrix. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2013, no. 7, pp. 257-272. DOI: 10.7463/0713.0590785 (in Russian).
8. Matveeva N.O., Gorbachenko V. I. Solution of systems of linear algebraic equations on GPUs using CUDA technology. Izvestiya Penzenskogo gosudarstvennogo pedagogicheskogo universiteta im. V.G. Belinskogo = Proceedings of Belinskii PSPU, 2008, no. 12, pp. 115-120. (in Russian).
9. Saad Y. Iterative Methods for Sparse Linear Systems. SIAM, 2003. 520 p. DOI: 10.1137/1.9780898718003 (Russ. ed.: Saad Yu. Iteratsionnye metody dlya razrezhennykh lineinykh sistemym. V2 t. T. 1. Moscow, MSU Publ., 2013. 326 p.).
10. Watkins D.S. Fundamentals of Matrix Computations. New York, John Wiley & Sons, 2002. 633 p. (Russ. ed.: Watkins D.S. Osnovy matrichnykh vychislenii. Moscow, Binom. Laboratoriia znanii Publ., 2006. 664 p.).
11. Shevtsov G.S., Kryukova O.G., Myznikova B.I. Chislennye metody lineinoi algebry [Numerical methods of linear algebra]. St. Petersburg, Lan' Publ., 2011. 496 p. (in Russian).
12. Verzhbitskii V.M. Osnovy chislennykh metodov [Fundamentals of numerical methods]. Moscow, Vysshaya shkola Publ., 2009. 848 p. (in Russian).
13. Tyrtyshnikov E.E. Matrichnyi analiz i lineinaia algebra [Matrix analysis and linear algebra]. Moscow, Fizmatlit Publ., 2007. 480 p. (in Russian).
14. Tyrtyshnikov E.E. Metody chislennogo analiza [Methods of numerical analysis]. Moscow, Akademiya Publ., 2007. 320 p. (in Russian).
15. George J.A., Liu J.W.-H. Computer Solution of Large Sparse Positive Definite Systems. Prentice-Hall, 1981. (Russ. ed.: George J.A., Liu J.W.-H. Chislennoe reshenie bol'shikh razrezhennykh sistem uravnenii. Moscow, Mir Publ., 1984. 333 p.).
16. Osterby O., Zlatev Z. Direct Methods for Sparse Matrices. Springer-Verlag, 1983. (Russ. ed.: Osterby O., Zlatev Z. Pryamye metody dlya razrezhennykh matrits. Moscow, Mir Publ., 1987. 120 p.).
