CC BY
69
УЧЕНЫЕ ЗАПИСКИ ЦА Г И
Т о м VI
197 5
№ 2
УДК 518:512.35
629.735.33.015.4:533.6.013.422
ПОЛНАЯ ПРОБЛЕМА СОБСТВЕННЫХ ЗНАЧЕНИЙ МАТРИЦ В РАСЧЕТАХ НА ФЛАТТЕР
В. Г. Буньков
Улучшены известные методы решения полной проблемы собственных значений произвольных матриц, что позволило рассчитывать полный спектр собственных колебаний неконсервативных систем с числом степеней свободы до тридцати. Достигнуты следующие улучшения: упрощены формулы в методе Гивенса, служащего для получения характеристического полинома матрицы; разработан новый метод вычисления корней полинома на основе комбинации известного метода спуска и метода квадратичной интерполяции; по-новому осуществлено определение собственных векторов по найденным значениям корней полинома.
1. Свободные колебания в условиях полета летательного аппарата, идеализируемого как линейная неконсервативная система с г степенями свободы, описываются следующим матричным уравнением:
где (?, С, В, А й* — заданные матрицы г-го порядка [1], 11 — вектор обобщенных координат, V — параметр потока, х — показатель.
Уравнение (1) сводится к алгебраической проблеме собственных значений X матрицы А, имеющей порядок п — 2г:
Комплексные частоты 5 = 8 + ш и формы колебаний и (имеется в виду и = иен) связаны с собственными числами X и векторами X матрицы А следующим образом:
Нормирующий множитель <1>с — характерная частота системы —вводится с целью избежать переполнения при вычислении коэффициентов полинома /(> ) = |Х£" — А |. Для степени р в выражении (2),
<?£/ + С О + Vх В і/ + юОІ! + Б* О = 0,
(1)
определяющем нормирующий множитель с»с, опытным путем найдено наиболее подходящее значение, равное 3.
Приведем два примера из практики расчета колебаний самолетов в полете. Оба расчета проводились по методу многочленов [1, 2].
Первый пример. Расчет самолета, схематично изображенного на фиг. 1, характерен тем, что здесь проверены возможности предлагаемой ниже методики на решении полной проблемы собственных значений матрицы 60-го порядка. Это наибольший порядок
500
400
100
200
100
О
- х-\ \ \ | 3
\ \ \ ¿6 шЙ- из гид 0 руля
1 / У са6
/ > \ \ ър '(¿¡¿-ЗУ теран
ч ^V гг-01 °>3 ш,-6раще-
МН| \ °>1 а>0-,и> 'а
30 -20 -10 О
Фиг. 3
10 &
несимметричных матриц с комплексными собственными значениями, для которых удалось получить с шестью — девятью знаками все корни и векторы. При этом весь расчет на БЭСМ-6 ¡занимал пять минут.
Второй пример. Расчет самолета, изображенного на фиг. 2, содержит расчет комплексных частот колебаний 5 — 8 +/<в в широком диапазоне скоростей полета V. На фиг. 3 показаны годографы комплексных частот при вариации скорости V от нуля. Вариация скорости задается по арифметической прогрессии из 10 -т~ 30 шагов. На каждом шаге решается полная проблема собственных значений для матрицы 24-го порядка (здесь учитывались 12 степеней свободы), что занимает около пяти минут. Читатель, знакомый с теорией флаттера, сможет увидеть на фиг. 3 развитие рулевой и элеронных форм флаттера.
2. Среди различных методов решения полной проблемы собственных значений несимметричных матриц заслуживают внимания I/?- и С?1-алг©ритмы [3], родственные итерационным степенным методам, а также методы, приводящие к раскрытию характеристического полинома матрицы — это методы Гивенса и Хаусхольде-ра [4]. Из двух последних мы отдаем предпочтение методу Гивенса [5], поскольку он оказывается для матриц А специального
вида (2), экономичнее метода Хаусхольдера, чего нельзя сказать в отношении матриц общего вида.
Не описывая целиком метод Гивенса, отметим, что этот метод позволяет путем подобного преобразования исходной матрицы А с помощью некоторой легко определяемой ортогональной матрицы Т привести ее к нижней почти треугольной матрице С:
Двойным рекуррентным формулам для раскрытия характеристического полинома можно придать более изящный вид, чем в оригинале [5], если определить полином не традиционно: \А — ХЕ\ * а наоборот: /(X) = | ХЕ — А \ . Напишем эти формулы:
Вычисление корней полинома /(X) составляет основную часть настоящей статьи и приведено ниже. Здесь же показано, как определяются собственные векторы У по найденным собственным значениям X. Начинать вычисление векторов можно, руководствуясь [5]* т. е. решение однородной системы
В большинстве случаев для простых (некратных) корней и небольших порядков (га<С20) описываемая процедура дает приемлемое решение. При этом, говоря о небольших порядках, мы имеем в виду ЭЦВМ БЭСМ-6 с 40-разрядной мантиссой (относительная точность одиночных операций з — 210 —12). Получение полинома методом Гивенса и большая часть операций в расчете корней полинома должны быть запрограммированы с двойной точностью — только при этом можно говорить о п = 20 как о „небольшом“ порядке и достичь максимально реализуемого в расчетах порядка п = 60.
Обобщим процедуру (3) на случай кратных корней, а также на случай малых разностей, возникающих в процессе вычислений.
АХ = XX]
С—Т АТ\ (Т' Т = Е)\ С = СГ = ХГ; X — ТУ;
си ¿>,0 0. . . 0
с21 с22 Ь2 0 . . . 0
Ъп-1
I £п 1 Сп 2 сп 3 • • • С'пп )
Л = х/й-х + 5*_1; к = 2, 3, . . . , щ /!=Х — си\ Б] = Ь] Sj—\ Сь% у+1 fу = 1, 2, . . . , к 1; 50 = •
получить по рекуррентным формулам:
к
к
(3>
ук+1= — к= 1> • ■ • > л_1-
я — 1 .
Кратность корней определим следующим образом. Расположим корни Хк в порядке возрастания модулей |ХА|^|Хй+1| и обозначим среднеарифметическое всех модулей |ХЙ| через Хср. Будем считать корень ХЛ кратным, если | ХА — ХА+11 <е"Хср, и трехкратным, если |ХЙ—Хй+1| + + | >*+1 — ^А+2 | < в2 К? (пропуская сопряженные корни).
Поскольку, как будет показано ниже, простые корни вычисляются с девятью точными знаками, кратные (двойные)—с четырьмя знаками, а трехкратные — с тремя знаками, назначим для е" и е” величины 10-4 и 10—3 соответственно.
Вернувшись к процедуре (3), обозначим
ёк = шах {| ск}у! |, |Ху*|}.
В случае \bk\^>zlgk (е1 = 210— 12) процесс вычислений идет обычно. При возникают два случая:
1) I ^к I ~ числитель конечный, знаменатель—„нуль“. Про-
цесс может быть продолжен лишь после придания нулевых значений всем только что найденным компонентам: _у} = • • • =_уА = 0 и назначения ук+\ = 1. Для е2 берется либо е'2 — 10-в для простых корней, либо е’=10-4 для кратных, либо е'" = 10—3 для трехкратных корней;
2) —неопределенность „нуль, деленный на нуль“. В этом случае для простого корня оставшаяся часть вектора определяется решением системы
где Ух = (у1, у2, . . ., ук)' — уже найденная часть вектора, У2 = = {Ук+1, • • УпУ — оставшаяся часть вектора, С0 С2 и С12--квадратные и прямоугольный блоки матрицы С — ХЕ. Легко видеть, что оставшаяся часть вектора определяется так:
Но для кратных корней обратная матрица Сг”1, очевидно, не существует: |С2| = 0. Об этом можно судить хотя бы по тому, что для кратного корня собственный вектор неединствен. Поэтому для кратных корней поступаем следующим образом. При первой же встрече с неопределенностью „нуль, деленный на нуль“ назначим в качестве собственного вектора уже найденные компоненты (К!), придав нулевые значения оставшимся (Г2 = 0). Другой собственный вектор будет иметь в качестве первого ненулевого компонента ук+1 = 1, и процедура (3) продолжается. Такой подход возможен только тогда, когда кроме условия | С2 [ = 0 выполняется условие С12У 1 = 0, в противном случае мы имеем так называемый „жорда-новский ящик“, т. е. случай потери одного из собственных векторов. Таким образом, выявление неопределенности „нуль, деленный на нуль“ позволяет получить векторы, линейно-независимые для кратных и тройных корней (для трехкратных корней такая неопределенность должна встретиться дважды).
По-видимому, из-за слишком высокого критерия для диагностики „нуля“ в знаменателе (\bk\^.e1gk, е1 = 210 —12) зачастую
или
Сі У, = о,
Сп К1 + С2У2 = 0,
У2 — — С2 С^2У^.
будут упускаться возможности описываемого выше приема и для корня ХА, диагностируемого как кратный, процедурой (3) определится лишь одни вектор. Но тогда другой вектор легко можно получить, „отменив“ кратность, т. е. находить его для нового значения Х*+1 =^=ХА обычным путем. Например, на практике для простых корней матрицы 48-го порядка, отличающихся лишь на две единицы в четвертом знаке, были получены векторы с такой же хорошей линейной независимостью как и для изолированных корней (такая линейная независимость оценивалась перемножением всех векторов прямой А и транспонированной матрицы А', в результате чего перекрестные произведения оказались на шесть порядков меньше одноименных).
3. Ни процедура (3), ни описанный выше прием не позволяют получить приемлемый результат для векторов матриц среднего и высокого порядка (га = 4060), а для некоторых типов матриц — даже и для небольшого порядка (га~20). В лучшем случае получаем лишь промежуточные данные для дальнейшего уточнения. В связи с этим мы предлагаем алгоритм уточнения собственных векторов, перекликающийся с методом обратной итерации и некоторыми другими приемами из книги Уилкинсона [4].
Пусть относительная ошибка, характеризующая точность вычисленного /-го тона, т. е. совокупности значения Х;. и вектора X¡-.
_ \AXj-\jXjI
М*у1 ’
не удовлетворяет нас, т. е. больше некоторой максимально допустимой величины о': 0у>з'. В качестве примера укажем, что для матрицы 60-го порядка вектор Х} с ошибкой а,—0,02 может не иметь ничего общего с истинным вектором. Зададимся целью уменьшить ошибку до величины а"<о'. Обычно принималось а'=10-3,, о" = Ю-6.
Составим новую матрицу Н, связанную с исходной следующим образом:
Н={А -(1+7)Х,£■]-!,
где Ху — имеющееся приближенное собственное значение, у—малый вещественный добавок (на практике у = 10~3).
Вещественная и мнимая части комплексной матрицы Н — Нх-(-+ Шг могут быть получены по следующем формулам:
а + ¿р = (1 + т)Ху; # = .4—а£;
(Я+ №)-'; я2 = 5Л,.
Эти вычисления можно делать по обычным стандартным программам обращения матриц.
Легко убедиться, что наибольшее собственное число X матрицы Н по модулю на несколько порядков больше всех остальных чисел (если только (1+т) Ху не попало в окрестность Ху+1) и связано с точным значением собственного числа Х;. матрицы А следующим образом:
*7 = ^/0 + Т) + ~ •
Соответствующие собственные векторы матриц А и Н совпадают.
Применяя степенной метод, т. е. произведя итерации:
Zft = // 6=1, 2, . . I,
где в качестве исходной итерации Е0 берется приближенный вектор X] (вообще можно брать произвольный вектор), мы получим сильную СХОДИМОСТЬ отношения компонентов векторов Zk И Zk-1 к собственному значению X, а самого вектора Zk — к уточненному собственному вектору XПопутно уточняется и корень: Ху вместо Х;-. Критерием сходимости служит ошибка а у, которая должна вычисляться на каждой итерации, пока не достигнет а" или пока число итерации не превысит заданное число I (мы считаем, что при / = 8^-10 попытку уточнения следует прекратить, как неудав-шуюся).
Опыт расчетов на флаттер показал, что предложенный метод уточнения собственных векторов дает хорошие результаты вплоть до « = 60. Недостатком его является большое время расчета; например, для матрицы порядка « = 48 (расчет на флаттер с 24 степенями свободы) метод Гивенса с раскрытием полинома занимает 20 с, вычисление корней полинома — 2 с, процедура (3)—1 с, а уточнение каждого вектора — по 2 с. При п = 48 уточнять приходится 14—20 старших тонов (30—40%), при « = 24 — 2 — 4 тона, а при п = 60—практически все. В то же время оказалось, что собственные значения во всех случаях вычисляются с 8—9 знаками и без уточнения.
Интересно, что с уменьшением числа верных знаков в коэффициентах исходных матриц О, С, В, /), точность определения комплексных частот системы стремительно падает. Так, специальными расчетами для системы с 24 степенями свободы было получено, что уменьшение числа верных знаков в матрицах О, С, В, О, /)* (24-го порядка) с 12 до 9 (что соответствует уменьшению числа двоичных разрядов в мантиссе с 40 до 30) сокращает точность определения Ху с 8—9 знаков до 4—5; матрицы с шестью верными знаками (20 двоичных разрядов в мантиссе) дают едва 1—2 верных знака в частотах.
4. Определение корней полинома. Хорошая асимптотическая сходимость многих методов: Ньютона, Мюллера, Берстоу и др. [4] недостаточна для создания надежной вычислительной процедуры определения корней полинома. Для этой цели наиболее подходящим: методом нам представляется метод спуска [6], обладающий не только асимптотической квадратичной, но и общей сходимостью.
Ниже предлагается комбинированный метод. Комбинируется метод спуска по касательной для модуля \f^z)\ с интерполяцией квадратной параболой для самой функции /(г) =^гп-\-ах г"~1 + • • • + + ап. Комбинация осуществляется на каждой итерации, причем метод спуска гарантирует общую сходимость, а интерполяция параболой ускоряет сходимость и, кроме того, выявляет комплексные корни независимо от того, начинаются итерации из вещественной или из мнимой точки.
В основе процедуры скорейшего спуска по модулю лежит формула
2*+1 А; Л = Г{г)ф0,
Пч)
которая при ¿0=1 совпадает с формулой Ньютона. В работе [6] предлагается подбирать параметр ^ последовательным дроблением
пополам ts—Ц^-^, Д° тех П0Р> пока /(гк+\) по модулю не
окажется меньше, чем /(гк). Однако, этот прием оказывается мало эффективным для кратных корней, поскольку вблизи корня 2* кратности р имеет место приближенное равенство
* * ргы
и, следовательно, для исходного значения параметра ¿„==1 хотя и
/(г \
будет реализоваться условие спуска | {* | <С 1, где [а = н0 ПОД"
ход к корню будет проходить весьма медленно: по геометрической прогрессии со знаменателем ^ ^. Поэтому мы предлагаем допол-
нить указанный прием наращиванием шага на величину исходного: ^ = ^_1 + 1 до тех пор, пока модуль функции не примет наименьшего значения. Такое наращивание проводится только для кратных корней. Признаком кратного корня является выполнение в его окрестности следующего неравенства:
, (4)
4 е
где ¡а0 — отношение значений функций в точках и гк для исходного значения параметра *0=1, т. е. в точках, полученных просто по формуле Ньютона.
Докажем неравенство (4), для чего запишем полином /(г), имеющий в точке 2* корень кратности р, в следующем виде:
/(*) = (г —2,У»«р (2).
Обозначив точку гк вблизи корня как гк = г* + й, получим
/(**)=<*’?(**);
2к+1 = гк-£^- = гк----------------------= 2* + ^ + 0 (42)-,
/(**) /> + <¥.(**)/?(**) р 1 ^ '
'р-\ 1,« Дг*+0 —1р-
/(2*4.1) = (---------- й + О (й2)] у (г*+х); Нт
/ /(^а)
Таким образом, неравенство (4) выполняется только в пределе. Чтобы не потерять случай кратных корней, достаточно несколько расширить интервал, определяемый этим неравенством, например, назначить в качестве критерия для наращивания шага выполнение неравенства
0,2<|!х0|<0,4.
Применение приема наращивания при [н-0|>0,4 не оправдано из-за отдаленности корня, а при | ¡л01 <■ 0,2 — из-за близости простого корня.
Покажем, что процедура спуска позволяет без дополнительного вычисления функции или ее производной реализовать также еще и квадратичную интерполяцию. В самом деле, при спуске из точки гк, во-первых, вычисляются две величины f(zk), /' (гк), во-вторых, как только получена новая точка 2*4-1, она немедленно проверяется на снижение, т. е. еще раз вычисляется функция /(2) теперь уже
В ТОЧКе ¿£+1.
Пусть после ряда проб (изменение параметра (3) будет получена точка г£+1, для которой выполняется условие спуска |Н<1 (верхний индекс „с“ означает, что эта точка получена процедурой спуска). Обозначим три имеющихся в нашем распоряжении величины: /(«*), /' (2Й), /(2£+1) через /, /' и /с. С помощью этих трех величин построим квадратную параболу и один из ее корней, ближайший к точке 2£+1, будем считать за новую итерацию гк+\-, в этом и состоит квадратичная интерполяция. Однако при реализации следующего шага итерационного процесса значение функции в новой точке /(2й+1) прежде всего сравнивается по модулю с /(г|+1) и если интерполяция оказалась неудачной, т. е. |/(2й+1) | >|/(2£+1) |, то в качестве новой точки остается
Л + 1
Формулу для интерполирующей параболы запишем в виде, наиболее удобном для использования результатов процедуры спуска:
/(и) = [1 ~~ 0 + и)]/+ а 4- и)]2' (5)
Корни квадратного уравнения /(и) = 0 являются теми относительными добавками к параметру ts, которые позволяют определить новое приближение для 2*4.1:
2*4-1 = 2й — ^(1+и)Л = 4+1 — и^к.
Отметим, что переменная ts — вещественная, а и — в общем комплексная. Неизвестная константа а в уравнении (5) определяется из условия /(0)=/с, т. е.
.0 — **)/ + — А-
Парабола (5) может быть преобразована к более удобному виду, поеле чего уравнение /(и) = 0 примет следующий вид:
Ы2 + {2% + ^)и- |1 = 0, где [х =/с //; 5 = 1 — ц — гз.
Из двух корней этого уравнения при
щ,, = у(-£-1- +]/е + -^)
должен быть взят корень с меньшим модулем. Случай 1 = 0 (при этом гк+1 =гк — К) должен быть отброшен как не представляющий интереса. .
5. В предлагаемом комбинированном методе определения корней полинома каждый новый итерационный шаг начинается с выбора
лучшего приближения из двух имеющихся: г\ и гк, причем для одного из них функция уже вычислена —/(г£). Если в качестве исходных значений задать 2£ = 0 и 20 = А — любое число, то такой ритм итерационного процесса будет выдержан, начиная уже с нулевого шага, поскольку значение /(0) известно без вычисления. В качестве А можно брать А = 1, А=1 или только что найденный корень. Последнее предпочтительнее, так как позволяет ускорить нахождение близких и кратных корней.
Определенный интерес представляют полиномы с вещественными коэффициентами, поскольку в этом случае возможно проводить все вычисления в вещественном поле, в том числе и для комплексных корней (схема Горнера) [7]. Однако, поскольку в этом случае найденные комплексные корни исключаются парами, т. е. опасность исключения несуществующей пары комплексных корней вместо одного простого вещественного корня, это может произойти, когда итерационный процесс, начавшись из комплексной точки,, не дошел до вещественного корня X на мнимую величину к, лежащую в пределах требуемой точности: 2* = X + ге. Формально этот случай может рассматриваться как случай комплексных корней.
Легко понять, что вопрос о вещественности корня связан с ситуацией, возникающей для пары близких или кратных корней. Известно, что, если для простых корней достигнутая относительная точность оценивается величиной а, то корень р-й кратности
р,-
может быть определен с точностью не лучшей, чем у О [8].
Для диагностики вещественности корня предлагаем следующий способ. Пусть 2* = X + — найденный корень. При £>г = 10/а
будем считать найденный корень комплексным, потому что невозможно, чтобы вещественный корень 2* = X определялся с такой плохой точностью (имеется в виду простой вещественный корень). При е<г проверим корень г*=Х-}-/е на вещественность, для чего пройдем по полуокружности
С==Х + гег<Р; 0
Если 1т/(С) будет при этом оставаться величиной постоянного знака, то внутри окружности радиуса г может быть не более одного корня. Для корня кратности р будем иметь р— 1 смен знака. Практически достаточно сделать подсчет функции /(С) в конечном числе
и т,
точек на полуокружности, в частности в пяти точках при = ;
т = 1, . . .,5. При таком малом числе точек е,сть вероятность не заметить кратности корней внутри окружности и ошибочно посчитать пару корней Х + /г за один X, но такая ошибка не портит всего хода расчета. Обратный же случай невозможен, т. е. при простом действительном! корне невозможна смена знаков у 1ш/(С) в упомянутых пяти точках и, значит, устраняется опасность исключения несуществующей пары корней вместо одного действительного.
На основе изложенного метода и с привлечением приемов преодоления затруднений в окрестности стационарных точек, где /' (2*) = 0 [6, 9], была написана стандартная программа для БЭСМ-6. В программе предусмотрен переход на двойную точность при вычислении /(2) вблизи корней (производная же вычислялась только с простой точностью). При исключении корней сохранялась двойная точность в коэффициентах полинома.
В качестве минимального шага &г, служащего критерием сходимости итерационного процесса |Л|-<ег, бралась величина вг —
= 2/?-10~12, где /? = шах (|/|аА|+аг])—теоретическая граница корней полинома [8, 10].
Кроме критерия по аргументу г, вводился также критерий по функции, т. е. г считался корнем при выполнении неравенства |/(г)|<10-13|а„|, где ап — свободный член полинома (исходного или преобразованного в результате исключения уже найденных корней).
Результаты расчетов оказались следующими. Корни полиномов порядка п = 60 вычислялись за 3 с. При этом для простых корней получалась точность 9 знаков, для кратных 5—6 знаков, для трехкратных 3—4 знака.
Характерные примеры:
2) /(г) = г60 + 24г55 + 253гъо + 1540г45 + 5985г40 + 15 504г35 +
+ 27 132г30 + 31 824г36 + 24 ЗЮг20 + 11440г15 + ЗООЗг10 + 364г5 + 13;
этот полином взят из работы [6] с заменой 2 на г5; его корни:
3) /60(г)=П 22 -)- - ; 60 корней этого полинома распола-
*=л ^
гаются на мнимой оси (задача, характерная для флаттера);
10, в отличие от предыдущего примера, в котором все 60 корней вычислены с 9 знаками, были определены всего лишь с 7 знаками — это пример плохой обусловленности полинома [4].
Преодоление затруднений около стационарных точек проверялось на следующем полиноме:
/, (г) = г8 - 8г7 + (24 + г) - (32 + 6е)г5 + (14 + Не) г4 + (8 - 4е) г8 -- (24 4 5в) г2 + (32 + 2е) г - 15(1 + е);
его корни: 1 + г тА ; —1; 4 /; 2+г; 3.
Около стационарной точки г= 1 этого полинома располагается пара близких комплексных корней. Задавалось два варианта: е = = 10~8 и е=10-10; в первом варианте получилось семь верных знаков: = 1 + 0,0001000/, во втором — шесть знаков = 1 +0,000010/.
В полиноме /40 (г) = (г16 — было задано е=10-4, так что
среди 40 его корней три оказались трехкратными: —1, +/, — / и четверка близких: 1, 1, 1 + 0,01/. В этом примере трехкратные корни вычислялись с четырьмя знаками, а четверка близких корней получилась такой: 1,0002, 0,9998, 1,000000+0,010002/ (из двух расчетов с исходной итерацией 20=1 и 20 ==/ оставлен менее точный результат).
1)
/(г) = (г3® - 1)*; ./(г) = (*18-1)8;
30
20 / £, Ч
4) /20(2)= П (г-у)
1. Буньков В. Г. Особенности свободной схемы летательного аппарата при решении задач аэроупругости. Труды ЦАГИ, вып. 1166, 1969.
2. Буньков В. Г. Учет деформации сдвига при расчете колебаний крыла малого удлинения методом многочленов. .Ученые записки ЦАГИ“, т. III, № 4, 1972.
3. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. М., Физматгиз, 1963.
4. Уилкинсон Дж. X. Алгебраическая проблема собственных значений. М., „Наука“, 1968.
5. G i v е п s W. The characteristic value-vector problem. J. Assoc. Coraput. Machinery, 1957, 4, N 3.
6. В о e в о д и н В. В. Численные методы алгебры. М. .Наука“, 1966.
7. Березин И. С., Жидков Н. П. Методы вычислений. Т. II, М., Физматгиз, 1960.
8. Островский А. М. Решение уравнений и систем уравнений. М., Изд. иностр. лит., 1963.
9. К а ц И. С., Маергойз М. Д. Решение нелинейных алгебраических и трансцендентных уравнений в комплексной области. „Журн. вычисл. матем. и матем. физики“, т. 7, №3, 1967.
10. Lagrange I. L. Oeuves completes, т. VIII, Париж, 1867-92. ■
Рукопись поступила 26/IV ¡974
