CC BY
15_ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА_
2023 • Математика. Механика. Информатика • Вып. 1(60)
Научная статья
УДК 531.01+004.94
DOI: 10.17072/1993-0550-2023-1-30-46
Итерационный метод решения систем линейных алгебраических уравнений с положительно полуопределенными матрицами системы
Владимир Николаевич Иванов
Пермский государственный национальный исследовательский университет, Пермь, Россия [email protected], https://orcid.org/0000-0002-2634-4296
Аннотация. Рассматривается задача нахождения нормального псевдорешения возмущенных систем линейных алгебраических уравнений с положительно полуопределенными матрицами системы. Данная задача возникает при численном моделировании динамики систем связанных твердых тел, на движение которых наложены дополнительные связи. Задача решается итерационным квазиньютоновским методом, основанным на SR1-формуле поиска минимума квадратичных функций. Доказаны основные свойства алгоритма. Приведен пример исследования динамики одной механической системы, на котором показана работоспособность метода. Ключевые слова: системы линейных алгебраических уравнений; итерационные методы; квазиньютоновские методы; SR1-формула; системы абсолютно твердых тел; уравнения движения; дополнительные связи; численное интегрирование
Для цитирования: Иванов В. Н. Итерационный метод решения систем линейных алгебраических уравнений с положительно полуопределенными матрицами системы // Вестник Пермского университета. Математика. Механика. Информатика. 2023. № 1(60). С. 30-46. DOI: 10.17072/1993-0550-2023-1-30-46.
Статья поступила в редакцию 06.11.2022; одобрена после рецензирования 26.01.2023; принята к публикации 15.03.2023.
Research article
An Iterative Method for Linear Algebraic Equations Systems Solving With Positive Semidefinite Matrices
Vladimir N. Ivanov
Perm State University, Perm, Russia [email protected], https://orcid.org/0000-0002-2634-4296
Abstract. The normal pseudo-solution for perturbed linear algebraic equations systems with positive sem-idefinite matrices finding problem is considered. This problem is actual for numerical simulation process of a multibody systems dynamics taking into account additional constraints on its motion. The problem is solved by the iterative quasi-Newtonian method based on the SR1-formula for quadratic functions minimum finding. The main algorithm properties are proved. An one mechanical system dynamics study example is given, which shows the method efficiency. Keywords: systems of linear algebraic equations; iterative methods; Quasi-Newton methods; SR1-formula; multibody system; equations of motion; additional constraints; numerical integration
For citation: Ivanov V. N. An Iterative Method for Linear Algebraic Equations Systems Solving With Positive Semidefinite Matrices. Bulletin of Perm University. Mathematics. Mechanics. Computer Science. 2023;1(60):30-46. (In Russ.). DOI: 10.17072/1993-0550-2023-1-30-46.
The article was submitted 06.11.2022; approved after reviewing 26.01.2023; accepted for publication 15.03.2023.
Эта работа © 2023 Иванов В. Н. лицензируется под СС BY 4.0. Чтобы просмотреть копию этой лицензии, посетите http://creativecommons.Org/licenses/by/4.0/.
Введение
Рассматривается задача определения нормальных псевдорешений возмущенных систем линейных уравнений (СЛАУ) с положительно полуопределенными матрицами системы. Задача возникает при численном моделировании динамики систем абсолютно твердых тел (САТТ), на движение которых наложены дополнительные связи.
Будем считать, что в САТТ все связи голономны, склерономны и идеальны. Введем
обозначения: q = q (I) = (^,..., ^ ) - матрица-
столбец обобщенных координат, </ = —, t -
ей
время, N - число тел в системе, V = со! | V ,..., | = Ш[ - матрица-столбец абсолютных скоростей всех тел системы, В = 5 = В ^) - матрица касательного локального базиса многообразия . Пусть на
механическую систему наложено п дополнительных связей:
е (q ) = о. (1)
Будем считать, что их избыточное число, т.е. часть из них либо является линейно зависимыми, либо тождественно равными 0. Такая ситуация возникает при образовании замкнутых кинематических цепей с помощью замыкания части вершин первичного графа механической системы со структурой дерева или, другими словами, наложения дополнительных шарниров между телами.
В этом случае уравнения движения механической системы можно выписать в обобщенных координатах в форме уравнений Ла-гранжа 1 рода [1-3], учитывая дополнительные связи с помощью множителей Лагранжа
М = (м1,...,мп):
(2)
= б 4+/1^,4) = 0,
* * / \ Т где М = М (q) = В МВ - симметричная, положительно определенная инерционная матрица, М = diag (м^,..., М^) - матрица масс и
моментов инерции тел системы, / = / (</,</,/)
- вектор слагаемых уравнений движения, не зависящих от ускорений, включающий актив-
ные, потенциальные, гироскопические и пр.
8g ( q)
силы, G = G (q ) = -
8qT
- n x m матрица, ор-
тогональная к дополнительным связям (1),
Я ( ?lrw( п\ \
h = h(q,q) =
dq1
Мйл
dq7
q - вектор-столбец
правых частей дважды продифференцированных по времени уравнений связей.
Пусть rank(G) = r < n - ранг матрицы
G в силу избыточности связей (1).
Уравнения (2) представляют собой систему m + n дифференциально-алгебраических уравнений (ДАЕ), линейных относительно обобщенных ускорений q и множителей
Лагранжа ц . При этом матрицы M (q) и
G(q), входящие в систему, являются переменными во времени. Для численного интегрирования подобных ДАЕ на каждом шаге алгоритма необходимо разрешать их относительно ускорений q .
Для этого из (2) необходимо исключить множители Лагранжа. Выражая из первого уравнения системы (2) вектор q и подставляя его во второе уравнение системы, можно получить явное уравнение для вектора р в виде А р = b, (3)
где A = G (M * )_1 GT = G (BT MB )_1
GT -
сим-
метричная положительно полуопределенная п х п матрица ранга г, совпадающего с рангом матрицы О [4], Ь = 0(ВТМВ)* / + Н -
вектор правых частей системы уравнений (3),
причем ранг расширенной матрицы (А \ Ь) так
же равен г .
Векторное уравнение (3) при указанных выше условиях имеет бесконечно много решений. Однако в силу исходной механической постановки задачи, нас интересует только нормальное псевдорешение [4] этой системы уравнений, т.е. решение минимальной длины, доставляющее минимум квадратичной
функции (А м ~ Ь)2 . Это следует из того, что множители Лагранжа м однозначно определяют реакции дополнительных связей (1), и любая механическая система удовлетворяет принципу наименьшего принуждения Гаусса [3], в соответствии с которым уравнения дви-
жения (2) получаются в результате минимизации по вектору q квадратичного функционала
q-ff -/(Gq+h).
Нормальное псевдорешение СЛАУ обычно ищут методом сингулярного разложения матрицы системы [4], что требует больших вычислительных затрат по сравнению с методами Гаусса или Холецкого, которые могут быть применены в случае положительной определенности матрицы системы. Таким образом, при интегрировании системы уравнений (2) необходимо сначала сформировать и решить систему уравнений (3), а затем, найденные множители Лагранжа подставить в первую группу уравнений (2) для определения вектора ускорений q . Видно, что решение уравнений (2) требует множества матричных операций, число которых, в общем случае, возрастает по кубическому закону с ростом размерности задачи.
Заметим, что при малой величине шага интегрирования изменения матриц систем (2) и (3) малы. Это позволяет для снижения вычислительных затрат применить итерационные алгоритмы решения системы уравнений (3), основанные на использовании решений, полученных на предыдущих шагах интегрирования.
В настоящей статье для решения этой задачи предлагается использовать итерационный квазиньютоновский метод, основанный на SR1-формуле поиска минимума квадратичных функций. Особенность метода - возможность находить решение СЛАУ с положительно полуопределенными матрицами системы. Ниже доказываются основные свойства предлагаемого метода.
1. Постановка задачи
Сформулируем следующую задачу: требуется определить вектор Ц = Ц {t) из СЛАУ (3) в момент времени t j = tJ._l + At с
использованием решения в предыдущий момент времени tj_1, где At - малое приращение по времени.
Введем стандартные обозначения СЛАУ. Для этого переобозначим вектор ц = p{t) через «-мерный вектор-столбец неизвестных x . Тогда задачу уточнения решения системы уравнений (3) при малых изменениях коэффициентов уравнений можно переформулировать следующим образом.
Пусть задана СЛАУ Ax = b,
(4)
где А = А + еЩ, А, Щ - известные симметричные и положительно полуопределенные матрицы п-го порядка; Ь = Ь0 + е/0 е Яп - известный вектор-столбец правых частей; еЩ, е/0 - матрица и вектор-столбец возмущений; х е Яп - вектор-столбец неизвестных; е -
малый
f
= O
b
параметр; Щ = О(||^Ц)
Система уравнений (4) совместна. Кроме того, будем считать, что, если матрица возмущений еЩ = цА, то е/0 = /иЬ0. В общем случае положим, что матрица возмущений еЩ непропорциональна матрице А.
Пусть ранг матриц А, А, Щ равен г < п .
Квадрат нормы матрицы:
ПТ (АТА) п
{A\ b)
\\A\\^ =| И 2 = sup
T
n n
В силу симметрии и положительной полуопределенности матрицы А все ее собственные значения являются действительными неотрицательными числами [4, 5]. Будем считать, что Л >... >Л ненулевые собственные значения матрицы А . Оставшиеся " — г собственных значений матрицы А равны 0. Матрица А допускает каноническое (спектральное) разложение на произведение 3-х матриц [4, 5]:
А = и"Х"А"Х"и"Х" = ип/,гКхДтхп , (5)
где А"Х" = (Л,...,Л), Агхг = (Л,-.Д) -диагональные матрицы собственных значений, ипхЛ =(и | ... | ип), ипхГ =(и \ ... | иг) -ортогональные матрицы нормированных собственных векторов матрицы А, соответствующие собственным значениям, ит и = Е и ит и = Е . Заметим, что
пхп пхп пхп гхп пхг гхг '
для положительно полуопределенных матриц 2-нормы и собственные значения связаны соотношением ||А | = Л.
Обозначим через Рг с Я" линейное многообразие векторов с базисом {и,.., иг}, а через Р"—г с Я" подпространство векторов с базисом {иг+1,...,ип} . Рг и Рп—г взаимоорто-
гональны в том смысле, что скалярное произведение любого вектора из Рг на любой вектор из Р^-г равно 0, причем Рг ^ Р^ = Яп.
Из равенства рангов матриц системы А, А и расширенной матрицы системы (А \ Ь) следует, что Ь, Ь0, / е Рг, т.е. вектор Ь можно представить в виде линейной комбинации базисных векторов
Ь=иПг в, (6)
где в - матрица-столбец координат вектора Ь в базисе {и1,...,иг} .
Будем считать, что все собственные векторы матрицы А, соответствующие ненулевым собственным значениям, принадлежат подпространству Рг .
Требуется найти приближение к нормальному псевдорешению х « х0 = А+Ь системы уравнений (4), если известна псевдообратная матрица Н0 = А+ . Приближение к решению рассматриваем в среднеквадратиче-ском смысле, т.е. считаем, что в точке х « х0 выполняется неравенство
(Ах * - Ь)Т (Ах * - Ь)<*б,
где - допустимое абсолютное отклонение. Особенность поставленной задачи заключается в том, что в алгоритме нахождения вектора х* используется матрица Н в качестве начального приближения к псевдообратной матрице А+.
Заметим, что, если известно разложение (5), то псевдообратную матрицу можно найти по формуле [4]:
А = ипхгКхЛхп , (7)
где А-Хг = diag(^1-1,...,Я;1) . Кроме того, нормальное псевдорешение х0 е Рг, т.е. принадлежит многообразию Рг, что следует из следующей формулы, полученной последовательным применением формул (6) и (7):
х0 = АЬ = ипхЖИпипхГР = Р', (8)
где р' = А-^гр - матрица-столбец координат вектора х0 в базисе {и,...,иг} .
Дополнительно, легко показать, что
х0 = А+ Ах0 = АА+ х0. (9)
Например,
А+Ах0 = и А-1 и и А иТи р' = х0.
пхг гхг гхп пхг гхг гхп пхг*
2. Метод решения
Как уже было сказано выше, нормальное псевдорешение системы линейных уравнений (4) есть решение минимальной длины задачи минимизации квадратичной функции
/ (х) = 1 (А х - Ь)2 .
Введем обозначения: V/ (хк ) = гк = Ахк - Ь - градиент / (х) в точке хА, ^ = х^+1 - х^ - направление спуска,
у к = гк+1 - г = ч (хк+1)- v/ (хк) =
= А( хк+1 - хк ) = а:*к.
Пусть матрица А (матрица Гессе функции / ) положительно определена. В любом методе переменной метрики (квазиньютоновском методе) последовательность приближений {хА | к точке минимума задается следующей процедурой:
як = -а1 А-'гк = -акНктк,
хк+1 = хк + 'к, где Нк = А,-1 оценка обратной матрицы А-1 в
точке х . Для сходимости метода длина шага
*
ак в направлении вектора ^ должна удовлетворять условиям Вульфа [6, 7]. Обычно выбирают а* = а^тт/(х^ +аНкгк), т.е. как
решение задачи одномерной минимизации / (х) в направлении ^ . Разные методы переменной метрики различаются формулами пересчета матриц Н или Д.. При этом гарантируется при определенных условиях сходимость последовательности {Нк} к А-1.
В квазиньютоновских алгоритмах, основанных на 8Я1-формуле [8-31], обновление матриц Н выполняется по формуле
„ „ , ('к -НкУк)(Ч -НкУк)Т (10)
Нк+1 = Нк +-:--—г?-. (10)
(*к - Нк У к) У к
Существуют исследования [8-10], которые показывают, что в ряде неквадратичных экстремальных задач 8Я1-метод демонстрирует более высокие результаты, чем другие, более сложно устроенные методы переменной метрики (такие как ББ08, ББР).
Но формула (10) имеет и свои недостатки, которые отмечаются рядом авторов [1114]. В частности, к недостаткам относится возможность обращения в ноль знаменателя
[15] и сложность процедур предотвращения подобного явления.
В настоящее время продолжаются исследования, направленные на улучшение этого метода. Так, в работах [8-10] описываются свойства локальной и глобальной сходимости метода при решении различных задач. Ряд его модификаций предлагается в [17-29].
Применим формулу (10) для решения поставленной задачи решения систем линейных уравнений с симметричными положительно полуопределенными матрицами. Для этого сделаем в SRl-методе следующие преобразования. Зафиксируем параметр а* = 1. Тогда st -Hkук =—Hkrk+l. В соответствии с этим преобразуем формулу (10) пересчета матриц H и переставим местами порядок вычисления H и sk. В результате получим следующий алгоритм.
Алгоритм:
Н0 = А, *i = s0 = H0b, у0 = Axi, ri = у0 — b, к = 0 while rk+^r+i > s2a6c
s = Нк rk+i
if
s у > s
I Jk\ абс
(
Hk+i = Hk -
sl 'I у к
j _ (H к Гк+i ) Гк +i (Hk rk+i У у к ( Hk Гк+1 )( Hk Гк +1 )
Hk Гк +1
(11)
(Hk rk+i) у к
Hk+i = Hk
else
Sk+i = ' end if
k = к +1
rk+1 = Axk+1" Ук = rk+i - rk end while
x0 = xk+\, = Hk-
В дальнейшем покажем, что указанная модификация позволяет сохранить для СЛАУ с симметричными положительно полуопределенными матрицами системы основное свойство методов переменной метрики - сходимость к решению за конечное число шагов. Причем, если известно начальное приближение H0 = A, то число шагов в данном методе не превосходит ранга матрицы возмущений F системы линейных уравнений (4). В результате этого, в применении к поставленной
во введении задачи многократного разрешения уравнений движения механических систем относительно ускорений, возникающей в процессе их численного интегрирования, данный метод сходится быстрее других методов, в которых не вычисляется псевдообратная матрица системы или не используется начальная оценка этой матрицы.
Кроме того, на каждой итерации алгоритма (11) выполняется только две операции умножения матрицы на вектор, требующие порядка п2 арифметических операций. В 8Я1-методе вычисляются три таких произведения.
Также в 8Я1-методе требуются дополнительные
*
операции для вычисления длины шага а* .
Конечно, сокращение числа арифметических операций не прошло бесследно. Ниже показано, что это привело к добавлению еще одной итерации для достижения псевдонормального решения СЛАУ.
Свойства алгоритма (11) в применении к решению СЛАУ с симметричными положительно определенными матрицами системы рассмотрены в работе [30].
Пусть теперь А - симметричная и положительно полуопределенная матрица ранга г < п . В первую очередь заметим, что вектор
является нормальным псевдорешением системы уравнений Лйк = ук, т.е. ^ = Л+ ук .
Покажем, что предлагаемая в (11) схема обеспечивает линейную независимость векторов \ и сходимость последовательности
{хк } к решения СЛАУ (4) за конечное число
итераций.
Докажем пять вспомогательных лемм.
Лемма 1. Если хк е Рг, як е Рг, то
У е Рг .
Доказательство
х^+1 = х + ^ е Рг в силу линейности
многообразия Рг . С использованием формул (5), (6) и (8) построим следующие последовательности равенств (, = к или к + 1):
r =
: Лх, -Ь = (ипхгАгхги1п^в- ипгр =
= (АгХг в'- в ) = ипХг в, ,
где в = Агхгв'-в - матрица-столбец координат вектора г в базисе {«[,...,иг} .
Таким образом, гк, гк+! е Рг, у к = Гк+! - Гк е Р. Лемма доказана.
Дополнительно из леммы 1 следует, что все векторы жк, у^, генерируемые алгоритмом (11), принадлежат линейному многообразию Рг . Это можно обосновать следующей последовательностью рассуждений.
Разложение (7) представим в виде
Н = А = и а- 1 иТ =
0 ^ пхг гхг гхп
г п
= V Л-1и.иТ + 0 V ииТ. ^^^ i i i ^^^ i i
i=l i=г+1
Отсюда, для любого х е Рг следует, что Н0 х е Рг и для любого х е Р±п-г следует, что Н0х = 0 . Из условия задачи ь е Рг. Тогда, проводя вычисления по алгоритму (11), получим, что х1, ж0, У0, ж, Г, Н0г е Рг. Так как матрица Н формируется из матриц Н0 прибавлением к ней матрицы диадного произведения вектора Н0г самого на себя, то матрица Н обладает теми же проективными свойствами, что и матрица Н0. Можно продолжить аналогичные рассуждения на последующих шагах алгоритма (11).
Таким образом, получаем, что все векторы жк, Ук е Р".
Следующие леммы 2-5 в работе [30] доказаны для положительно определенных СЛАУ. Здесь доказательства лемм обобщены на случай положительно полуопределенных СЛАУ.
Лемма 2. Матрицы Нк+Х (к = 0, г -1),
вычисляемые по формулам (11), удовлетворяют квазиньютоновским условиям
Нк+lУi = А У, = Ж, i = (12)
Доказательство
При г = к (к = 0, г -1) выполнение
условия (12) можно проверить непосредственной подстановкой (11) в (12):
( н к Гк+1 )(Н к Гк+1 ) Ук _
(НкГк+1 )Т Ук (13)
= Нк Ук - Нк гк+1 = Нк Ук + Эк - Нк Ук = А+У к.
Дальнейшее доказательство проводим методом математической индукции. Пусть
Нк у, = А У, , , = 0,1-1. Тогда из формулы (11) после умножения Н+1 на У получим
(Нк гк+1 )Г у, = (Ж - Нк У к )г у, =
= У, - УкТНкУ, = (А+ Ук )Т У, - УкТА+Уг = 0.
Hk+i у к - Hk У к --
Следовательно, Нк+1 у1 = Нку1 = А у1, , = 0, к -1. Лемма доказана.
Лемма 3. Последовательность направлений поиска {ж^ } удовлетворяет условиям
'к+1 =ак (Е-НкА)'к, (к = 0,1,...,г -2). (14)
Доказательство
Для начала напомним, что ук = Ажк.
Тогда из формул (11) следует цепочка равенств:
{ / ТТ А
^ (Нк Гк+1 ) Гк+1
, (Нк гк+1 )Т у к = ак (жк - Нк у к ) = ак (Е - НкА) жк,
Hk rk+1-
где ак - —
j _ ( H к Гк+1 ) Гк+1 ( нк гк+1 )Г У к
Лемма доказана.
Лемма 4. Векторы жк и
(к = 0,1,...,г -2) линейно независимы. Доказательство
Сначала покажем, что вектор где
s - - H г
sk+1 н к+1гк+1 ,
Гк+1 - АХ*+1- Ь ■
* * -ш- -ш- * *
sk — а Нк гк — а sk — хк+1 — хк ■
а* = arg min f (xt + aHkrk ), можно предста-
aeR
вить в виде линейной комбинации векторов
sk и sk+1:
s*+1 - —нк+1 (Ах*+1 — b)- —Hk+1 (As* + гк)--— нк+1 (As* + гк+1 — у к )--—aHk+1 Ask — Нк+1Гк+1 + Hk+1 Ук -
- (1 — а) нк+1 Ук + sk+1 - (1 — а ) sk + sk+1.
*
Но вектор sk+1 A -сопряжен с sk, т.к.
ЖкАжк+1 = ЭкАНк+1Гк+1 = (Нк+1 Ажк ) Гк+1 =
= -( Нк+1 Ук )Т гк+1 =- жТ V/ ( хк+1 ) = 0.
Равенство нулю имеет место в силу того, что в точке х*+1 V/ (х*+1 )± як . Сопряженные направления линейно независимы [4].
Следовательно, используя доказанное выше свойство, что вектор '* представим в виде линейной комбинации векторов ' и ' , приходим к заключению, что Е - НА Ф ХЕ и векторы жк и линейно независимы. Лемма доказана.
Хк+1 Хк + Sk '
Лемма 5. Векторы , к = 2,г -1 удовлетворяют условиям Л - сопряженности
яТ.Мк = 0, , = 0, к - 2 .
(15)
Доказательство
Равенства (15) доказываются последовательным применением соотношений (14) и (12) с учетом равенств (9):
к = У'г Ь = ак-1УТг (Е - Нк-1Л) = ак-1 [(Е - ЛНк-1) у, у sk-1 =
^ -1 =
= а
^"ЛЛ+У,) 1= Лемма доказана.
Из лемм 4 и 5 следует, что все направления поиска ^, к = 0, г -1 линейно незави-
симы, следовательно, векторы у к, к = 0, г -1
также линейно независимы.
Из леммы 2 следует, что
Нгуг = Л+у,, V/ = 0, г -1.
(16)
Из рассуждений, выполненных после леммы 1, следует, что, с одной стороны у, Нгу е Рг, а с другой, для любого собственного вектора и е Р,п-г Н и = Л+и = 0.
г г ± г г г
Построим матрицу ^ =(у0..у-1 иг+l,...,ип) . Тогда равенства (1 6) можно записать в виде
Н X = Л+X .
г пх п пх п
Матрица Хихи по построению имеет полный ранг. Отсюда следует, что Нг = Л+.
Таким образом, если все вычисления проводились точно, последний шаг алгоритма
х , = х - Н г = х - Л V
г+1 г г г г г
приведет нас в точку х 0 , являющуюся нормальным псевдорешением системы уравнений (4). Заметим, что в точной арифметике для достижения точки х 0 в поставленной задаче решения СЛАУ с симметричной положительно полуопределенной матрицей системы достаточно, чтобы Н0 = Л+ .
Мы доказали, что рассматриваемый метод является конечно сходящимся. Однако для того чтобы использовать его в качестве итерационного метода решения СЛАУ, необходимо обеспечить монотонность убывания последовательности погрешностей решения
{||х - х*|Ц. В этом случае решение в средне-
квадратическом смысле может быть достигнуто за меньше, чем г +1 итерацию к (к < г +1).
Найдем, при каких условиях метод обеспечивает монотонную сходимость последовательностей {Нк} и {хк} к Л+ и х0 соответственно.
Введем обозначения:
Л = Н+, Нк = Л+ , Л-Лk, еАНк = Нк -Н , где £ - малый параметр, е¥к - симметричная
матрица возмущений (е¥0 = е¥) .
Обозначим через к(Л ) = |||-||Л+|| =
А
Лк|| -||Нк\\ = — - эффективное число обу-
А
словленности матрицы Л .
Сделаем предварительные выкладки. Используя свойство псевдообратных матриц, зафиксированное в равенстве (9), и матричный ряд Тейлора [4] для матрицы полного
ранга (Е + еЛ+¥) , можно получить следующую последовательность равенств:
Н = Л+ = (Лк + е¥к )+ = = (Лк (Е + еЛ+ ¥к ))+ = (Е + еЛ+ ¥к )" Л+ = = (Е-еНк¥к + е2 (Нк¥к)2 -...)Нк = = Нк-(Е-еНк¥к +...)еНк¥кНк = = Нк-(Е + еНк¥к )-1 еНк¥кНк.
Таким образом, возмущения прямой и обратной матрицы связаны следующим матричным соотношением:
еАНк = Н - Н = (Е + еНк¥к )-1 еНк¥кНк. (17)
Из равенства (17) следует, что, если матрицы Н , ¥ положительно полуопределены и имеют ранг г < п , а Е + еНк¥к - положительно определенная матрица, то ДНк также положительно полуопределенная ранга г < п матрица как произведение симметричных определенной и полуопределенной матриц.
Это свойство доказывает следующая лемма.
Лемма 6. Пусть С = Сихи - симметричная определенная матрица ранга п и О = Оихи - симметричная полуопределенная матрица ранга г < п . Тогда СD - симметричная полуопределенная матрица ранга г < п .
Доказательство
В соответствии с (5) матрицы С и Б допускают канонические разложения на произведение 3-х матриц:
С = и Ас иТ Б = V А° V'
С и пхпАпхпи пхп , Б г.....А— -г -
пхг"" "гхг гхп -
где Ас = А°пхп, АБ = А^хг - диагональные матрицы положительных собственных значений, и = иихИ, V = - ортогональные матрицы нормированных собственных векторов и , V . Для простоты изложения будем считать, что первые г векторов и., V принадлежат
линейному многообразию Рг . Пусть Ж = и^ - ортогональная матрица координат щ векторов V в базисе {и,} . Тогда матрица
СD = UАCUTVАCVТ = = иАс (ЖАБЖТ )иТ
положительно полуопределена, если Ас (ЖАБЖТ) - положительно полуопределенная матрица. Разложим Ас (ЖАБЖТ) на сумму матриц ранга 1:
Ас (ЖА°ЖТ ) = ]Г (аХ ) .
1=1
Матрицы АсХБ положительно определены. Легко видеть, что
Ас (ЖА°ЖТ )щ = Щ (АсХБ )щ > 0
для всех к = 1, г и
г
иТАс (ЖАБЖТ )и = иТ]Г (АХ )щ (Щик ) = 0
для всех к = г +1, п . Отсюда следует положительная полуопределенность матрицы СD в пространстве Яп и положительная определенность на линейном многообразии Рг . Лемма доказана.
Найдем явную зависимость матрицы Н+1 от матриц возмущений еЕ^ , еДН .
Для этого заметим, что
Н г = Н (г + у ) = -ж + Н Аж =(Н А - е) ж =
к к+1 к \ к к} к к к \ к ! к
= ((Н + еДНк ) А - Е) -к
= едн а® ,
к к к
Н г = (Н А - е) ж =( Н (а + еЕ )-е) ж =
к к+1 \ к /к\к\к к) ) к
к к к к к к к к к к к к
Здесь использовано свойство псевдообратных матриц НкАкНк = Нк .
Дополнительно заметим, что г е Рг, и
используя свойство (9), можно получить обратное равенство:
г = еЕ ж . (18)
к+1 к к у '
Кроме того,
(н г )Т у = жТ (еАДН А) ж . \ к к+1! ■'к к \ к ! к
Подставим полученные соотношения в
формулу (11) для Н . Окончательно придем
к следующему выражению
(ДН Аж )(ДН Аж )Т
к к)\ к к)
Н = Н -е-
к+1 к
жТ (АДН А) $
кк
(19)
к \ к ) к Кроме того, из теории известно [31, 32], что любой метод переменной метрики может быть сформулирован как для обратной матрицы Н = А 1, так и для прямой исходной матрицы системы А = А +еЕ0 . Алгоритм пересчета матрицы приближений А к прямой матрице А можно получить из формул (10) заменой: Нк ^ А, жк ^ У к, У к ^ ^ (У к - Ак жк)(У к - Ак жк)Т
3 к ■
\Т
Ак+1 = Ак + ■
(У к - Ак жк)Т жк
_л , гк+1 гк+1 _ а , с.
= Ак + Т = Ак + е Т"
(20)
к+1 жк
жк Ек жк
Здесь дополнительно использованы равенства (18). Легко проверить непосредственной подстановкой, что, если Н0 = Д+, то последовательности приближений удовлетворяют условию взаимной обратимости матриц Н и
Нк+1 Ак+1Нк+1 Нк+1
Ак +1 :
Ак+1Нк+1Ак+1 = Ак+1 и квазиньютоновским условиям (см. лемму 1 и формулы (12) и (9)):
Нк+1У,= ж,=А + Уг.
А, ж = у= Аж , , = 0, к.
к+1 ,' '
(21)
Теперь сформулируем и докажем достаточное условие монотонности убывания последовательности погрешностей решения
Теорема. Если выполнены все условия, зафиксированные в постановке задачи, и
1 К || = ^, (22)
11еЕ0|| <-
к
(А))'
Н
то алгоритм (11) обеспечивает выполнение следующих свойств:
1=1
1. {Нк} - последовательность положительно определенных матриц по отношению к любым векторам ц е Рг, т.е. цТНкц > 0
для всех ц е Рг, к = 1, г .
2. Последовательности {|£¥| 1},
{||Н - Нк|} и {||хк - х0|} монотонно убывают.
Доказательство Из условия (22) следует, что 1
К
( Л )
= шт Л = Л .
л( л ) г
То есть спектральный радиус матрицы е¥0 меньше, чем минимальное собственное значение матрицы Л . Используя это неравенство и каноническое (спектральное) разложение матриц легко показать (аналогичными рассуждениями, приведенными в лемме 6), что, если матрицы ¥ и Л имеют ранг г и положительно
определены на линейном многообразии Р г , то минимальное ненулевое собственное значение матрицы Л + е¥0 положительно и
цТЛц = цТ (Л +е¥0) п > 0
для всех ц Ф 0 е Рг и цТЛц = 0 для всех
ц е Я" \ Рг, т.е. условие (22) является достаточным условием положительной определенности на Рг матрицы Л ранга г при любом знаке малого параметра £ . Кроме того, в силу неравенства (22), положительной полуопределенности матрицы Н0¥0 (по лемме 6) и свойств 2-норм матриц, можно выписать следующую цепочку неравенств:
1 >11НЛ • ИII ^ ||еН0¥01| =|£|. шах л = |4лп„,
л(Н0)
где Лтж > 0 - максимальное собственное значение матрицы Н0¥0.
Следовательно, повторяя рассуждения, приведенные при доказательстве положительной определенности матрицы Л + е¥0 на
линейном многообразии Рг , заключаем, что цТ (Е-£|Н0¥0)ц > 0 для всех ц Ф 0 е Рг. Это означает, что матрица Е + еН
положительно определена на Р г при любом знаке малого параметра £ .
Отсюда следует, используя равенство (17), что матрица ДН0 также является положительно определенной на Рг .
Дальнейшее доказательство можно провести методом математической индукции. Сначала проверим справедливость теоремы при к = 1.
Заметим, что по доказанным леммам все направления поиска , к = 0, г -1 принадлежат линейному многообразию Рг и линейно независимы. Следовательно, любое подмножество этих векторов можно рассматривать в качестве базиса некоторого линейного многообразия. Обозначим через Рг-к - линейное многообразие с базисом {«г, sг _1,..., sk}. Причем Р1 с Р2 с... с Рг.
Докажем положительную определенность матрицы Н по отношению к любому
вектору ц е Рг, используя равенство (19).
Нам необходимо показать, что для всех ц Ф 0 е Рг справедливо неравенство
цтн1 ц = цТн0 ц
- £
цТ (ДН As )(ДН As
' \ 0 к/\ 0 к
( АДН 0 А) *
)Т ц (23) —>0.
По доказанным выше условиям знаменатель в формуле (23) положителен. Тогда, если е < 0 , то
цТн, ц = цТн0 ц + \е\
> 0.
^ (ЛДН Л)s
0 V 0 0
где и = s ЛДН ц .
0 0 '
Следовательно, в этом случае матрица Н положительно определена на Рг .
Пусть теперь е > 0 . Приведем неравенство (23) к общему знаменателю и исследуем знак выражения, получаемого в числителе
/ = ( цТЯц)( ^ )
-еи2 =
:( цТНц)(
s1 ЛДН Лs ) +
0 0
)■
+е (цТ ДН ц) ( sToЛАHoAso ) -ецТ (ДНоЛ^ )(Д^0Л^ )Т ц.
0
0
2
и
В силу симметрии и того, что цТ ДН^п > 0 для всех п £ РГ и цт ДН^п = 0 для
всех ц £ Я" \ РГ, матрицу ДНо можно разложить (с помощью спектрального разложения) на произведения двух матриц [4, 5]:
ДН = V Л ¥т =(>Л1/2 )(Л12ГТ ) =
0 "хг Гхг Гх" \ /\ /
= С ст = ССт.
"хг Гх"
Тогда функцию / можно представить в виде: / = (цТНц) (атА% )2 + е (Стц)2 (СтА% )2 -
~8\(сТп)Т (^0)_
Выражение
(Ст п)2 (стА*0 )2-Г(Ст ц )Т (стА*0 )
> 0
в силу справедливости для всех а, Ь £ Я" не-
2 2 / т \2
равенства Коши-Буняковского а Ь > (а Ь) .
Таким образом, окончательно получаем следующую цепочку неравенств:
/ >(пНп)(СтА*о )2 >0, для всех п Ф 0 £ Рг и ж Ф 0, доказывающих требуемые неравенства птН п > 0 для всех п £ Рг.
Теперь покажем, что птЕ п > 0 д
всех п £ Рг и птЕ п > 0 для всех п £ Рг Используя равенства (20), выведем явную за-
висимость матрицы возмущений еЕ от
0 '
еЕ = А - А1 =( А - А0)-(А1 - А0 ) =
= е
Е -
Е0 Ж0 (Е0 Ж0 ) ^
Жт Е0 Ж0 у
(24)
Во-первых, еЕж0 = е(ЕЖ -Ежо) = 0 , т.е. для всех п £ Рг \ Рг 1 выполняется птЕ п = 0. Во-вторых, матрица Е является симметричной
положительно определенной матрицей полного ранга на линейном многообразии Рг и ее можно разложить на произведение двух матриц (аналогично матрице Н ):
^ =С Ст =ССТ.
ептЕп = еГптЕп - п'Е0Т(/0*0)т п^
у Ж0 Ж0 у
Е0 Ж0
= е
Так как векторы п и ж0 не параллельны, а матрица С имеет полный ранг на линейном многообразии Рг , то скалярное произведение
(СТП) (Сгж0) строго меньше произведения длин \Стч\ и |е~'7 .ч0|. Поэтому выражение
(С тп)2(С т*0)2 -((С тп)Т (С т*0))2 >0.
Т
Следовательно, п Е1п > 0 при любом
тлТ -1
п £ Р , т.е. матрица положительно определена на линейном многообразии Рг 1 и положительно полуопределена на Рг . Кроме того, из этого неравенства следует и неравенство
для норм еЕ < еЕ .
Дополнительно заметим, что из доказанных свойств матриц Е , Н , соотношений
(17) и леммы 6 следует, что матрица возмущений ДНХ должна обладать теми же свойствами, что и матрица возмущений Е .
Теперь докажем, что
еН Е < еН Е
1 1 0 0
1 и Н-Н < Н-Н
1 0
Из (19) следует, что еДЯ = Я -Н = (я -Но)-(н-Но) =
\Т \
= е
ДН -
0
(ДН Аж )(ДН Аж )
0 0 0 0
( АДН 0 А) ж
(25)
Легко заметить совпадение структур равенств (24) и (25) для матриц возмущений еЕ и е ДН^. Поэтому, если повторить для
матрицы еДЯ^ все рассуждения, выполненные при доказательстве свойств матрицы , то можно окончательно доказать, что
т Г 1
пт дн п > 0 при любом п £ Р , т.е. матрица возмущений ДНХ также положительно опре-Тогда при любом п £ Рг 1 равенство (24) делена на линейном многообразии Рг 1 и можно преобразовать к виду: птДН1 п = 0 для всех п £ Рг \ Рг-1.
"х г гх."
2
Кроме того, из равенства (25), по аналогии с равенством (24), следует и неравенство
сти
для
норм
£Ш
<
еДЯ
или
И - И < И - И
1 о
Вернемся к равенству (17). Его можно преобразовать к виду еАНк = ЕИЕкИк = = ЕИкЕкИ. Подставляя последнее в полученное выше неравенство для норм получаем еще одно требуемое неравенство:
еИЕ < еИ Е < 1
1 1 о о
Теперь
II << IIХд X
Заметим, что
осталось о II
показать,
что
x х° ~ x н г х°
= х0 - хо - И0 (Лх0 - Ь) = = х0 - х0 -И0 (Лх0- Ах0) = = (Е - И0А)(х0 - х0). Следовательно, ||х - х0|| < ||х0 - х0 , если ||Е - И0А| < 1. Последнее неравенство выполняется, так как в силу (9)
ПТ (Е - И о (Ао + е^ ))2 г,
||Е - И о А|| = 8ир-
П^о
Т
П П
= 8Ир
П^о
ПТ (Е - И о Ао )2 п - ПТ (£И о Е )2 Ч
т
п п
= \\еиоЕо||2 < 1. Таким образом, утверждения теоремы доказаны при к = 1.
Кроме того, мы попутно доказали следующий ряд утверждений, которые необходимы для дальнейшего доказательства теоремы по индукции. При к = 1:
1. ПТ¥1П > о при любом п е РГ~к и ПТ¥1П = о
при любом п е РГ \ РГ к, т.е. матрица положительно определена на линейном многообразии
Р
г-к
и положительно полуопределена на Рг. < 1.
2. еИЕ < еИ Е
1 1 о о
Если теперь принять, что эти два условия справедливы при ] = к , то очевидно, что мы можем повторить все рассуждения теоремы при ] = к +1, 1 < к < п -1.
Таким образом, условие (22) является достаточным условием монотонной сходимо-
последовательностей норм {Це-Е 11} , {IИ -Ик\} и {||х - х0|}. Теорема доказана.
Замечание. Дополнительно из равенства (25) следует то, что в процессе выполнения итераций по формулам (11) знаменатель
(Ик гк+1 )Т У к = 4 (ААИкА) sk никогда (кр°ме
последнего шага) не обращается в ноль при соблюдении всех условий, указанных в доказанной теореме и постановке задачи.
Приведенные выше рассуждения позволяют установить еще один факт. Оказывается, что если матрица возмущений еЕ имеет ранг меньше г, то описываемый итерационный метод сойдется меньше, чем за г итераций.
Пусть матрица возмущений
еЕ = А - А имеет ранг г1 (г1 < г < п). Обозначим через Рг1 с Рг с Яп линейное многообразие векторов с базисом {иг,...,и, }, состоящим из собственных векторов, соответствующих ненулевым собственным значениям матрицы еЕ , а через РП-г1 с Яп подпространство векторов с базисом {иг+1,..., ип},
включающим собственные векторы, соответствующие нулевым собственным значениям
матрицы еЕ . Тогда (А - А ) И = еЕщ = о и
(А + - И0) и = о для всех и е Р^ г' .
В этом случае вычисления, выполненные по формулам алгоритма (11), в соответствии с доказанными леммами приведут к тому, что все векторы У принадлежат многообразию Рг1, к = о, г -1 и линейно независимы. Кроме того, согласно рассмотренной теореме, все матрицы И , к = 1, г1 будут удовлетворять следующим условиям: пТИ п > о для всех п е , И п = Ик ]Ц для всех
п-г
п е Р 1 и квазиньютоновским условиям
Ик У г = А+У г, г = о, к -1 .
Отсюда следует, что столбцы матрицы И г можно найти из системы линейных уравнений
(А+ - ИГ1) Хпп = о, (26)
с матрицей ХпХп =( У,...У, -! ил +1,..., ип ) .
Матрица системы (26) Хихи по построению имеет полный ранг. Следовательно, система (26) имеет единственное решение
ИГ1 = А+.
Таким образом, можно заключить, что итерационный алгоритм (11) при выполнении всех условий, указанных в постановке задачи, и доказанной теореме, обеспечивает сходимость последовательности приближений {хк }
к точному решению х о системы уравнений (4) за конечное число шагов, не превышающих ранг матрицы возмущений еЕ = А - А.
Заметим, что при численном решении уравнений движения механических систем вида (2) обычно на каждом шаге интегрирования эффективный ранг матрицы возмущений много меньше ранга матрицы системы. Поэтому следует ожидать, что приближенное решение системы уравнений (2) рассматриваемым итерационным методом в практических задачах будет достигаться за небольшое число итераций.
4. Пример моделирования механической системы с дополнительными связями
Применение рассмотренного итерационного метода проиллюстрируем на примере моделирования динамики одной механической системы с двумя замкнутыми кинематическими цепями, представляющей собой модель весельной лодки с механическим приводом (рис. 1).
Рис. 1. Модель весельной лодки
Модель включает в себя два весла, закрепленные в уключинах, обеспечивающих угловые перемещения весел в двух плоскостях, и двигатель, задающий вращательное движение, с двумя кривошипами, соединенными сферическими шарнирами с веслами. Система имеет одну степень свободы.
Уравнения движения построим в форме уравнений (2). Для этого расцепим соединения весел с кривошипами. В результате получим расчетную схему системы с пятью степенями свободы, состоящую из трех тел (двух весел и двигателя с кривошипами), представленную на рис. 2. Обобщенные координаты указаны на рис. 1 (т = 5).
Порядок обобщенных координат:
ч = (P, г^^ Уз, аз).
Рис. 2. Расчетная схема лодки
Введем обозначения: г0 - расстояние между кривошипом и уключиной вдоль Оу\ и высота кривошипа, гй = ОС - расстояние вдоль весла от уключины до крепления с кривошипом, ОС = (о, го, го), ОС = (о, - го, го),
^ - осевой момент
= 2г0,
OO =-2Го,
инерции механизма вращения с кривошипами, У2 = У3 - осевой момент инерции весла, т2 = т - масса весла, г2с = гЪс - расстояние от уключины до центра тяжести весла.
Вращение кривошипов обеспечивается внешним упруго-демпфирующим моментом:
Л Л
Мр= cp
+dr.
T„
2л
t - k„
T
2л
-sin-
2л T,
t
-Р
T
V ТР
л i 2л
1 - k„ cos—t
T
Р
Р
\ \ -р
У
+
где Tp - период оборота кривошипа, 2л/Tp -максимальная круговая частота вращения, Cp, dp - коэффициенты упругости и демпфирования. Сомножители 1 — kpCOs2лt|Tp и t — kp Tpj2л • sin 2лt|Tp определяют переменную скорость вращения и введены для имита-
ции торможения скорости вращения за счет сил сопротивления при погружении весел в воду, кр< 1 - коэффициент торможения.
Переход к исходной одностепенной модели будем осуществлять наложением шести дополнительных связей, обеспечивающих совпадение координат точек С = С ' и С = С '
соединения весел с кривошипами (n = 6). Очевидно, что две из них избыточны.
Уравнения связей имеют вид: g = -r0 • sin f + r • cos« • siny = 0,
g2 = Г \ 1 - \ - Г • C0S «2 • C0S ^2 = 0 ]¡ Г
g3 = r0 • cosf -rh • sin« = 0, g = -r0 • sin f — r • cos« • siny = 0,
gs = -ru \ 1 - \ + Г •Cos «з •Cos Гз = 0,
V rh
ge = -r0 • cosf + rh • sin« =
Компоненты системы дифференциально-алгебраических уравнений вида (2) в нашем случае имеют вид:
M* = J, M* = J2 cos2 «, M* = J2, M* = J3 cos2 «, M5*5 = J3, Mj = 0, i Ф j ,
fl=M{J> f2=J2- SÍn2«2 • Г2«2,
/3 = • • C0S «2 " ■ SÍn 2«2 • /I,
/4 = J4 -sin2« -y«,
f5 = m3g ■ r3c ■ cos«з -^ J3 • sin2a3 ■ f3,
Gn = -r0 • cosf, Gu = ^ • cos« • cosy,
G13 = -rh • sin « • sin y2, G14 = G15 = 0,
G22 = rh • cos« • siny, G23 = rh • sin« • cosy,
G21 = G24 = G25 = 0
Gl = ^ • sin f, G33 = • cos «, С —С —С — п
W32 ^^34 ^^35 0,
G41 = -0 • cos A, G44 = -rh • cos« • cos y3,
G45 = rh • sin« • sin y3, G42 = G43 = 0,
G54 = -rh • cos« • sin y3, G55 = -rh • sin« • cos y3,
G51 = G52 = G53 = 0
G61 = r0 • sin f3, Ges = rh • cos «2,
Ge2 = Ge3 = Ge4 = 0,
I\=-r0- sin p- p2 +rh- cos « ■ sin y2 ■ (y2 + á\ ) +
+2rh ■ sin« • cosy • «у,
h2 = -rh ■ cos a2 ■ cos у • (y + á\ ) +
+2rh ■ sin« • siny • «у,
h3 =—r0 -cos P~ p —rh -sin« •«,
/г4 = -r0 - sin p- p2 -rh- cos « • sin у • (y2 + «2 j
-2гЛ • sin« • cosy • «у,
h5=rh ■ cos « • cos у • (у2 + «2) -
-2гЛ • sin« • siny • «у,
h6 = -r0 ■ cos p-p2 +rh- sin« • «2.
В уравнения движения в форме Ла-гранжа 1 рода вида (2) для их замыкания используются дважды продифференцированные уравнения дополнительных связей g{q) = 0.
Для предотвращения ухода со связей добавим в эти уравнения стабилизирующие слагаемые по методу Баумгарта [2]:
g + 2 Р 2 + к 2 = 0,
О Г 5
где Д, к^ - коэффициенты демпфирования и
жесткости.
Приведем значения всех конструктивных параметров:
r = 0.5, r = i, J = 1, m = m = 4,
Г0 „ = Г3 с = 1.25, J2 = J3 = 13, Д = 5 ,
к = 10000 Г„ = 5 « = 10000 = 100
* > д > д > д >
кД= 0.7, g = 9.81. Начальные условия:
q(0) = ^0,0, "|,0,0, ,
g(0) = (0, о, о, о, о, 0).
Отрезок интегрирования t е [0,12].
Уравнения движения в форме (2) интегрировались в САВ Mathematica. Система линейных уравнений (3) для определения множителей Лагранжа решалась итерационным методом (11). Начальное приближение псевдообратной матрицы = А вычислялось встроенной процедурой PseudoInverse.
Результаты численного моделирования представлены на рис. 3 а, б, в.
а)
б)
в)
Рис. 3. Результаты моделирования
На графиках угловых перемещений весел хорошо заметно сымитированное замедление их движений при погружении в воду.
На рис. 4, для примера, приведены графики отклонений дополнительных связей между одним веслом и кривошипом.
Графики на рис. 4 показывают, что представленный алгоритм решения задачи обеспечивает хороший уровень точности в удовлетворении дополнительных геометрических связей на длительном промежутке интегрирования.
Рис. 4. Точности удовлетворения дополнительных связей
В заключение заметим, что в рассмотренной механической задаче на промежутке моделирования среднее число шагов в итерационном алгоритме (11) составило к = 1.994, что в два раза меньше ранга матрицы системы уравнений дополнительных связей и в три раза меньше их числа.
Таким образом, результаты моделирования подтверждают описанные в статье свойства итерационного алгоритма решения систем линейных уравнений с положительными полуопределенными матрицами и возможность его применения при моделировании динамики механических систем с дополнительными связями.
Заключение
В статье предложена модификация обратного итерационного метода Пауэлла-Бройдена, основанного на симметричной формуле ранга один (8Я1), предназначенная для решения возмущенных СЛАУ с симметричной положительно полуопределенной матрицей системы. Метод используется для решения относительно ускорений уравнений движения механических систем с дополнительными связями при их численном интегрировании. Рассмотрены основные свойства алгоритма.
1. Показано, что все направления поиска решений - векторы я , кроме двух последних, являются А -сопряженными. Два последних направления на каждой итерации линейно независимы. Тем самым доказано, что метод глобально сходится за конечное число итераций и точное решение может быть найдено не более чем за г + 1 шаг, где г - ранг матрицы системы.
2. Найдены условия, при которых метод локально сходится, т.е. последовательность приближений монотонно стремится к точному решению.
3. Показано, что если матрица возмущений имеет ранг r\ < r, и известно начальное приближение псевдообратной матрицы системы, то алгоритм (\\) сходится за r\ + \ шаг. Таким образом, приближенное решение возмущенной системы n линейных уравнений, может быть найдено с помощью алгоритма (\\) много меньше, чем за n + 1 шаг.
4. На одном примере проиллюстрировано применение метода при интегрировании уравнений движения механических систем с избыточным числом дополнительных связей. Показано, что количество итераций не превосходит числа линейно независимых связей и остается небольшим на всем промежутке интегрирования.
Список источников
\. Величенко В.В. Матрично-геометрические методы в механике с приложениями к задачам робототехники. М.: Наука, 1988. 280 с.
2. Лилов Л.К. Моделирование систем связанных тел. М.: Наука, 1993. 272 с.
3. Суслов Г.К. Теоретическая механика. М.: Гостехиздат, 1946. 656с.
4. Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, \999. 548с.
5. Шевцов Г.С. Численные методы линейной алгебры: учеб. пособие / Г.С. Шевцов, О.Г. Крюкова, Б.И. Мызникова. М.: Финансы и статистика: ИНФРА-М, 2008. 480 с.
6. Wolfe Р. Convergence conditions for ascent methods // SIAM review, Apr. \969, Vol. \\, № 2. P.226-235.
7. Wolfe Р. Convergence conditions for ascent methods. II: Some corrections // SIAM Review, Apr. \97\, Vol. \3, № 2. P. \85-\88.
8. Conn A.R., Gould N.I.M., Toint Ph.L. Convergence of quasi-Newton matrices generated by the symmetric rank-one update // Mathematical Programming. \99\. 50, № \. Р. \77-\95.
9. Byrd R.H., Khalfan H.F., Schnabel R.B. Analysis of a symmetric rank-one trust region method // SIAM J. on Optimization. \996. 6. Р.\025-\039.
\0. Khalfan H.F., Byrd R.H., Schnabel R.B. A theoretical and experimental study of the symmetric rank-one update // SIAM J. on Optimization. \993. 3. Р. \-24. \\. Аоки М. Введение в методы оптимизации.
М.: Наука, \977. 344 с.
\2. Дэннис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М.: Мир, 1988. 440 с.
\3. Поляк Б.Т. Введение в оптимизацию. М.: Наука, 1983. 384 с.
\4. Химмельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. 536 с.
\5. Nocedal J., Wright S.J. Numerical Optimization. Berlin: Springer, 2006. 664 р.
\6. Фиакко А., Мак-Кормик Г. Нелинейное программирование. Методы последовательной безусловной минимизации. М.: Мир, 1972. 240 с.
\7. Farzin M., Malik Abu H., Wah J.L. Multisteps symmetric rank-one update for unconstrained optimization // World Applied Sci. J. 2009. 7, № 5. Р. 6\0-6\5.
\8. Farzin M., Malik Abu H., Wah J.L. Memor-yless modified symmetric rank-one method for large-scale unconstrained optimization // American J. of Applied Sciences. 2009. 6, № \2. Р.2054-2059.
\9. Farzin M., Malik Abu H., Wah J.L. Convergence of symmetric rank-one method based on modified quasi-Newton equation // J. of Math. Res. 20\0. 2, № 3. P. 97-\02.
20. Yueting Y., Chengxian X. A switching algorithm based on modified quasi-Newton equation // Numer. Math. A J. of Chinese Universities. 2006.\5, N 3. 257-267.
2\. Wah J.L., Malik Abu H. Convergence of a positive definite symmetric rank-one method with restart // Advanced Modeling and Optimization. 2009. \\, № 4. P. 423-433.
22. Wah J.L., Malik Abu H. A restarting approach for the symmetric rank one update for unconstrained optimization // Comput. Optim. Appl. 2009.42. P.327-334.
23. Al-Baali M., Fuduli A., Musmanno R. On the performance of switching BFGS/SR\ algorithms for unconstrained optimization // Optimization Methods and Software Vol. \9, № 2, April 2004, Р. \53-\64.
24. Oztoprak F.& §., ilker Birbil S.I. A Symmetric Rank-One Quasi-Newton Line-Search Method Using Negative Curvature Directions // Optimization Methods and Software. Vol. 26, 20\\. P. 455-486.
25. Farzin M., Malik Abu H., Wah J.L. A symmetric rank-one method based on extra updating techniques for unconstrained optimization // Computers and Mathematics with Applications. 62 (20\\). Р. 392-400.
26. Farzin M., Wah J.L. On the performance of a new symmetric rank-one method with restart for solving unconstrained optimization problems // Computers and Mathematics with Applications. 64 (20\2). Р. 2\4\-2\52.
27. Shu-Zhen Lai, Hou-Biao Li, Zu-Tao Zhang. A Symmetric Rank-One Quasi-Newton Method for Nonnegative Matrix Factorization // ISRN Applied Mathematics. Vol. 20\4, Article ID 846483. \\ p.
28. Benson H.Y., Shanno D.F. Cubic regulariza-tion in symmetric rank-\ quasi-Newton methods // Math. Prog. Comp. 20\8. \0. P. 457-486.
29. Indrapriyadarsini S., Mahboubi S., Ninomiya H., Kamio T., Asai H. Accelerating Symmetric Rank-\ Quasi-Newton Method with Nesterov's Gradient for Training Neural Networks // Algorithms 2022, \5, 6. \-\6p.
30. Иванов В.Н. Основные свойства обратного итерационного алгоритма решения систем линейных уравнений с положительно определенными матрицами // Вычислительные методы и программирование. 2012. Т. 13, № 2. С. 366-376.
3\. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. М.: Мир, 1985. 509 с.
32. МинуМ. Математическое программирование. Теория и алгоритмы. М.: Наука, 1990. 488 с.
References
\. Velichenko V.V. Matrichno-geometricheskie metody v mekhanike s prilozheniyami k zadacham robototekhniki. M.: Nauka, \988. 280 s. (In Russ.).
2. Lilov L.K. Modelirovanie sistem svyazannyh tel. M.: Nauka; \993. 272 s. (In Russ.).
3. Suslov G.K. Teoreticheskaya mekhanika. M.: Gostekhizdat; \946. 656 s. (In Russ.).
4. Golub Dzh., Van Loun Ch. Matrichnye vychisleniya. M.: Mir; \999. 548 s. (In Russ.).
5. Shevcov G.S. Chislennye metody linejnoj al-gebry: ucheb. posobie / G.S. SHevcov, O.G. Kryukova, B.I. Myznikova. M.: Finansy i statistika: INFRA-M; 2008. 480 s. (In Russ.).
6. Wolfe Р. Convergence conditions for ascent methods. SIAM review, Apr. \969;(\\(2)):226-235.
7. Wolfe Р. Convergence conditions for ascent methods. II: Some corrections. SIAM Review, Apr. \97\;(\3(2): \85-\88.
8. Conn A.R., Gould N.I.M., Toint Ph.L. Convergence of quasi-Newton matrices generated by the symmetric rank-one update. Mathematical Programming. \99\;(50(\):\77-\95.
9. Byrd R.H., Khalfan H.F., Schnabel R.B. Analysis of a symmetric rank-one trust region method. SIAM J. on Optimization. \996;(6): \025-\039.
\0. Khalfan H.F., Byrd R.H., Schnabel R.B. A theoretical and experimental study of the symmetric rank-one update.SIAM J. on Optimization. \993;(3): \-24.
\\. Aoki M. Vvedenie v metody optimizacii. M.: Nauka; \977. 344 s. (In Russ.).
\2. Dennis Dzh., Shnabel' R. Chislennye metody bezuslovnoj optimizacii i resheniya nelinejnyh uravnenij. M.: Mir; \988. 440 s. (In Russ.).
\3. Polyak B.T. Vvedenie v optimizaciyu. M.: Nauka; \983. 384 s. (In Russ.).
\4. Himmel'blau D. Prikladnoe nelinejnoe pro-grammirovanie. M.: Mir; \975. 536 s. (In Russ).
\5. Nocedal J., Wright S.J. Numerical Optimization. Berlin: Springer; 2006. 664 р.
\6. Fiakko A., Mak-Kormik G. Nelinejnoe pro-grammirovanie. Metody posledovatel'noj bezuslovnoj minimizacii. M.: Mir; \972. 240 s. (In Russ.).
\7. Farzin M., Malik Abu H., Wah J.L. Multisteps symmetric rank-one update for unconstrained optimization. World Applied Sci. J. 2009;(7(5):6\0-6\5.
\8. Farzin M., Malik Abu H., Wah J.L. Memor-yless modified symmetric rank-one method for large-scale unconstrained optimization. American J. of Applied Sciences. 2009;(6(\2):2054-2059.
\9. Farzin M., Malik Abu H., Wah J.L. Convergence of symmetric rank-one method based on modified quasi-Newton equation. J. of Math. Res. 20\0;(2(3):97-\02.
20. Yueting Y., Chengxian X. A switching algorithm based on modified quasi-Newton equation. Numer. Math. A J. of Chinese Universities. 2006;(\5(3):257-267.
2\. Wah J.L., Malik Abu H. Convergence of a positive definite symmetric rank-one method with restart. Advanced Modeling and Optimization. 2009;(\\(4):423-433.
22. Wah J.L., Malik Abu H. A restarting approach for the symmetric rank one update for unconstrained optimization. Comput. Optim. Appl. 2009;(42):327-334.
23. Al-Baali M., Fuduli A., Musmanno R. On the performance of switching BFGS/SR\ algorithms for unconstrained optimization. Optimization Methods and Software. 2004; (\9(2): \53-\64.
24. Oztoprak F.& §., ilker Birbil S.I. A Symmetric Rank-One Quasi-Newton Line-Search Method Using Negative Curvature Directions. Optimization Methods and Software. 2011;(26):455—486.
25. Farzin M., Malik Abu H., Wah J.L. A symmetric rank-one method based on extra updating techniques for unconstrained optimization. Computers and Mathematics with Applications. 2011;(62):392-400.
26. Farzin M., Wah J. L. On the performance of a new symmetric rank-one method with restart for solving unconstrained optimization problems. Computers and Mathematics with Applications. 2012;(64):2141-2152.
27. Shu-Zhen Lai, Hou-Biao Li, Zu-Tao Zhang. A Symmetric Rank-One Quasi-Newton Method for Nonnegative Matrix Factorization. ISRN Applied Mathematics. 2014;(2014): 11. Article ID 846483.
28. Benson H. Y., Shanno D. F. Cubic regulariza-tion in symmetric rank-1 quasi-Newton methods. Math. Prog. Comp. 2018;(10):457-486.
29. Indrapriyadarsini S., Mahboubi S., Ninomiya H., Kamio T., Asai H. Accelerating Symmetric Rank-1 Quasi-Newton Method with Nesterov's Gradient for Training Neural Networks. Algorithms. 2022;(15(6): 1-16.
30. Ivanov V.N. Osnovnye svojstva obratnogo it-eracionnogo algoritma resheniya sistem linejnyh uravnenij s polozhitel'no opredelen-nymi matricami. Vychislitel'nye metody i programmirovanie. 2012;(13(2):366-376. (In Russ).
31. Gill F., Myurrej U., Rajt M. Prakticheskaya optimizaciya. M.: Mir; 1985. 509 s. (In Russ).
32. Minu M. Matematicheskoe programmiro-vanie. Teoriya i algoritmy. M.: Nauka; 1990. 488 s. (In Russ.).
Информация об авторе:
В. Н. Иванов - кандидат физико-математических наук, доцент, доцент кафедры высшей математики механико-математического факультета Пермского государственного национального исследовательского университета (614068, г. Пермь, ул. Букирева, 15), AuthorlD 13318.
Information about the author:
V. N. Ivanov - Candidate of Physical and Mathematical Sciences, Associate Professor, Associate Professor of the Department of Higher Mathematics, Faculty of Mechanics and Mathematics, Perm State University (15, Bukireva Street, Perm, Russia, 614068), AuthorlD 13318.
