Содержание к диссертации
Введение
Глава 1. Матричные уравнения 11
1. Матричное уравнение Сильвестра 11
2. Матричное уравнение Стсйна 16
3. Квадратичные матричные уравнения 17
Выводы 21
Глава 2. Матричные уравнения типа Сильвестра и Стейна. Условия разрешимости и численные алгоритмы 22
1. Уравнения типа Сильвестра 23
2. Уравнения, сопряженные уравнениям типа Сильвестра 51
3. Уравнения типа Стейна 66
Выводы 82
Глава 3. Решение матричных уравнений в самосопряженном случае 84
1. Уравнения Сильвестра и Стейна 85
2. Уравнения типа Сильвестра 90
3. Уравнения типа Стейна 94
Выводы 108
Глава 4. Квадратичные и полуторалинейные матричные уравнения 109
Выводы 124
Заключение 125
Список литературы
- Матричное уравнение Стсйна
- Уравнения, сопряженные уравнениям типа Сильвестра
- Уравнения типа Сильвестра
- Уравнения типа Стейна
Матричное уравнение Стсйна
Матричные коэффициенты А и В уравнения Сильвестра АХ + ХВ = С (1.1) являются квадратными матрицами, вообще говоря, различных порядков тип. Матрицы X и С имеют размер т х п.
Условия однозначной разрешимости. Как и для всякого линейного уравнения, условия однозначной разрешимости уравнения (1-1) эквивалентны условиям, при которых однородное уравнение
Мы заменили матричное уравнение (1.2) уравнением (1.6) того же вида, но в котором матричные коэффициенты имеют жорданову форму. В соответствии с квазидиагональным видом матриц А и В разобьем матрицу Y на блоки: Y = (Yal3) (a = l,2,...,«;/3 = l,2,...,v) (здесь Ya/3 — прямоугольная матрица размером ра х qp). Используя это разбиение, произведем умножение матриц в левой части уравнения (1.6). Тогда это уравнение распадется нага матричных уравнений (XJPa + JPa{0))Y«P + УарЫяр + V0)) = где 1ра единичная матрица порядка ра: a JPo(Q) — жорданова клетка, отвечающая собственному значению 0. Перепишем эти уравнения следующим образом: (Ла + iJLp)Yap = -JpMYctP - Yaf}Jqp(0). (1.7) Возьмем какое-нибудь из уравнений (1.7). Возможны два случая: 1. \а + 1Л/з т 0. Проитерируем равенство (1.7) (обе части равенства (1.7) умножаем на Аа + /і# и в каждом члене правой части заменяем (\а + Иіз)УаіЗ
Если здесь взять г ра + Q/з - 1, то в каждом члене суммы, стоящей справа, выполняется по крайней мерс одно из соотношений а ра: т qp, &, значит, и по крайней мере одно из соответствующих соотношений Jpa(0) = 0, JL(Q) = 0. Следовательно, правая часть равна нулю, а потому это уравнение допускает единственное решение Yap = 0.
Нетрудно видеть, что это уравнение допускает нетривиальные решения: таковым будет, например, всякая матрица, у которой отличен от нуля лишь один угловой элемент в позиции (l,qp). Таким образом получены условия однозначной разрешимости уравнения Сильвестра: Теорема 1 Уравнение Сильвестра однозначно разрешимо для любой правой части С, если
Численный алгоритм. Опишем теперь один из эффективных численных алгоритмов для решения уравнения Сильвестра — алгоритм Бартелса-Стьюарта. Предполагается, что выполнены условия однозначной разрешимости, формулируемые теоремой 1. Алгоритм состоит из следующих этапов:
Приведение матриц А и В к форме Шура. Для определенности, будем использовать верхнюю форму Шура матрицы А и нижнюю форму матрицы В. Тогда содержанием данного этапа является вычисление унитарных матриц U и V таких, что матрица нижнетреугольная. Средством вычисления матриц (1.9) и (1.10) может служить, например, комплексный QR-алгоритм, интерпретируемый не как метод нахождения собственных значений, а именно как метод приведения к треугольному виду (см. об этом подробнее в [2, 4]).
Решение уравнения (1.11). Это матричное уравнение представляет собой треугольную систему линейных уравнений относительно тп коэффициентов yij матрицы Y. Вот как решается эта система: приравнивая в (1.11) элементы в позиции (т,п): имеем
Элементы атт и Ьпп суть собственные значения соответственно матриц А и В, поэтому атт + Ьпп 0 и уравнение (1.13) определяет коэффициент утп- Возвращаясь к уравнению (1.11), приравниваем элементы в позиции (т,п — 1):
Продолжая таким образом, определим все коэффициенты последней строки матрицы У. После этого вычисляем коэффициенты предпоследней строки последовательным приравниванием в (1.11) элементов в позициях (т —1, j), где j, убывая, меняется от п до 1; затем переходим к вычислению коэффициентов (т — 2)-й строки, и т. д. Если вычислены все элементы yij: для которых і к, и все г/fcj, для которых j /, то формула для вычисления элемента у и выглядит следующим образом:
Возврат к исходной матрице X. Обращая соотношение (1.12), получаем X = UYV . Сложность описанного алгоритма в случае, когда все матрицы уравнения (1.1) квадратные одного и того же порядка п, находится из формулы (1.14) и составляет 0(п ) арифметических операций.
Матричное уравнение Стейна Матричные коэффициенты А и В уравнения Стейна Х + АХВ = С\ (2.1) как и у уравнения Сильвестра, являются квадратными матрицами, вообще говоря, различных порядков т и п, а матрицы X и С имеют размер тхп. Условия однозначной разрешимости Теорема 2 Уравнение Стейна однозначно разрешимо для любой правой части С, если \i(A)\j(B) ф -1 Vi,j. (2.2) Численный алгоритм. Опишем алгоритм Голуба-Нэша-Ван Лоана для уравнения Стейна. Как и алгоритм Бартелса-Стьюарта, он состоит из четырех этапов: 1. Приведение матрицы А к верхней форме Хессенберга, а матрицы В к нижней форме Шура. А = U AU, В = V BV. (2.3) Форма Хессенберга может быть вычислена последовательностью т — 2 подобий с отражениями в качестве трансформирующих матриц (см. [2, 3]). 2. Преобразование правой части. На этом этапе вычисляется матрица
Итак, установлено соответствие между решениями уравнения (3.1) и n-мерными инвариантными подпространствами Сn матрицы М такими, что в любой их базисной матрице верхняя п х п-подматрица невырожденна. Специальный выбор базисной матрицы (3.12) показывает, что это соответствие взаимно однозначно: различным решениям отвечают различные инвариантные подпространства, и обратно.
Численный алгоритм. Опишем алгоритм построения какого-либо решения уравнения (3.1). Алгоритм состоит из трех этапов. 1. Составление матрицы М. По формуле (3.2) из коэффициентов исходного уравнения строится матрица М. 2. Приведение матрицы М к форме Шура. Это достигается QR алгоритмом. Пусть U — унитарная матрица, трансформирующая М к верхней форме Шура S:
Уравнения, сопряженные уравнениям типа Сильвестра
По предположению, пара (Ai,Bi) удовлетворяет условиям теоремы 1. Это означает, в частности, что хотя бы одна матрица в указанной паре невырожденна. Пусть, например, невырожденна матрица В\. Тогда (1.33) можно переписать в виде уравнения Стейна: Все собственные значения жордановой клетки j +1 равны нулю. Таким образом, выполнены условия теоремы 2 главы 1 и единственным решением уравнения (1.34) является Удп = 0. Но тогда и Y\N = 0. Пусть теперь невырожденна матрица А\. В этом случае уравнению (1.33) можно придать вид уравнения Сильвестра
При вырожденной матрице В\ у этого уравнения существуют ненулевые решения, так как условия однозначной разрешимости из теоремы 1 главы 1 нарушены наличием нулевого собственного значения у обоих матричных коэффициентов. Однако теперь важную роль приобретает соотношение A\z = 0, которому подчинен последний столбец z матрицы Y/vi- Поскольку А\ невырожденна, то последний столбец искомого решения Y/vi должен быть нулевым.
Разобьем матрицу W на полосы, ширина которых соответствует порядкам жордановых клеток в матрице J. Согласно [1, гл. VIII, 2], те полосы, которые отвечают ненулевым собственным значениям в J, целиком состоят из нулей. Однако в общем решении W уравнения (1.37) полосы, соответствующие нулевому собственному значению матрицы J, могут быть ненулевыми. Ненулевая зона каждой такой полосы ширины s есть верхнетреугольная теплицева матрица порядка s, стоящая в последних s столбцах. Все прочие элементы данной полосы равны нулю.
Так как последний столбец искомого решения W нулевой, то нулевыми являются и описанные выше теплицевы матрицы. Тем самым W = 0, а потому (см. (1.36)) Удгі = 0. Но тогда и Y\N = 0.
Итак, все блоки в последних блочной строке и блочном столбце матрицы Y нулевые. Поэтому Y = 0, и теорема 3 доказана.
Теорема 4 Пусть в уравнении матричные коэффициенты А и В получены из матриц (1.29) следующим образом: А есть результат окаймления матрицы А сверху t нулевыми строками; приписывание к В слева t нулевых столбцов дает матрицу В. Если регулярный пучок А1 + ХВ , содержащийся в матрицах (1.29), удовлетворяет условиям, теоремы 1, то уравнение (1.38) имеет тюлько нулевое решение.
Согласно теореме 3, уравнение (1.39) имеет только тривиальное решение І2 = 0. Отсюда следует, в частности, что правое ядро пучка А + ХВ состоит только из нулевого вектора. Это означает, что общим решением уравнений (1.40) может быть лишь нулевая матрица Y1. Итак, Y = 0, что и требовалось доказать.
Теорема 4 дает (в терминах канонической формы пучка А + ХВ ) наиболее общие условия однозначной разрешимости уравнения (1.1) при любой правой части С. Численный алгоритм. Пусть в уравнении (1.1) все матрицы квадратные.
Средством доказательства теоремы 1 о единственности решения уравнения (1.1) является каноническая форма Вейерштрасса пучка (1.2). Будучи удобным теоретическим инструментом, эта форма мало пригодна в практических расчетах. Хорошо известно, что задача вычисления формы Вейерштрасса (как и ее частный случай — вычисление жордановой формы квадратной матрицы) численно неустойчива. Этот раздел посвящен именно практическому решению уравнений вида (1.1). Предложен прямой ортогональный алгоритм для решения таких уравнений, который можно рассматривать как аналог алгоритма Бартелса-Стьюарта для уравнений Сильвестра, описанного в 1 главы 1. Применением QZ-алгоритма исходное уравнение приводится к уравнению того же типа с треугольными матричными коэффициентами А и В. Полученное матричное уравнение эквивалентно последовательности систем линейных уравнений малого порядка относительно коэффициентов искомого решения. Посредством численных примеров моделируется ситуация, когда "почти" нарушены условия однозначной разрешимости. Прослежено ухудшение качества вычисленного решения в этой ситуации.
Предположим, что матричный пучок (1.2), связанный с уравнением (1.1), регулярен и выполнены условия единственности решения из теоремы 1. Предлагаемый алгоритм для решения уравнения (1.1) имеет ту же структуру, что и алгоритм Бартелса-Стьюарта.
Приведение матриц А и В к форме Шура. Для определенности, будем приводить пучок (1.2) к верхней форме Шура. Тогда содержанием данного этапа является вычисление унитарных матриц U и V таких, что матрицы
Коэффициенты уп,п-2 И 2/п-2,п Продолжая таким образом, находим все коэффициенты последней строки и последнего столбца матрицы Y. После этого вычисляем коэффициенты предпоследней строки и предпоследнего столбца (кроме уже известных коэффициентов уп,п-1 и Уп-1,п) последовательным приравниванием в (1.42) элементов сначала в позиции (п — 1, п — 1), а затем в парах позиций (п — 1, j) и (j, п — 1), где j, убывая, меняется от п — 2 до 1. Вслед за этим переходим к вычислению коэффициентов (п — 2)-й строки и (п — 2)-го столбца, и т. д. Пусть вычислены все у , для которых либо і / либо j /, и все г/j/ и у/j, для которых і к (к I). Тогда коэффициенты уik и уи вычисляются из следующей системы
Только что описанный алгоритм ABST (Algorithm of the Bartels-Stewart Type) был реализован в виде функций BST и BSTstar языка Matlab. Обсуждаемые ниже эксперименты использовали именно эти функции.
Поскольку этапы 1, 2 и 4 алгоритма ABST основаны на ортогональных процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку.
Были проведены две массовые серии расчетов для уравнений десятого порядка. В первой серии коэффициенты А, В и С уравнения (1.1) генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10. Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного решения X оценивалось посредством (евклидовой) нормы невязки
Уравнения типа Сильвестра
Возвращаясь из С в пространство МП(С), построим из блоков вещественной матрицы A Y комплексную матрицу Таким образом, еще раз показано, что в пространстве Мп(С) сопряженным к оператору левого умножения на фиксированную матрицу А является оператор левого умножения на А . различных собственных значения Х{ и Xj этого пучка не удовлетворяют соотношению Теорема 5 Пусть в уравнении (2.6) т = п и пучок (2.13) регулярен. Тогда для однозначной разрешимости этого уравнения при любой правой части С необходимо и достаточно выполнение следующих условий: 1) хотя бы одна из матриц А и В невырожденна; 2) никакое из собственных значений пучка (2.13) не удовлетворяет соотношению за тем исключением, что приТі = Т допускается собственное значение А = 1 кратности 1; 3) никакие два различных собственных значения Х{ и Xj этого пучка не удовлетворяют соотношению В отличие от теоремы 1, в теореме 5 вместо обобщенного пучка используется конкретный пучок (2.13). Связано это именно с тем, что под сопряженным к уравнению (2.6) следует понимать уравнение
Численный алгоритм. В дальнейшем предполагается, что уравнение (2.6) удовлетворяет приведенным выше условиям однозначной разрешимости. Найдем типы преобразований, сохраняющих вид уравнений (2.6). Одно из допустимых преобразований очевидно: умножая обе части (2.6) слева на невырожденную матрицу Q7 приходим к уравнению того же вида
Теперь мы можем приступить к описанию алгоритма численного решения уравнения (2.6), основанного на формулах (2.31). Он имеет структуру, сходную со структурой метода Бартелса-Стьюарта для матричных уравнений Сильвестра, и, подобно последнему, состоит из четырех основных этапов: 1. Приведение матриц А и В к (верхней) треугольной форме. Соответствующие матрицы Р и Q выбираем так, чтобы они были унитарными. Средством их вычисления может быть хорошо известный QZ-алгоритм.
Это уравнение эквивалентно блочно-трсугольной системе уравнений относительно неизвестных yij: причем диагональные блоки матрицы коэффициентов имеют порядки 1 либо 2. Вот как решается эта система: приравнивая в (2.33) элементы в позиции (п,п) и полагая А = (dij)}B = (bij): имеем
Отношение Хп = f есть одно из собственных значений пучка (2.13). Привлекая первое и второе условия теоремы 5, делаем вывод, что уравнение (2.34) однозначно определяет элемент упп.
Определитель этой системы CLn-i n-\CLn,n—bn-i n-i n,n не равен нулю, так как в противном случае было бы нарушено либо первое, либо третье условие А Ап т 1 теоремы 5. Тем самым элементы уп-\,п и уп,п-\ однозначно определены. Переходя к позициям (п, п — 2) и
Продолжая указанным образом, находим все элементы последней строки и последнего столбца матрицы Y. После этого вычисляем элементы предпоследней строки и предпоследнего столбца (кроме уже известных элементов уп,п-1 и Уп-1,п) последовательным приравниванием в (2.33) элементов сначала в позиции (п — 1, п — 1), а затем в парах позиций (п — 1, j) и (j, п — 1), где j, убывая, меняется от п — 2 до 1. Вслед за этим переходим к вычислению коэффициентов (п — 2)-й строки и (п — 2)-го столбца, и т. д. Пусть вычислены все yij: для которых либо і I либо j /, и все г/j/ И у/j, для которых і к (к I). Тогда коэффициенты yik и Возврат к исходной матрице X. Эта матрица выражается соотношением (2.28) при R = Рп. Можно показать, что общая трудоемкость этого алгоритма составляет 0(п ) арифметических операций.
Описанный алгоритм ABST был реализован в виде функций BSTaxbx и BSTaxbxstar (для уравнений (2.4) и (2.5) соответственно) языка Matlab. С этими функциями были проведены эксперименты, аналогичные экспериментам предыдущего параграфа. Проведенные эксперименты подтверждают численную устойчивость алгоритма.
Были проведены две массовые серии расчетов для уравнений десятого порядка. В первой серии коэффициенты А, В и С уравнения (2.6) генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиуса 10. Качество вычисленного решения X оценивалось посредством (евклидовой) нормы невязки
Для уравнений (2.6) с псевдослучайными матричными коэффициентами условия однозначной разрешимости, перечисленные в теореме 5, выполнены, как правило, со значительным запасом. Поэтому никаких проблем с численной устойчивостью в этих двух сериях не было выявлено. Типичное значение нормы невязки (полученное усреднением по 100000 уравнений) составляло в первой серии 1.6545 х 10 12 и 4.2689 х Ю-11 соответственно для функций BSTaxbx и BSTaxbxstar. Для второй серии типичными значениями ошибок (при том же усреднении) были
Еще две серии экспериментов были проведены с тем, чтобы проследить ухудшение качества вычисленного решения по мере приближения к ситуации, когда условия однозначной разрешимости уравнения (2.6) нарушены. В первой из них рассматривается последовательность матриц At и Btl нормализованных таким образом, чтобы одно из собственных значений пучка (2.13) равнялось соответственно для функций BSTaxbx и BSTaxbxstar. Графики eabs и rei как функций от t представлены на рис. 9, 10 (BSTaxbx) и 11, 12 (BSTaxbxstar). Эти графики также получены усреднением по 10 однотипным уравнениям.
Моделированию нарушения третьего условия теоремы 5 посвящена последняя серия экспериментов. Последовательность матриц At и Bt строится таким образом, чтобы для некоторых двух собственных значений ХІ и Xj пучка (2.13) выполнялось соотношение
Размер всех матриц X, Д В и С в уравнениях (3.1) и (3.2) одинаков, в общем случае он равен m х п. Размерности матриц в уравнении (3.3) те же, что и в уравнении Стейна: X и С — матрицы размера m х n, а матрицы А и В — квадратные порядков соответственно тип. Условия однозначной разрешимости уравнений (3.2) и (3.3) найдены в [7]. С этими уравнениями связываются уравнения Стейна, что позволяет сразу же предложить численный алгоритм на основе алгоритма из 2 главы 1. Не так просто обстоят дела с уравнением (3.1): в [7] для него найдены лишь условия разрешимости, условий однозначной разрешимости авторам этой статьи использованным ими методом получить не удалось. И хотя с уравнением (3.1) также связывается уравнение Стейна, говорить о надежности построенного на его основе алгоритма уже не приходится. В этом параграфе мы исправляем ситуацию, выводя условия однозначной разрешимости уравнения (3.1) и предлагая прямой численный алгоритм.
Уравнения типа Стейна
Описанные только что алгоритмы решения самосопряженных уравнений (1.1) и (1.2) были реализованы в виде функций языка Matlab BSss и Stcinss. Мы провели сравнение времени работы этих функций с временем работы стандартных процедур Matlab a lyap и dlyap.
Для этого сравнения был взят случай т = п. Матричные уравнения для наших тестов конструировались следующим образом: матричные коэффициенты А: В и С генерировались как п х n-матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10. Кроме того, при построении матриц А и В учитывались условия теорем 1 и 2.
Для каждого значения п решались десять самосопряженных уравнений типа (1.1), используя вначале lyap, а затем BSss. Время решения, усредненное по этим десяти уравнениям, обозначим через t\{n) в случае lyap и через t2(n) для BSss. Начальное значение порядка п равнялось 50. В дальнейшем п возрастало с шагом 100 до значения 2950. На рис. 1 отношение ti(n)/t2(n) показано как функция от порядка п.
Аналогичные тесты были проведены для самосопряженных уравнений типа (1.2). При этом использовались функции dlyap и Steinss. Усредненное время решения обозначается через t (n) в случае dlyap и через t±{n) для Steinss. На рис. 2 показано отношение tz{n)/t±{n) как функция от п. Это отношение несколько меньше отношения ti(n)/t2(n): что можно объяснить высоким быстродействием функции dlyap. По крайней мере, в наших тестах эта функция всегда работала быстрее, чем lyap. 500 может рассматриваться как система из т2 линейных (или полулинейных) уравнений относительно тп неизвестных Xij. При тфп такая система или недоопределена, или переопределена. В этом случае наиболее правильным подходом к его решению является метод наименьших квадратов. Далее мы рассматриваем случай т = п.
Напомним, что в 1 главы 2 алгоритм для решения уравнения (2.1) был реализован в виде функций языка Matlab BST и BSTstar. Алгоритм для решения самосопряженного уравнения (2.1), описанный выше, таким же образом был реализован в виде функций BSTss и BSTstarss.
В наших экспериментах мы интересовались ускорением, которое можно получить, решая самосопряженное уравнение (2.1) алгоритмом данного параграфа вместо использования алгоритма для уравнения общего вида.
Для самосопряженного уравнения (2.1) были проведены четыре серии расчетов, в которых использовались функции BST, BSTss, BSTstar и BSTstarss. В каждой из серий матрицы А и С уравнения (2.5) генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10 (при этом матрица А строилась как самосопряженная). Порядок п этих матриц возрастал от 50 до 950 с шагом 100. Для каждого из этих значений п было измерено время работы программ BST, BSTss, BSTstar и BSTstarss: i(n), (п), t (n) и t±(n), с усреднением по 10 однотипным уравнениям. Интерес представляют отношения -(Ч и щ, зависимость которых от порядка п показана на рис. 3 и 4.
Хотя все алгоритмы имеют сложность О(п ) арифметических операций, на графиках наблюдается линейная зависимость. Она объясняется неоптимизированностью третьих этапов алгоритмов BST и BSTstar для 100 200 300 400 500 Є0О 700 800 900 1000
В этом параграфе мы предложим для самосопряженных случаев уравнений X + АХТВ = С, X + АХ В = С и X + АХ В = С модифицированные численные алгоритмы. Для последних двух уравнений помимо прямого численного алгоритма возможна модернизация основного алгоритма (алгоритма на основе сведения к уравнению Стейна). По отношению к первому уравнению рассмотрение этого способа мы посчитали нецелесообразным ввиду возможной неоднозначной разрешимости соответствующего уравнения Стейна (см. 3 главы 2).
Напомним, что алгоритмы для решения уравнения (3.1), описанные в 3 главы 2, были реализованы в виде функций языка Matlab BSTxaxb (алгоритм типа Бартелса-Стьюарта) и dlyapT (сведение к уравнению Стейна). Алгоритм для решения самосопряженного уравнения (3.1) таким же образом был реализован в виде функции BSTxaxbss.
Были проведены три серии расчетов, в первой из которых использовалась функция BSTxaxb, во второй — dlyapT, в третьей — BSTxaxbss. В каждой из серий квадратные матрицы Аи С уравнения (3.4) генерировались как матрицы с псевдослучайными элементами, равномерно 100 200 300 400 500 Є00 700 800 900 1000 dlyapT, BSTxaxbss: f-(n) распределенными в круге радиуса 10. Порядок п этих матриц возрастал от 50 с шагом 100. Для каждого из этих значений п было измерено время работы программ BSTxaxb, dlyapT и BSTxaxbss: ti(n), (п) и t (n), с усреднением по 10 однотипным уравнениям. Зависимость отношений jj4 и jj-т от порядка п показана на рисунках о и о.
На рис. 5 наблюдается линейная зависимость, хотя оба алгоритма BSTxaxb и BSTxaxbss имеют сложность 0(п ) арифметических операций. Эта зависимость объясняется так же, как и в аналогичном случае с функциями BST и BSTstar в 2: неоптимизированность третьего этапа алгоритма BSTxaxb для среды Matlab неблагоприятным образом сказывается на времени работы этого метода в несамосопряженном случае.
Теперь, вместо того чтобы решать уравнение (3.16) стандартными алгоритмами, заметим, что его коэффициенты являются эрмитовыми матрицами, а, значит, уравнение (3.16) удовлетворяет условиям теоремы 2. Для решения этого самосопряженного уравнения будем использовать алгоритм, разработанный в 1.
Алгоритм для решения уравнения (3.8), описанный теоремой 6 в 3 главы 2, был реализован в виде функции язвіка Matlab dlyapStar. Алгоритмві для решения самосопряженного уравнения (3.8) бвіли реализованы в виде функций SteinStarss (сведение к уравнению Стейна) и Starss (прямой алгоритм).
Для самосопряженного уравнения (3.8) были проведены три серии расчетов, в первой из которых исполвзоваласв функция dlyapStar, во второй — SteinStarss, в третьей — Starss. Расчеты проводились так же, как в предыдущем разделе. Интерес представляют отношения JJ4 и JJ4, зависимость которых от порядка п показана на рисунках 7 и 8.
