CC BY
966
Обзор методов решения обобщённой проблемы собственных значений плотной матрицы, элементы которой нелинейно зависят от параметра
Баскаков А.Е. , Борисов Н.И.
ИТАС, МИЭМ [email protected]
Постановка задачи
Обобщённой проблемой собственных значений нелинейного оператора Г0: С -г £пНп называется задача вычисления пар (Л*,я*) £ £ X £Т1\{0), удовлетворяющих уравнению
П1*К = о- (1)
Как и в случае стандартной проблемы собственных значений А* £ £ называется собственным значением (1), х% - правым собственным вектором, отвечающим А*. Аналогично определяется левый собственный вектор. Множество всех собственных значений Т(■) называется спектром. Алгебраическая кратность собственного значения Аопределяется алгебраической кратностью нулевого собственного значения обыкновенной проблемы собственных значений Т(Х)х = цх. Под размером проблемы понимается порядок матрицы Т(Х). Под Т (А) понимается —
' ' ' СТА
В данном обзоре рассматриваются методы решения обобщённой проблемы собственных значений (далее ОПСЗ) (1) плотных несимметричных матриц, элементы которых нелинейно зависят от параметра. Важно отметить, что методы решения плотных проблем широко используются в проекционных методах решения ОПСЗ с разреженными матрицами большого порядка.
Несмотря на значительный объём работ по методам решения ОПСЗ, существует большое количество прикладных задач, для решения которых возможностей современных методов недостаточно. Применяемые методы зачастую не соответствуют уровню сложности современных задач. Поэтому, как отмечается в [4], требуется глубокое исследование в данной области.
Для решения обыкновенной проблемы собственных значений Ах = Азе или ОПСЗ вида Ах = КЕх доступны многие, хорошо зарекомендовавшие себя методы, включающие в себя оценки точности и обусловленности собственных значений и собственных векторов. Эти методы способны адекватно решать большинство задач, возникающих на практике - от небольших до огромных размеров [5, 6, 7, 8, 9].
Однако для более общих случаев проблемы (1) не существует методов, сравнимых по эффективности и возможности оценки точности решений с методами решения обыкновенной поблемы собственных значений [4].
Возникновение обобщённой проблемы собственных значений в приложениях
Обобщённая проблема собственных значений (1) возникает во многих задачах из различных предметных областей. Теоретические и практические результаты, изложенные в различных статьях зачастую получают рассматривая конкретный вид и
вытекающие свойства Г (Л). Наиболее широко изученным вариантом ОПСЗ является проблема в котором Т зависит от А квадратично
(Я2М + АС + К)х = 0. (2)
где М, С, К £ М*. Такие случаи широко распространены в области динамического анализа механических конструкций [1]. Здесь обычно матрица жёсткости К, матрица масс М являются симметричными положительно определёнными (полуопределёнными). Ещё одним важным источником проблем такого вида является моделирование вибраций вращающихся конструкций [2], где К = Кт, М = Мт - положительно определённые, а матрица демпфирования € = —Ст является действительной кососимметричной матрицей. В большинстве указанных приложений требуется вычислить собственное значение с наименьшей действительной частью.
Другой недавно изученной ОПСЗ является проблема, возникающая в задаче оптимизации акустической эмиссии высоко-скоростных поездов [13, 14]. Модель вибраций рельс приводит к проблеме
^АМХ (со) + М0(со) +1 М[ (со)'} х = 0, (3)
в которой М(ь М1 являются большими разреженными матрицами, зависящими от частоты возбуждающих воздействий со. Здесь М1 (й>) - матрица неполного ранга,
М(>(й>) - комплексная симметричная. Собственные значения возникают парами Л,-,
X
причём большинство сосредоточено в окрестности 0 и эо. Требуется вычислить конечные ненулевые собственные значения и отвечающие им собственные векторы. Важно отметить, что ни один из классических методов не сработал и лишь специальные методы, учитывающие вид Т(А) для данной задачи смогли дать достаточно точные приближения к собственным значениям.
Существуют и другие важные приложения, в которых Т зависит от параметра А квадратично. Подробный обзор такого вида проблем наряду с методами их решения дан в [3].
В более общем случаее зависимость Г от А является полиномиальной
Г®* = (Е5Ц> = 0. (4)
Группы методов решения обобщённой проблемы собственных значений
Методы решения ОПСЗ можно разделить на 3 группы:
1. Методы линеаризации для ОПСЗ с
полиномиальной или рациональной матрицей. В результате линеаризации получают ОПСЗ с пучком матриц большого порядка. Такие методы нацелены на поиск небольшого количества интересующих собственных значений и собственных векторов.
2. Методы уточнения. В данный класс
попадает большинство существующий методов решения ОПСЗ (метод Кублановской, метод обратной итерации и его разновидности и др.). Суть методов состоит в последовательном поиске собственного значения из начального приближения.
3. Методы отделения корней. Сюда входят
методы, основанные на поиске корней детерминантного уравнения. Методы данной группы используются для поиска приближений к
собственным значениям.
Линеаризация
Теоретический анализ и численные методы решения ОПСЗ с полиномиальной или рациональной матрицей основаны на идее линеаризации, то есть на приведении (4) к обобщённой проблеме
(А ~ КЕ)х = 0, (6)
имеющей больший порядок по сравнению с исходной проблемой [ 14]. Собственные значения и собственные векторы (4) вычисляются через (6).
Недостатком такого подхода является значительное увеличение размера проблемы, а так же возможное увеличение числа обусловленности. Это означает, что собственные значения (6) могут быть более чувствительны к изменениям матриц проблемы, чем собственные значения исходной проблемы [10]. Важно так же отметить, что при применении классических методов линеаризации симметричность матриц А и Е в (6) не следует из симметричности матриц М= 0 ...к (4). Разработке новых методов линеаризации, отражающих свойства матриц исходной проблемы посвящены работы [11, 12, 13].
Наряду с методами линеаризации, сохраняющими свойства матриц исходной проблемы, были разработаны методы решения (6), использующие эти свойства. Исследование такого рода методов проводилось многие годы [15, 16] и продолжается до сих пор.
Как показывает практика, линеаризация является на сегодняшний день наиболее простым и широко используемым подходом к решению ОПСЗ с полиномиальной или рациональной матрицей.
Методы уточнения собственных значений
Метод Кублановской
Для решения ОПСЗ в [20] было предложено использовать ^-разложение. Пусть
адред = еотад '
такое разложение, где F(Л) - матрица перестановки, выбираемая таким образом, что диагональные элементы Г)ДА} матрицы Жд) убывают по модулю, т.е. ки(А)| > |г22(Л)| > ••• > \гппт Так как
то А является собственным значением тогда и только тогда, когда
гпп(Х) = 0. (7)
Применяя метод Ньютона к (7) можно получить формулу для вычисления следующего приближения к собственному значению
4+1 = 4 -
Приближениями к левым и правым собственным векторам являются последние г столбцов матриц (£(/.) и /?(А)й(А)-1 соответственно.
Улучшенная версия, данного алгоритма была представлена в [21], где так же была доказана квадратичная скорость сходимости. Трудоёмкость метода достаточно велика - несколько С1 (к3) разложений для вычисления одного собственного значения. Несмотря на это он часто применяется для уточнения найденных собственных значений.
Важно отметить, что на практике скорость сходимости метода Кублановской является приемлемой только в случае наличия простого собственного значения. С целью преодоления этого недостатка в [22] был предложен метод, аналогичный
- 225 -
методу Кублановской, но требующий меньшее количество итераций для вычисления кратного собственного значения.
Метод обратной итерации
Суть метода обратной итерации состоит в применении метода Ньютона для решения системы нелинейных алгебраических уравнений относительно х и Я
(л -1) = (8)
где V - некоторый нормирующий вектор, выбираемый таким образом, чтобы = 1. Применив к (8) шаг Ньютона легко получить формулу для итерационного вычисления х:
XX ЛгттттттлоО (1 ТГТЛ ^, . . =
Умножив (9) на щ и учитывая, что получим
з - 3
лк-1-1 “ Лк ~ н». '
где 11^+1 = Т(Ак)~1Т' (Ак)хк ~ направление к новому приближению.
Метод имеет локальную квадратичную сходимость [17].
Разновидностью метода обратной итерации является метод, основанный на формуле Ланкастера
3 _з УьШкУхк ППЧ
*+1 к
получаемой в результате применения шага Ньютона для решения уравнения
йдам = о
относительно А. Данная формула была предложена Ланкастером [1, стр. 71] для решения ОПСЗ с полиномиальной Т(Х).
Узким местом метода обратной итерации является необходимость вычисления большого количества разложений матрицы Т (Я) при решении СЛАУ 'ик+1 = 'ПА\)~1С целью снижения трудоёмкости метода обратной итерации Ньюмайером [18] был предложен метод остаточной обратной итерации, позволяющий провести несколько или все итерации с использованием одного разложения.
Предположим, что в методе обратной итерации Т7к — г?, а Т(Я) дважды непрерывно дифференцируема. Тогда
<1х = хк - хк.+1 = + 0((лк+1 - лку).
Отбросив второе слагаемое, получим
хк+1 = хк~ Т(Ак)~1Т(Ак+х)хк‘ (11)
Ньюмайер показал, что в (11) вместо Т(Л^) можно использовать 1"(Лр}, где А(> остаётся постоянным в процессе работы метода
**+1 = хк- Т(Л0)~1Т(Лк+1)хк.
На практике, однако, при падении скорости сходимости рекомендуется выбирать новое значение Л&, более близкое к собственному значению.
Скорость сходимости метода определяется следующей формулой. Если А -
/*' А и и и
простое собственное значение, х - отвечающий ему собственный вектор, такой, что Vйх = 1, то для А(> достаточно близкого к А
= 0(д» - Л), |Лк+1 - Я| = 0(\хк -х\ч
где 1 = 2, если Г (Л ) Эрмитова и Л £ К, иначае 1=1.
Метод последовательных линейных приближений
Метод последовательных линейных приближений был предложен в [19]. Пусть Afe - приближение к собственному значению, тогда если Т(А^) является дважды непрерывно дифференцируемой, то
П4 + м) = Ш*) + /<Г'(ад + (12)
где < sup £|«; ft |Г”(ДЬ + 01 Д™ решения исходной ОПСЗ требуется
найти такие х и fj., что
rtffe + = а (13)
Отбрасывая й(Лй, ft) в (12) и учитывая (13), получаем ОПСЗ
TQk}x = -hTXWx, которая решается одним из стандартных методов. Далее берётся
Afe-t-i = 4 +
Метод имеет квадратичную скорость сходимости к простому собственному значению, если Т(А) дважды непрерывно дифференцируема.
Методы отделения корней
Ещё одним подходом к решению (1) является поиск корней детерминантного уравнения
/(А) в det!F(A) = 0, (14)
которые и являются конечными собственными значениями Г(Д}. Существует множество методов поиска корней функций комплексного переменного. Их принято разделять на методы нулевого порядка, использующие лишь значения функции в некоторых точках, методы первого, второго и более высокого порядка, требующие так же вычисления значений производной к-ото порядка /(А). Важной особенностью методов данной группы является теоретически обоснованная возможность поиска нескольких или всех собственных значений.
Значение /(А) при заданном А обычно вычисляется через нормализованное QR-разложение Т(_А). Пусть дано такое разложение
ТШ№ = Q(A}Д(Д), (15)
где /?(А) верхне-треугольная, Ж/1) - ортогональная, причём элементы Ту (А)
матрицы Й(Д) удовлетворяют неравенствам
\гцЩ > |г22(Д)| > - > кп(А)|,|гДЛ)| < |гй(Д)| для/ > i,i = 1,2,... ,п - 1. Тогда
Ключевыми аспектами при анализе алгоритмов поиска корней (14) являются:
1. Предельная сорость сходимости или
асимптотическое поведение скорости сходимости.
2. Критерий достижения корня.
3. Исчерпывание найденных корней.
4. Влияние ошибок округления и ограниченной точности вычислений на указанные выше аспекты.
Скорость сходимости методов сильно варьируется в зависимости от наличия кратных или патологически близких корней. Важно так же, что различные методы ведут себя по-разному в зависимости от того насколько начальное приближение близко к собственному значению. По указанным причинам трудно гарантировать
эффективность того или иного метода для решения широкого класса задач. На практике вместо использования одного метода применяют комбинации методов. Более того, методы, входящие в комбинацию могут быть заменены другими на этапе вычислений.
Подробное рассмотрение всех указанных выше аспектов выходит за рамки данного обзора и может быть найдено в [25].
Методы последовательной линейной и квадратичной интерполяции
Наиболее известными методами нулевого порядка, применяемыми для решения (14), являются методы последовательной линейной и квадратичной интерполяции (см. [25, стр. 435]). В таких методах очередное приближение к корню f(X) находится ПО к предыдущим приближениям Д;_и..., как корень полинома степени к, проходящего через точки (А*л/(А*)),:/(Л_ 1)),...Из корней такого полинома выбирается тот, который расположен ближе всего к А;. Обозначим
л=тг
В случае линейной интерполяции к = 1 и следующее приближение к собственному значению вычисляется по формуле
1 _ 1 1 ■л;—1)
*+1"* + _л^Г
Метод квадратичной интерполяции предполагает к =2. Мюллер [24] предложил эффективный метод вычисления Х(+1.
Недостатком метода линейной интерполяции является невозможность вычисления комплексных корней в случае если £ К и /(А) £ К. при А £ К.
Метод квадратичной интерполяции наоборот может заходить в комплексную плоскость даже при действительных начальных приближениях.
Метода Мюллера для решения обобщённой проблемы собственных значений с использованием формулы Ланкастера (26) для поиска нескольких корней рассмотрен в [23].
Интерполяционные методы порядка 3 и выше применяются редко, так как они требуют решения полиномиального уравнения для вычисления следующего приближения к корню.
Метод Ньютона
Метод Ньютона - классический представитель класса методов первого порядка. Применяя шаг Ньютона к (14), получим
э э
1+1 1
Вычисление ~гг~т возможно проводить по формуле Ланкастера [26].
/ & О
<р сад = уЩ) = гтасе{таг1га)), (I8)
где А не является собственным значением Т(А). При построении методов второго порядка требует вычисления <р-! (А;.). Можно показать, что
<р' (Я,) = гтасе{Т(Х)-1Т"т - (Г(Л)-1Г'(Л))2}. (19)
Вычисление /'(А) зачастую сравнимо по трудоёмкости с вычислением /(А). На основании такого предположения Уилкинсон [25] проводит сравнение предельной скорости сходимости метода Ньютона и методов квадратичной и линейной
интерполяции и показывает, что метод Ньютона проигрывает им как в случае кратного, так и в случае простого корня.
Методы второго порядка
Среди методов второго порядка можно выделить метод Халлея [27]
4+1 = 4 - фЩУ-'р'т/гид ^
и Лагерра [26] ^ ^
где (р{4) = Знак в знаменателе (21) выбирается так, чтобы
" Г
минимизировать \А( — 1. Метод Лаггера имеет кубическую сходимость к простому
корню для любого действительного К > 1.
Метод Лагерра имеет кубическую скорость сходимости к простому корню и широко используется на практике. Наиболее строгое исследование эффективности этого метода проведено Парлеттом [29].
Поиск нескольких собственных значений
Для предотвращения сходимости к уже найденным корням уравнения (14) обычно используют подход, изложенный Ланкастером [1]. Если найдены к корней А11... ,4? т0 производится замен /(Я) на новую функцию /(А)
Л» - <2б)
Процесс такого исключения найденных корней называется исчерпыванием. На практике полученные корни содержат ошибки, и в процессе исчерпывания возникают дополнительные ошибки. Как отмечает Уилкинсон, это может привести к катастрофическим результатам. В частности, если корни полинома сильно отличаются по абсолютной величине и найдены приближения для одного из больших корней, то проведение исчерпывания перед нахождением меньших корней может привести к значительной потере точности. Детальный анализ этой проблемы проведён в [28].
В случае, если производится решение ОПСЗ методами, основанными не на решении детерминантного уравнения так же имеется возможность поиска нескольких собственных значений, однако, в таком случае не существует стратегии, гарантирующей отсутствие сходимости к уже найденным корням и производящей вычисление значительной части спектра.
Один из подходов предполагает использование в качестве начальных приближений решения обобщённой проблемы собственных значений Ах = цВх, возникающей в результате линеаризации Г (Л) в окрестности А(>
[П4) + <Д - 4)П4)1* = о
, где А = 7Т(Л&), В = —Т*(к(у). После того как найдено наименьшее по модулю собственное значение данной проблемы приближение к следующему собственному значению Г(Л) вычисляется согласно формуле 4 = 4 + №■ Для собственных значений не лежащих в окрестности Л^ требуется заново решить линеаризованную ОПСЗ для получения начальных приближений. Такой подход, однако, не гарантирует того, что какие-либо собственные значения не будут найдены повторно.
Другой стратегией, позволяющей предотвратить сходимость к уже найденному собственному значению, является преобразование проблемы с целью сдвинуть
найденные собственные значения в не рассматриваемые области, например к нулю или бесконечности. Преобразование Т(Л) н» Т(Х), при котором
л(г(о) = д (?(•)) \{,уи {ос}, где Ль - уже найденной собственное значение ТО) может рассматриваться как один из видов исчерпывания [30].
Список литературы
1. P. Lancaster. Lambda-matrices and Vibrating Systems. Dover Publications. Mineola, New York, second edition, 2002.
2. R. J. Duffin. The rayleigh-ritz method for dissipative and gyroscopic systems. Quart. Appl. Math., 18:215-221, 19б0.
3. F. Tisseur and K. Meerbergen. The quadratic eigenvalue problem. SIAM Rev., 43:234--28б, 2001.
4. V. Mehrmann and H. Voss. Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods. GAMM-Mitteilungen, 27:121-152, 2004.
5. E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, and D. Sorensen. LAPACK Users' Guide. SIAM, Philadelphia, PA, second edition, 1995.
6. Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst. Templates for the Solution of Algebraic Eigenvalue Problems. SIAM Publications, Philadelphia, 2000.
7. L.S. Blackford, J. Choi, A. Cleary, E.D. Azevedo, J. Demmel, I. Dhillon, J. Dongarra, G. Henry, A. Petitet, K. Stanley, D. Walker, and R. Whaley. ScaLAPACK Users' Guide. SIAM, Philadelphia, PA, second edition, 1997.
8. R. B. Lehoucq, D. C. Sorensen, and C. Yang. ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM, Philadelphia, 1998.
9. MATLAB, Version 7.0. The MathWorks, inc., 24 Prime Park Way, Natick, MA 017б0-1500, USA, 199б.
10. F. Tisseur. Backward error analysis of polynomial eigenvalue problems. Linear Algebra Appl., 309:339--3б1, 2000.
11. D.S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Vector spaces of linearizations for matrix polynomials. Technical report, TU Berlin, Inst. f. Mathematik, Berlin, Germany, 2004.
12. N. Higham, D.S. Mackey, N. Mackey, and F. Tisseur. Linearizations of symmetric polynomials by symmetric pencils. Technical report, Univ. of Manchester, 2004.
13. D.S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Palindromic polynomial eigenvalue problems: Good vibrations from good linearizations. Technical report, TU Berlin, Inst. f. Mathematik, Berlin, Germany, 2004.
14. Y. Su and Z. Bai. Solving Rational Eigenvalue Problems via Linearization. Technical report, UC Davis, 2008.
15. H. Fafi bender. Symplectic Methods for the Symplectic Eigenvalue Problem. Kluwer Academic, New York, 2000.
16. V. Mehrmann. The Autonomous Linear Quadratic Control Problem, Theory and Numerical Solution. Number ІбЗ in Lecture Notes in Control and Information Sciences. Springer-Verlag, Heidelberg, July 1991.
17. P.M. Anselone and L.B. Rall. The solution of characteristic value-vector problems by Newton's method. Numer. Math., 11:38—45, 19б8.
18. A. Neumaier. Residual inverse iteration for the nonlinear eigenvalue problem. SIAM J. Numer. Anal., 22:914-923, 1985.
19. A. Ruhe. Algorithms for the nonlinear eigenvalue problem. SIAM J. Numer. Anal., 10:б74-б89, 1973.
20. V.N. Kublanovskaya. On an application of Newton's method to the determination of eigenvalues of Я-matrices. Dokl. Akad. Nank. SSR, 188:1240 — 1241, 1969.
21. N. K. Jain and K. Singhal. On Kublanovskaya's approach to the solution of the
generalized latent value problem for functional Я-matrices. SIAM J. Numer. Anal., 20:1062-1070, 1983.
22. R.-C. Li. Compute multiple nonlinear eigenvalues. J. Comput. Math., 10:1-20,
1992.
23. Т. Я. Конькова, В. Н. Кублановская, Л. Т. Савинова, “К решению
нелинейной спектральной задачи для матрицы”, Численные методы и автоматическое программирование, Зап. научн. сем. ЛОМИ, 58, Изд-во “Наука” , Ленинград. отд., Л., 197б, 54-бб.
24. Muller, David E. A method for solving algebraic equations using an automatic computer. MTAC, 10 (195б), 208-215.
25. J.H. Wilkinson. The Algebraic Eigenvalue Problem. Oxford University Press, New York, 1988.
26. P. Lancaster. Lambda-Matrices and Vibrating Systems. Pergamon Press, Oxford, 19бб.
27. W. Gander. On Halley's iteration method. Am. Math. Monthly, 92:131—134,
1985.
28. J.H. Wilkinson. Rounding errors in algebraic processes. Dover Publications, Incorporated, 1994.
29. B. Parlett. Laguerr's method, applied to the matrix eigenvalue problem. Math Computation, 18:4б4-485, 19б4.
30. J.-S. Guo, W.-W. Lin, and C.-S. Wang. Numerical solutions for large sparse quadratic eigenvalue problems. Linear Algebra Appl., 225:57-89, 1995.
