close

Вход

Забыли?

вход по аккаунту

?

26.Применение матричных методов для расчета частот и форм свободных колебаний динамических моделей силовых передач колесных машин с конечным числом степеней свободы

код для вставкиСкачать
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Московский государственный технический университет
имени Н.Э. Баумана
М.Г. Лахтюхов
ПРИМЕНЕНИЕ МАТРИЧНЫХ МЕТОДОВ
ДЛЯ РАСЧЕТА ЧАСТОТ И ФОРМ
СВОБОДНЫХ КОЛЕБАНИЙ ДИНАМИЧЕСКИХ
МОДЕЛЕЙ СИЛОВЫХ ПЕРЕДАЧ
КОЛЕСНЫХ МАШИН С КОНЕЧНЫМ
ЧИСЛОМ СТЕПЕНЕЙ СВОБОДЫ
Рекомендовано редсоветом МГТУ им. Н.Э. Баумана в качестве
учебного пособия по курсу «Динамика колесных машин»
Под редакцией А.А. Полунгяна
Москва
Издательство МГТУ им. Н.Э. Баумана
2007
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
УДК 629.113(075.8)
ББК 39.33-04+22.213
Л29
Рецензенты: Е.А. Галевский, С.В. Аринчев
Лахтюхов М.Г.
Л29
Применение матричных методов для расчета частот и форм
свободных колебаний динамических моделей силовых передач
колесных машин с конечным числом степеней свободы: Учеб.
пособие по курсу «Динамика колесных машин» / Под ред.
А.А. Полунгяна. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2007. –
60 с.: ил.
ISBN 5-7038-2933-Х
Даны рекомендации по выбору матричных методов расчета частот
и форм свободных колебаний консервативных и неконсервативных
динамических моделей силовых передач колесных машин с конечным
числом степеней свободы в зависимости от типа матриц и их
размерности, объема решаемой задачи (расчета всех или части частот и,
возможно, форм свободных колебаний) и т. д. Показано, что решение
задачи нахождения частот и форм свободных колебаний математически
эквивалентно решению задачи на собственные значения. Рассмотрены
особенности расчета собственных значений и собственных векторов.
Приведены сведения о программных средствах для решения спектральных задач. Включены примеры расчета частот и форм свободных
колебаний динамических систем в среде MathCAD.
Для студентов старших курсов специальностей «Автомобилестроение» и «Многоцелевые гусеничные и колесные машины».
УДК 629.113(075.8)
ББК 39.33-04+22.213
Учебное издание
Михаил Георгиевич Лахтюхов
Применение матричных методов
для расчета частот и форм свободных колебаний
динамических моделей силовых передач
колесных машин с конечным числом степеней свободы
Редактор С.А. Серебрякова
Корректор Л.И. Малютина
Компьютерная верстка А.Ю. Ураловой
Подписано в печать 26.12.2006. Формат 60×84/16. Бумага офсетная.
Печ. л. 3,75. Усл. печ. л. 3,49. Уч.-изд. л. 3,25.
Тираж 100 экз. Изд. № 104. Заказ
Издательство МГТУ им. Н.Э. Баумана.
105005, Москва, 2-я Бауманская ул., 5.
ISBN 5-7038-2933-Х
 МГТУ им. Н.Э. Баумана, 2006
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Поведение силовой передачи колесной машины (КМ) при динамическом воздействии во многом определяется ее частотами и
формами свободных колебаний. При проведении исследований
расчет частот и форм собственных колебаний динамических систем является одной из важнейших задач. При этом можно выделить несколько задач, отличающихся по объему вычислений.
В зависимости от поставленных целей может возникнуть необходимость вычисления: нескольких низших собственных частот и,
возможно, соответствующих форм колебаний; всех собственных
частот и соответствующих форм колебаний; только собственных
частот (всех или нескольких); числа собственных частот, лежащих
в заданном диапазоне, и т. п.
1. УРАВНЕНИЕ СВОБОДНЫХ КОЛЕБАНИЙ
ДИНАМИЧЕСКОЙ СИСТЕМЫ
С КОНЕЧНЫМ ЧИСЛОМ СТЕПЕНЕЙ СВОБОДЫ
Теоретическому исследованию динамической нагруженности силовой передачи КМ предшествует ее схематизация. При исследованиях в низкочастотном диапазоне силовая передача КМ (рис. 1),
представляющая собой систему с рассредоточенными параметрами,
как правило, моделируется линейными системами с сосредоточенными параметрами (рис. 2), состоящими из конечного числа жестких
масс, совершающих вращательное движение и соединенных между
собой безынерционными упругими валами. Число степеней свободы
при этом не превосходит 150.
Дифференциальные уравнения голономной динамической системы с сосредоточенными параметрами и стационарными связями
могут быть получены с помощью дифференциальных уравнений
Лагранжа второго рода
d ∂T ∂T ∂Ф ∂П
−
+
+
= Qi ,
dt ∂qi ∂qi ∂qi ∂qi
i = 1, n,
где T , П – соответственно кинетическая и потенциальная энергии
системы; Ф – диссипативная функция системы; Qi – i-я обобщенная
3
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сила; qi , qi – соответственно i-е обобщенные координата и скорость;
t – время; n – число степеней свободы динамической системы.
При малых колебаниях около положения равновесия кинетическая и потенциальная энергии системы, а также ее диссипативная
функция могут быть выражены квадратичными формами от обобщенных координат:
T=
1 n n
1 n n
1 n n
aij qi q j , П = ∑∑ сij qi q j , Ф = ∑∑ bij qi q j ,
∑∑
2 i =1 j =1
2 i =1 j =1
2 i =1 j =1
(1)
где aij , cij – соответственно инерционные и жесткостные коэффициенты; bij – коэффициенты диссипации.
Рис. 1. Кинематическая схема силовой передачи колесной машины (8×8)
с раздачей крутящего момента по бортам
4
момента по бортам
Рис.2. Приведенная динамическая система силовой передачи колесной машины (8×8) с раздачей крутящего
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
5
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Уравнение, описывающее движение линейной динамической
модели силовой передачи с n степенями свободы, может быть записано в матричной форме
Aq + Bq + Cq = F ,
(2)
где A, C – матрицы инерционных и жесткостных коэффициентов
(размерности n × n ); B – матрица коэффициентов диссипации
(демпфирования) ( n × n ); q – n- мерный вектор обобщенных координат; F – n- мерный вектор внешних сил.
Если вектор внешних сил F равен нулевому вектору 0, то из
уравнения (2) получаем матричное уравнение
Aq + Bq + Cq = 0,
(3)
или, что одно и то же, систему однородных дифференциальных
уравнений,
a11q1 + a12 q2 + ... + a1n qn + b11q1 + b12 q2 + ... + b1n qn +

+ c11q1 + c12 q2 + ... + c1n qn = 0,


a21q1 + a22 q2 + ... + a2 n qn + b21q1 + b22 q2 + ... + b2 n qn +

+ c21q1 + c22 q2 + ... + c2 n qn = 0,


. . . . . . . . . . . . . . . . . . . . . . .
a q + a q + ... + a q + b q + b q + ... + b q +
nn n
n1 1
n2 2
nn n
 n1 1 n 2 2

+ cn1q1 + cn 2 q2 + ... + cnn qn = 0.

(4)
Уравнение (3) и система уравнений (4) описывают малые свободные колебания линейной неконсервативной системы.
Выражения для квадратичных форм (1) и соответственно
структура матриц A , B и C определяются выбором обобщенных координат q. На практике в качестве обобщенных координат q обычно принимают отклонения масс от положения равновесия. В этом случае выражения (1) для кинетической и потенциальной энергий и для диссипативной функции преобразуются к виду
6
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
T=
1 n
1 n −1 n
J i qi2 , П = ∑ ∑ kij qi − q j
∑
2 i =1
2 i =1 j =i +1
(
Ф=
1 n −1 n
∑ ∑ fij qi − q j
2 i =1 j =i +1
(
)
2
+
)
2
+
1 n
∑ ki qi2 ,
2 i =1
1 n
∑ fi qi2 ,
2 i =1
(5)
где J i – момент инерции i-й массы; kij , fij − соответственно коэффициенты жесткости и диссипации упругого участка, соединяющего i-ю и j-ю массы; ki , fi − соответственно коэффициенты
жесткости и диссипации упругого участка, соединяющего i-ю массу с заделкой.
Все три квадратичные формы в (5) положительны.
Уравнение движения произвольной i-й массы (рис. 3), связанной с j-й, l-й, m-й и т. д. массами и заделкой, имеет вид
(
)
J i qi + fij + fil + fim + … + fi qi − fij q j − fim qm − fil ql − … +
(
)
+ kij + kil + kim + … + ki qi
− kij q j − kim qm − kil ql − … = 0.
(6)
Рис. 3. Произвольная i-я масса динамической системы
Нетрудно убедиться, что матрица инерционных коэффициентов A в данном случае является диагональной и, следовательно,
симметричной. На ее главной диагонали располагаются моменты
7
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
инерции масс. Все внедиагональные элементы матрицы равны нулю. Так как все моменты инерции системы и, следовательно, все
диагональные элементы матрицы A положительны:
aii = J i > 0, i = 1, n,
то матрица A, кроме того, является положительно определенной.
Матрица жесткостных коэффициентов C также является симметричной. Ее внедиагональный элемент cij , расположенный на
пересечении i-й строки и j-го столбца ( i ≠ j ) , равен взятому со
знаком минус коэффициенту жесткости kij упругого участка, соединяющего i-ю и j-ю массы:
cij = − kij ,
(7)
если связь между массами есть, или равна нулю:
cij = 0,
если этой связи нет.
Элементы матрицы C , расположенные на главной диагонали,
равны сумме коэффициентов жесткостей, характеризующих упругие связи i-й массы с заделкой и другими массами:
cii = ki +
n
∑
j =1 (j ≠ i )
kij , i = 1, n.
(8)
Выражения для потенциальной энергии и диссипативной функции в
(5), как и коэффициенты при q и q в (6), одинаковы по структуре. Поэтому матрицы B и C также имеют сходную структуру, и, следовательно, матрица B также является симметричной. Формирование матриц B и
C подчиняется одним и тем же правилам.
На практике заполнение матриц B и C удобно проводить одновременно. Вначале все элементы обеих матриц приравниваются
к нулю. А затем в соответствии со связями в динамической системе коэффициенты жесткости kij согласно (7) и (8) прибавляются к
элементам cij , c ji , cii и c jj матрицы C , а коэффициенты fij с
8
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
такими же знаками, что и kij , − к элементам bij , b ji , bii и b jj матрицы B. Коэффициенты ki и fi , характеризующие упругий участок, который соединяет i-ю массу с заделкой, прибавляют соответственно к диагональным элементам матриц C и B, расположенным на пересечении i-й строки и i-го столбца.
2. ЭКВИВАЛЕНТНОСТЬ НАХОЖДЕНИЯ СОБСТВЕННЫХ
ЧАСТОТ И РЕШЕНИЯ ЗАДАЧИ НА СОБСТВЕННЫЕ
ЗНАЧЕНИЯ
2.1. Линейные консервативные системы
При отсутствии демпфирования в системе все коэффициенты
диссипации bij , являющиеся элементами матрицы B в уравнении (3), равны нулю:
bij = 0, i = 1, n; j = 1, n.
Тогда после упрощения уравнения (3) и системы уравнений (4)
приходим к матричному уравнению
Aq + Cq = 0
(9)
и системе однородных дифференциальных уравнений
a11q1 + a12 q2 + ... + a1n qn + c11q1 + c12 q2 + ... + c1n qn = 0,

a21q1 + a22 q2 + ... + a2 n qn + c21q1 + c22 q2 + ... + c2 n qn = 0,

. . . . . . . . . . . . . . . . . . . . . . . . .

an1q1 + an 2 q2 + ... + ann qn + cn1q1 + cn 2 q2 + ... + cnn qn = 0,
(10)
которые описывают малые свободные колебания линейной консервативной системы.
Частное решение однородной системы линейных дифференциальных уравнений (10) соответствует свободным колебаниям системы с s-й круговой частотой и имеет вид
qs = qsМ sin(ωs t + ϕs ),
(11)
9
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 q1Мs

 qМ
=  2s
 ...
 qМ
 ns



М
где qsМ
 − искомый амплитудный n-мерный вектор; q js –



искомая амплитуда колебаний j-й массы ( j = 1, n) при свободных
колебаниях с s-й круговой частотой; ωs , ϕs − искомые s-я круговая частота и соответствующий фазовый угол.
Колебания, описываемые уравнением (11), называются главными.
Общим решением системы линейных дифференциальных уравнений (10) является линейная комбинация частных решений (11)
n
q = ∑ qsМ sin(ωs t + ϕs ).
(12)
s =1
При подстановке частного решения (11) в (10) получим
(
(
)
)
(
(
)
)
(
(
) q = 0,
) q = 0,
 c11 − a11ω2s q1Мs + c12 − a12ωs2 q2Мs + ... + c1n − a1n ω2s

 c − a ω2 q М + c − a ω2 q М + ... + c − a ω2
 21 21 s 1s
22
22 s
2s
2n
2n s

. . . . . . . . . . . . . . . . . . . . . . .

2
М
2
М
2
 cn1 − an1ωs q1s + cn 2 − an 2 ωs q2 s + ... + cnn − ann ωs
(
)
(
)
(
М
ns
М
ns
. .
)q
М
ns
=0
или в матричной форме
(C − ω A) q
2
s
М
s
= 0.
(13)
Уравнение (13) можно преобразовать к виду
(A
−1
)
C − ω2s I qsМ = 0 ,
(14)
где A−1 − матрица, обратная матрице A ; I − единичная матрица
(размерности n × n ).
10
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Известно [7], что система (14) имеет ненулевое решение qsМ
только тогда, когда определитель матрицы
(A
−1
C − ω2s I
)
равен
нулю:
A−1C − ω2s I = 0.
(15)
Уравнение (15), называемое характеристическим, является алгебраическим уравнением n-й степени относительно ω2s и имеет в
общем случае n различных корней λ s = ω2s (s = 1, n ).
Матричное уравнение вида (13) представляет собой обобщенную форму записи алгебраической проблемы собственных значений. Корень характеристического уравнения (15) λ s = ω2s и ненулевой вектор qsМ , удовлетворяющие уравнению (13), в линейной
алгебре называются соответственно собственным значением и
собственным вектором. Совокупности собственных значений и
собственных векторов образуют соответственно спектры собственных значений и собственных векторов.
Из симметричности вещественных матриц A и C и из положительной определенности матрицы A следует, что все корни характеристического уравнения (15) или, что то же самое, все собственные значения λ s (s = 1, n ) являются действительными и неотрицательными [5]. Также действительными будут и все собственные векторы системы qsМ (s = 1, n ) .
При подстановке в систему (13) круговой частоты ωs , соответствующей произвольному s-му частному решению, любое из уравнений этой системы оказывается следствием остальных, а сама
система уравнений становится неопределенной, поскольку на n
независимых переменных приходится ( n − 1) независимых линейных уравнений.
Для устранения указанной неопределенности системы линейных уравнений на практике применяют два приема, основанные на
вычислении не самих амплитуд колебаний масс q М
js , а их относительных величин x js .
11
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В первом случае все амплитуды колебаний масс системы выражаются через одну из них, имеющую ненулевое значение, например, через q1Мs :
qsМ = q1Мs xs ,
(16)
где q1Мs − амплитуда колебаний первой массы с s-й круговой частотой, причем q1Мs ≠ 0; xs − вектор, содержащий относительные
амплитуды колебаний масс системы с s-й круговой частотой, где, в
М
(j = 2, n).
свою очередь, x1s = 1, x js = q М
js q1s
Подставив уравнение (16) в (13) и разделив обе части уравнения на ненулевой множитель q1Мs , получим систему линейных однородных алгебраических уравнений
( C − ω A) x
2
s
s
= 0,
которая является определенной, так как на
(17)
( n − 1)
неизвестную
приходится уже ( n − 1) независимых уравнений.
Во втором случае вектор, содержащий относительные амплитуды колебаний, нормируют:
xНт s xН s = 1,
(18)
где xНт s − транспонированный вектор xН s .
При этом векторы qsМ и xН s связаны между собой равенством
qsМ = qsМ xН s ,
где qsМ − длина вектора qsМ , qsМ =
(19)
(q ) + (q )
М 2
1s
М 2
2s
( ).
М
+ … + qns
2
После подстановки (19) в (13) и деления обеих частей уравнения на ненулевую величину qsМ , получаем
12
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
(C − ω A) x
2
s
Нs
= 0.
(20)
Теперь на систему линейных однородных уравнений (20) вместе с уравнением (18) на n неизвестных приходится n независимых уравнений.
Все три вектора qsМ , xs и xН s , различаясь по длине, задают
одно и то же направление, что следует из (16) и (19), и, следовательно, являются собственными векторами, соответствующими
s-му собственному значению λ s = ω2s .
Отличительной особенностью векторов xs и xН s является то,
что их координаты x js и xН js
( j = 1, n ) ,
представляющие собой
относительные амплитуды перемещений масс при свободных колебаниях, однозначно определяются свойствами системы. Абсолютные амплитуды колебаний q М
js
( j = 1, n ) зависят не только от
свойств системы, но и от масштабных множителей q1Мs в (16) и
qsМ в (19), которые пока остаются неизвестными.
Совокупность относительных амплитуд x js или xН js
( j = 1, n ) ,
характеризующая свободные колебания масс динамической системы с s-й круговой частотой ωs и зависящая только от свойств этой
системы, в теории колебаний называется s-й собственной формой,
а сами относительные амплитуды x js и xН js − коэффициентами
s-й собственной формы.
Наглядное представление о характере собственных колебаний
динамической системы дает графическое изображение собственных форм, которое строится следующим образом. На расчетной
схеме динамической системы (рис. 4) из центра каждой массы в
удобном для исследователя масштабе перпендикулярно отрезкам,
изображающим безынерционные упругие валы, строят отрезки,
длины которых пропорциональны относительным амплитудам колебаний соответствующих масс. Причем отрезок откладывают над
валом при положительном значении относительной амплитуды и
под валом − при отрицательном значении. Соединив друг с другом
концы полученных подобным образом отрезков, получают график,
13
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
изображающий формы колебаний. Точки пересечения графика с
безынерционными упругими валами соответствуют неподвижным
точкам этих валов и называются узлами.
Рис. 4. Формы колебаний динамической модели с 11 степенями свободы:
а – первая форма колебаний; б – вторая форма колебаний
Если на расчетной схеме динамической системы массы расположить друг от друга на расстоянии, обратно пропорциональном
коэффициентам жесткости соответствующих участков, соединяющих эти массы, то с помощью графика собственной формы может
быть выявлен участок системы, на котором при свободных колебаниях действует наибольший упругий момент. Наибольший момент будет действовать на участке, где график собственной формы
образует наибольший угол к отрезку, изображающему упругий
(
)
M
M
вал. В самом деле, величина qisM − q M
js kij , где qis , q js − макси-
мальные угловые отклонения i-й и j-й масс от положения равновесия при свободных колебаниях с круговой частотой ωs , с одной
стороны, равна упругому моменту, действующему на упругом
участке, соединяющем i-ю и j-ю массы, а с другой стороны, пропорциональна тангенсу угла наклона графика собственной формы
14
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
к отрезку, изображающему упругий вал. Поэтому больший угол
наклона графика будет соответствовать большему упругому моменту на участке, соединяющем массы.
Набор круговых частот ωs s = 1, n называется собственным
(
)
или, что то же самое, частотным спектром системы. Частоты и
формы свободных колебаний определяются исключительно параметрами динамической системы и ее структурой и не зависят от
внешних условий.
Для удобства последующих рассуждений независимо от применяемого метода, устраняющего неопределенность системы (13),
обозначим неизвестные масштабные множители q1Мs и qsМ через
qМs , а векторы xs и xНs − через x s . Тогда вектор qsМ будет определяться как
qsМ = qМs xs .
(21)
С учетом (21) частное решение системы уравнений (14), соответствующее s-й круговой частоте, можно записать в виде
qs = qМs xs sin(ωs t + ϕs ).
В свою очередь, общее решение (12) системы дифференциальных уравнений (13) примет вид
n
q = ∑ qМs xs sin(ωs t + ϕs )
(22)
s =1
или для произвольной j-й массы
n
q j = ∑ qМs x js sin(ωs t + ϕs ),
j = 1, n.
(23)
s =1
Если частоты и формы колебаний найдены в результате решения характеристического уравнения (15) или задачи на собственные значения (17) или (20), то 2n неизвестных постоянных
15
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
qМs и ϕs ( s = 1, n), входящих в (22), могут быть определены из
начальных условий − перемещений и скоростей всех масс в момент времени t = 0.
2.2. Линейные неконсервативные системы
Для линейной неконсервативной системы малые свободные
колебания описываются матричным уравнением (3). Частное решение уравнения (3) следует искать в виде
q (t ) = q М e λ t ,
(24)
где λ − искомое число; qМ − искомый амплитудный вектор.
Подстановка (24) в (3) приводит к однородной системе алгебраических уравнений относительно координат вектора qМ
( λ A + λB + C ) q
2
М
= 0,
(25)
которая имеет ненулевые решения qМ только тогда, когда опреде-
(
)
литель матрицы λ 2 A + λB + C равен нулю:
λ 2 A + λB + C = 0.
(26)
Полученное уравнение (26), называемое характеристическим,
является алгебраическим уравнением степени 2n относительно λ
и имеет 2n корней.
Поскольку система линейных уравнений (25) однородна, то
(
)
при подстановке в нее найденных корней λ s s = 1, 2n любое из
уравнений является линейной комбинацией всех остальных, и,
следовательно, компоненты искомого амплитудного вектора qsМ
могут быть найдены с точностью до постоянного множителя. Поэтому вектор qsМ , как и ранее, выразим через вектор собственной
формы колебаний xs
qsМ = qМs xs ,
16
(27)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1 


 x2 s 
где x s =  x3s  − искомый вектор собственной формы колебаний,


 …
 x 
 ns 
М
где, в свою очередь, x js = q М
js q1s (j = 1, n ); qМs − неопределен-
ный множитель, равный амплитуде перемещений первой массы
q1Мs . Предполагается, что q1Мs ≠ 0.
Тогда частное решение уравнения (3) примет вид
qs (t ) = qМs xs eλ st ,
(28)
и система линейных однородных алгебраических уравнений (25)
может быть заменена системой
( λ A + λB + C ) x = 0,
2
(29)
в которой на ( n − 1) неизвестное приходится ( n − 1) независимых
уравнений.
Общее решение уравнения (3) можно записать в виде
2n
q (t ) = ∑ qМs x s eλ st .
(30)
s =1
Для положительных квадратичных форм (5) с действительными коэффициентами корни характеристического уравнения (26)
могут быть действительными или комплексно-сопряженными.
Действительные корни всегда отрицательны или равны нулю.
Также отрицательны или равны нулю действительные части комплексных корней [8]. Причем действительные части всех комплексных чисел одновременно обращаются в нуль только тогда,
когда демпфирование в системе отсутствует, т. е. все элементы
матрицы B равны нулю ( bij = 0, i = 1, n, j = 1, n ). В этом случае
уравнение (26) может быть преобразовано в (15), а потому уравнение
17
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
(3), описывающее малые свободные колебания консервативной системы с нулевой матрицей B, следует заменить более простым уравнением (9).
Действительным корням λ s = −δ s ( δ s ≥ 0 ) характеристического уравнения (26) соответствуют действительные векторы x s ,
комплексно-сопряженным корням λ s = −δ s + iωs и λ*s = −δ s − iωs
( δs ≥ 0 )
− комплексно-сопряженные векторы собственных форм
колебаний
x s = Re xs + i Im xs = us + iv s ,
x *s = Re x s − i Im xs = us − iv s .
(31)
Каждому отрицательному действительному корню соответствуют апериодическое затухающее движение, описываемое уравнением вида
qs (t ) = qМs xs e −δst ( δs > 0 ) ,
(32)
а нулевому действительному корню − вращение всей системы как
единого целого.
При комплексных корнях система совершает затухающие колебания, описываемые уравнениями вида
qs (t ) = qМs xs e(
−δ s + iωs )t
.
(33)
Если имеются действительные и комплексные корни, то затухающие колебания накладываются на затухающее движение.
Неотрицательное действительное число δs в (32) и (33) характеризует демпфирование системы и называется коэффициентом
демпфирования, а положительное действительное число ωs в (33)
представляет собой s-ю круговую частоту колебаний системы с
демпфированием (3).
Чрезвычайно важным для практики является то, что, несмотря
на комплексный вид корней и форм колебаний, движение системы
может быть описано действительной функцией.
18
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для пояснения последующих рассуждений вспомним некоторые сведения из теории комплексных чисел.
На плоскости с прямоугольной системой координат произвольное комплексное число z = a + ib изображается точкой Z
(рис. 5). Абсцисса этой точки равна действительной части комплексного числа Re z = a , а ордината − мнимой части Im z = b.
Рис. 5. Изображение комплексного числа на плоскости
Сопряженное комплексное число z* = a − ib на числовой плоскости изображается точкой Z * , получающейся при зеркальном
отображении точки Z относительно горизонтальной оси.
Также комплексное число z может быть изображено вектором
OZ , выходящим из начала координат к точке Z . Вектор OZ (для
z ≠ 0 ) характеризуется длиной и углом наклона вектора к оси абсцисс. Длина вектора называется модулем комплексного числа.
Модули комплексно-сопряженных чисел равны
z = z * = a 2 + b 2 = r.
Действительная Re z и мнимая Im z части комплексного числа
z ≠ 0 могут быть выражены через его модуль z и угол ϕ между
осью абсцисс и вектором, изображающим это число,
19
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
a = r cos ϕ = z cos ϕ,
b = r sin ϕ = z sin ϕ.
Представление комплексного числа z в виде
z = z (cos ϕ + i sin ϕ)
(34)
называется его тригонометрической формой.
С учетом формулы Эйлера
cos ϕ + i sin ϕ = eiϕ
возможен переход от тригонометрической формы записи комплексного числа к экспоненциальной
z = z e iϕ
(35)
и наоборот.
Предположим, что неконсервативная система имеет только
комплексно-сопряженные корни уравнения (26), причем все корни
различны. Каждой паре комплексно-сопряженных корней
λ s = −δ s + iωs и λ*s = −δ s − iωs уравнения (26) соответствуют
комплексно-сопряженные формы колебаний (31) и комплексносопряженные неопределенные комплексные коэффициенты
qМs = Re qМs + i Im qМs = qМs eiα s
и
q*Мs = Re qМs − i Im qМs = qМs e−iα s .
Рассмотрим далее сумму частных решений уравнения (3), соответствующую комплексно-сопряженным корням λ s и λ*s ,
*
qs = qМs xs eλ st + q*Мs xs*eλ st ,
20
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
которую можно преобразовать к виду
q( s ) = qМs eiα s (us +iv s )e(
−δ s +iωs )t
= qМs e−δst (us +ivs )e (

i ωs t+α s )
+ qМs e−iα s (us −ivs )e(
−δ s −iωs )t
=
−i ω t+α
+(us −ivs )e ( s s )  =

= qМs e−δst (us +ivs )(cos(ωs t + α s )+i sin(ωs t + α s ))+

+(us −iv s )(cos(ωs t + α s )−i sin(ωs t + α s )) =

= qМs e−δst us cos(ωs t + α s )+iv s cos(ωs t + α s )+
+ius sin(ωs t + α s )−v s sin(ωs t + α s )+ us cos(ωs t + α s )−
−ivs cos(ωs t + α s )−ius sin(ωs t + α s )−v s sin(ωs t + α s ) =
= 2 qМs e−δst us cos(ωs t + α s )−v s sin(ωs t + α s ) =
= 2 qМs e−δst Re xs cos(ωs t + α s )− Im xs sin(ωs t + α s ) .
С учетом полученного результата комплексная векторная
форма (30), описывающая свободные колебания линейной неконсервативной системы, может быть заменена действительной векторной формой
n
q (t ) = 2∑ qМs e−δst  Re xs cos (ωs t + α s ) − Im xs sin (ωs t + α s ). (36)
s =1
Полученное уравнение (36) является общим решением уравнения (3) и содержит 2n неизвестных постоянных qМs и α s
21
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
(s = 1, n),
которые могут быть определены из начальных усло-
вий − перемещений и скоростей масс в момент времени t = 0.
Учитывая, что любое комплексное число z ≠ 0 и, в частности,
произвольная r-я комплексная координата вектора xs могут быть
представлены в экспоненциальной форме (35), общее решение
уравнения (3), описывающее движение произвольной r-й массы
системы, также может быть записано в скалярной действительной
форме. В самом деле,
n
qr (t ) = 2∑ qМs e−δst Re xrs cos(ωs t +α s )−Im xrs sin(ωs t +α s ) =
s=1
n
(
)
(
)
= 2∑ qМs e−δst Re xrs eiϕrs cos(ωs t +α s )−Im xrs eiϕrs sin(ωs t +α s ) =


s=1
n
= 2∑ qМs e−δst Re( xrs ⋅(cosϕrs +isin ϕrs ))cos(ωs t +α s )−

s=1
−Im( xrs (cosϕrs +isin ϕrs ))sin(ωs t +α s ) =

n
= 2∑ qМs e−δst  xrs ⋅cosϕrs cos(ωs t +α s )− xrs sin ϕrs sin(ωs t +α s )
s=1
или
n
qr (t ) = 2∑ qМs e−δst xrs cos(ωs t + α s + ϕrs ),
(37)
s=1
где xrs − модуль r-й комплексной координаты вектора x s ; ϕrs −
относительная фаза свободных колебаний r-й массы с s-й круговой
частотой.
Принципиальное отличие комплексной формы колебаний от
действительной, следующее из сравнения (23) и (37), заключается
22
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
в том, что массы неконсервативной динамической системы не одновременно достигают крайних положений и не одновременно
проходят положение статического равновесия, т. е. перемещения
масс происходят с относительными фазовыми сдвигами ϕrs
(r =1, n) [1]. По этой причине при свободных колебаниях точки
смены знака упругих смещений изменяют свое положение относительно масс динамической системы, и понятие узла формы колебаний теряет смысл.
Искомые параметры λ s и xs неконсервативной системы (3),
определяемые свойствами системы и не зависящие от внешних
условий, также могут быть найдены в результате решения задачи
на собственные значения, для чего система n дифференциальных
уравнений второго порядка должна быть приведена к системе 2n
дифференциальных уравнений первого порядка. С этой целью
расширим вектор неизвестных, добавив к обобщенным координатам q еще и вектор обобщенных скоростей v = q, который с учетом (24) равен
v = λq.
(38)
Тогда однородная система дифференциальных уравнений (3)
может быть заменена системой
q −v = 0,

Av + Bv + Cq = 0.
(39)
Умножение нижнего матричного уравнения в (39) на матрицу
A , являющуюся обратной к матрице A, и добавление к соответствующим множителям обоих уравнений единичных I и нулевых
O матриц (размерностью n×n ) позволяет перейти к системе
уравнений
−1
 Iq + Ov +
Oq −
Iv = 0,



−1
−1


 Oq + Iv + A Cq + A Bv = 0
или
23
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 I


O
O
I
 O

 Q+
 −1

 A C

−I 
 Q = 0,
A−1B
где Q; Q − 2n-мерный расширенный вектор неизвестных и его
q
q
производная по времени, причем Q = , Q = .
v
v
С учетом (24) и (38) вектор Q может быть выражен через вектор Q:
q q λq 
q 
q
Q = = = 2 = λ⋅ = λ⋅ = λQ.
v q λ q
λq
v
Отcюда получаем уравнение

O

− A−1C


 Q = λQ
− A−1 B
I
(40)
или
DQ = λQ ,
(41)
которое с учетом (27) можно записать в виде
DX = λX ,

O
где D =  −1
− A C
(42)

 − действительная матрица общего вида
− A−1B
x 
размерностью 2n×2n, X =  − 2n-мерный вектор.
λx
I
Уравнение (42) представляет собой стандартную форму записи алгебраической проблемы собственных значений. Число λ
есть собственное значение матрицы D, вектор X − собственный вектор. Собственный вектор X матрицы D состоит из
24
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
двух векторов x и λx , первый из которых является искомым
вектором собственной формы, а второй получается из первого
умножением на λ .
Таким образом, и для консервативной, и для неконсервативной
динамических систем задача нахождения частот и форм свободных колебаний математически эквивалентна решению задачи на
собственные значения.
3. ОСОБЕННОСТИ РАСЧЕТА СОБСТВЕННЫХ ЗНАЧЕНИЙ
И СОБСТВЕННЫХ ВЕКТОРОВ
Итак, нахождение собственных частот ω и соответствующих
им форм x линейной консервативной системы, малые свободные
колебания которой описываются уравнением (9), эквивалентно
решению обобщенной алгебраической проблемы собственных
значений линейной системы
(C −λA) x = 0
(43)
или, что то же самое, стандартной алгебраической проблемы собственных значений
Dx = λx ,
(44)
где D = A−1C − действительная несимметричная матрица, где, в
свою очередь, A−1 − матрица, обратная матрице A; λ, x − действительные собственное значение ( λ = ω2 ) и соответствующий собственный вектор.
Разложение матрицы A в виде
A= LLт ,
(45)
где L − нижняя треугольная матрица; Lт − транспонированная
матрица L, позволяет уравнение (43) привести к стандартному
виду с симметричной матрицей
Fy = λ y ,
(46)
25
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
т
( )
где F = L−1C L−1
− действительная симметричная матрица; y −
собственный вектор.
Матрицы D и F в (44) и (46) имеют одинаковые собственные значения, но разные собственные векторы, связанные соотношением
y = Lт x.
(47)
Корни λ характеристического уравнения (26) и векторы собственных форм x линейной неконсервативной системы (3) находят в результате решения стандартной алгебраической проблемы
собственных значений
DX = λX ,
 O
где D =  −1
− A C
(48)

 − действительная матрица общего вида
− A B
x 
размерностью 2n×2n; X =  − 2n-мерный собственный вектор.
λx
Собственные значения λ матрицы D (48) могут быть действительными и комплексно-сопряженными. Действительным собственным значениям соответствуют действительные собственные
векторы, комплексно-сопряженным собственным значениям –
комплексно-сопряженные собственные векторы.
Можно показать, что расчет корней λ характеристического
уравнения (26) и векторов собственных форм x линейной неконсервативной системы (3) также приводится к решению обобщенной алгебраической проблемы собственных значений
I
−1
(P −λS ) X = 0,
 O
где P = 
 −C
A 
 − действительная матрица общего вида разO 
 A
O
мерностью 2n×2n; S = 
 B
A
го вида размерностью 2n×2n.
26
(49)

 − действительная матрица обще
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Собственные значения λ и собственные векторы X , вычисленные в результате решения задачи (48), равны соответствующим
собственным значениям и собственным векторам в задаче (49).
Собственные значения произвольной матрицы T (n×n) являются корнями характеристического уравнения
T −λI = 0.
(50)
Левая часть характеристического уравнения (50) представляет
собой многочлен от λ степени n:
p(λ ) = λ n +t1λ n−1 +…+tn−1λ +tn .
(51)
Корни многочлена (51) являются аналитическими функциями
его коэффициентов. Вычисление корней многочлена (51) для матриц небольшого порядка ( n<5 ) не представляет сложности. Однако для матриц с n≥5 уравнение (50) в общем виде оказывается
неразрешимым в радикалах, т. е. корни уравнения не могут быть
выражены через его коэффициенты при помощи некоторой комбинации радикалов [7]. Поэтому для нахождения собственных
значений с большим числом степеней свободы следует применять
численные методы линейной алгебры.
Поскольку, с одной стороны, собственные значения матрицы
T являются корнями характеристического многочлена, а с другой
стороны, по любому многочлену p(λ ) легко составляется матрица, для которой p(λ ) является характеристическим многочленом,
задача вычисления собственных значений матрицы математически
эквивалентна задаче вычисления корней многочлена [4, 6].
На первом этапе в практической деятельности спектральные задачи решались методами, основанными на построении характеристических многочленов. С точки зрения выполнения минимального числа арифметических операций подобный подход был сравнительно
эффективен. Тем не менее, на рубеже 1940–1950-х годов произошла
переоценка применявшихся методов решения спектральных задач.
Было установлено [6, 14], что даже для благоприятных ситуаций, например, для многочленов с действительными хорошо разделенными корнями, корни многих многочленов весьма чувствительны к малым изменениям их коэффициентов. В то же время чувстви27
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
тельность собственных значений, являющихся гладкими функциями
элементов матрицы, к изменениям этих элементов совершенно
иная. Поскольку построение характеристических многочленов основывалось на неустойчивых алгоритмах, исходная спектральная
задача с приемлемой обусловленностью нередко заменялась на плохо обусловленную задачу вычисления корней многочлена, в которой корни многочлена оказывались чрезвычайно чувствительными
к незначительным изменениям его коэффициентов. На практике
указанная замена могла привести к существенным ошибкам при
нахождении собственных значений и собственных векторов [6].
Появление ЭВМ позволило существенно ускорить расчеты. В
свою очередь, усложнение решаемых задач потребовало уменьшения погрешности применяемых вычислительных методов. Поэтому разрабатываемые с тех пор методы решения задач на собственные значения не требуют построения характеристического многочлена. Реализация подобных методов, характеризующихся увеличением объема вычислений, позволила существенно повысить
точность решения и сделала возможным решение спектральных
задач значительно более высокого, чем прежде, порядка.
На данный момент разработано довольно большое число методов решения спектральных задач и их модификаций. Однако указать способ нахождения собственных значений и собственных векторов, который одинаково хорошо подходил бы для решения разнообразных практических задач, невозможно. Дело в том, что на эффективность и удобство решения спектральных задач влияет множество факторов: форма записи задачи на собственные значения,
тип матриц и их размерность, объем решаемой задачи (полная или
частичная проблема собственных значений), наличие кратных или
близких собственных значений, необходимая точность вычислений.
При выборе наиболее рациональных методов решения спектральных задач учитывались следующие требования:
1) обеспечение решения спектральной задачи для динамической системы произвольной структуры (в том числе систем с реактивными звеньями) без перепрограммирования алгоритма;
2) соответствие метода решения спектральной задачи размерности
динамической модели (n<150) , типу матриц, составляемых при анализе динамических систем трансмиссий КМ, и форме записи задачи;
3) высокая точность решения, достаточная для надежного вычисления собственных значений и собственных векторов динамических систем со многими степенями свободы;
28
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
4) высокая эффективность с точки зрения вычислительных затрат;
5) устойчивость алгоритма;
6) простота реализации алгоритма на ЭВМ;
7) наличие отлаженных программ для ЭВМ.
Метод цепных дробей и метод остатка, получившие широкое распространение в 1950-е годы для расчета собственных спектров динамических моделей силовых передач КМ, позволяют рассчитывать произвольное число низших собственных частот и соответствующих форм
колебаний, т. е. фактически решать не только частичную, но и полную
проблему собственных значений. Однако указанные методы неприменимы к динамическим системам со сложной структурой. Подробно с
методом цепных дробей можно ознакомиться в [10, 13].
Соответствие метода решения спектральных задач первому и
третьему требованиям привело к рассмотрению исключительно
матричных методов, обеспечивающих решение для произвольных
динамических систем и не требующих построения характеристического многочлена. Последнее, как было отмечено выше, позволяет получить необходимую точность расчетов для матриц высокого порядка.
Необходимость выполнения второго требования позволила исключить из рассмотрения методы решения больших спектральных задач с
n>500 и методы, предназначенные для работы с комплексными матрицами, и ограничиться методами, ориентированными на решение
стандартной и обобщенной проблемы собственных значений с действительными матрицами сравнительно небольшого порядка (n≤150) ,
допускающих их полное хранение в оперативной памяти ЭВМ.
Отметим несколько важных моментов, связанных с нахождением собственного спектра действительных матриц.
Решение обобщенной проблемы собственных значений (43) в
принципе труднее решения стандартной [9], поэтому обобщенная
проблема должна быть приведена к эквивалентной стандартной
форме. Однако подобное приведение требуется не всегда. Имеется
ряд методов, расширенных на случай обобщенной проблемы собственных значений. Их применение позволяет несколько упростить вычисление собственного спектра и сохранить разреженность исходных матриц A и C . Формирование матрицы A−1C
приводит к потере симметричности и разреженности матриц [9].
Из двух возможных эквивалентных стандартных форм (44) и
(46), к которым может быть преобразована (43), предпочтитель29
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ней является форма (46) с симметричной матрицей. Вычисление
собственного спектра действительных симметричных матриц
существенно проще решения той же задачи для матриц общего
вида, поскольку собственные значения первых всегда хорошо
обусловлены [6, 14].
Лучшие из имеющихся в настоящее время методов решения
стандартной алгебраической проблемы собственных значений для
неразреженных матриц требуют предварительного приведения к
более простой форме, например, к трехдиагональной форме для
симметричной матрицы и к верхней форме Хессенберга для несимметричных матриц. Приведение осуществляется при помощи
численно устойчивых преобразований подобия, поэтому ошибки
округления не являются серьезной проблемой.
Для удобства выбора рациональных способов решения спектральных задач в зависимости от формы записи задачи на собственные значения, типа матриц и объема необходимых вычислений
методы решения алгебраической проблемы собственных значений,
применимые к динамическим моделям силовых передач КМ, представлены в табл. 1.
Анализ потребностей, возникающих при динамических расчетах силовых передач КМ, позволил выделить четыре наиболее целесообразных метода решения алгебраической проблемы собственных значений.
При решении спектральных задач для динамических моделей
силовых передач КМ со многими степенями свободы требуется
знание лишь низших собственных частот и соответствующих
форм колебаний. Решение полной проблемы собственных значений в данном случае привело бы к увеличению объема и времени
вычислений. Поэтому для расчета низших собственных частот и
соответствующих форм колебаний консервативных линейных систем (9) и (46) рекомендуется метод Гивенса (метод деления пополам) или PWK-алгоритм, представляющий собой один из последних вариантов QL-алгоритма [11] и не содержащий операций вычисления квадратных корней, совместно с методом обратной итерации. Для решения аналогичной задачи применительно к неконсервативным линейным системам (3) и (48) следует воспользоваться QR-алгоритмом, предназначенным для действительных несимметричных матриц с комплексными в общем случае собственными значениями и собственными векторами, в сочетании с методом обратной итерации.
30
венные
значения
Часть собственных
значений и
собственных
векторов
Часть
собственных
значений
Все собственные значения и собственные
векторы
Все собст-
Объем
необходимых
вычислений
Метод Гивенса
–
–
–
Метод Гивенса Метод Гивенса
Метод Гивенса
или PWK-алгоритм
совместно
совместно с метосовместно с метос методом
дом обратных
дом обратных
обратных
итераций
итераций
итераций
QR-алгоритм
совместно
с методом
обратных
итераций
Метод Гивенса
Метод Гивенса
или PWK-алгоритм
Метод
Гивенса
QL-алгоритм
Метод Гивенса совместно с методом обратных
итераций
QL-алгоритм
QL-алгоритм
или PWK-алгоритм
QR-алгоритм
QL-алгоритм
QL-алгоритм
или
PWK-алгоритм
QL-алгоритм
СимметричСимметричные
ные ленточные трехдиагональные
матрицы
матрицы
QZ-алгоритм
QL-алгоритм или
PWK-алгоритм
Симметричные
матрицы
QR-алгоритм
Матрицы
общего вида
QL-алгоритм, PWKалгоритм или обобщенный метод Якоби
Симметричные
матрицы
QZ-алгоритм
Матрицы
общего
вида
Методы вычисления собственных значений и собственных векторов
Стандартная проблема собственных значений
Обобщенная проблема собственных
и собственных векторов Dx = λx
значений (C −λA) x = 0
Рекомендуемые методы решения спектральных задач при исследовании динамических процессов в силовых
передачах колесных и гусеничных машин
Таблица 1
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Отыскание нескольких собственных векторов применительно к
действительным симметричным и несимметричным матрицам
имеет свои особенности. В первом случае вычислению m собственных векторов предшествует расчет только тех m собственных
значений, для которых вычисляются собственные векторы, во втором − нахождение всего частотного спектра.
Оценку полного частотного спектра консервативных и неконсервативных систем рационально провести с использованием соответственно QL- (PWK-) и QR-алгоритмов.
Определенный интерес представляют методы, предназначенные для разреженных и ленточных матриц, поскольку динамические модели силовых передач КМ, как правило, описываются
разреженными матрицами. Например, цепная система при последовательной нумерации соединяемых масс в задаче (46) характеризуется трехдиагональной матрицей.
4. ПРОГРАММНЫЕ СРЕДСТВА ДЛЯ РЕШЕНИЯ ЗАДАЧ
НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ
Самостоятельная разработка программ для решения спектральных задач весьма сложна и трудоемка, поэтому в повседневной вычислительной практике целесообразно использовать стандартные подпрограммы. Применение стандартных подпрограмм,
предназначенных для решения наиболее часто встречающихся математических задач, позволяет существенно упростить и ускорить
процесс программирования, повысить надежность и эффективность расчетов. В этом случае разработка всей программы сводится к написанию и отладке сравнительно небольшой программы,
осуществляющей ввод исходных данных, формирование и преобразование матриц, обращение к стандартным подпрограммам или
функциям и вывод результатов расчета.
Решению задач на собственные значения посвящен классический
справочник [15], в котором приводятся исходные тексты 84 подпрограмм, написанных на языке АЛГОЛ-60 и охватывающих важнейшие
разделы линейной алгебры: преобразование матриц, решение систем
линейных уравнений и задач на собственные значения. В справочнике даются описания подпрограмм, краткие теоретические положения,
раскрывающие основные идеи предлагаемых алгоритмов, область
применения алгоритмов, оценки точности вычислений и примеры
решения тестовых задач. Однако с появлением новых алгоритмиче32
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ских языков и, в частности, языка ФОРТРАН, ставшего в последующем основным для решения научных и технических задач, АЛГОЛ
постепенно вышел из употребления.
Обширная библиотека стандартных и научных подпрограмм была разработана и составлена на языке ФОРТРАН, который продолжает широко применяться в наши дни. К сожалению, большинство математических программ, поставлявшихся вместе с ЭВМ, например с
ЭВМ единой серии (ЕС ЭВМ), хранилась в виде объектных кодов,
т. е. на языке команд самой ЭВМ. С переходом с ЕС ЭВМ на ПЭВМ
доступ к библиотекам стандартных научных подпрограмм, хранившихся в объектных кодах, оказался практически невозможным. Тем
не менее, с исходными текстами некоторых из программ, предназначенных для решения спектральных задач, можно ознакомиться в
сборнике [12], опубликованном на русском языке.
Подпрограммы NROOT и EIGEN [12] позволяют решать полную проблему собственных значений методом Якоби для задачи
(44) с одновременным определением всех собственных значений и
собственных векторов. Машинная реализация метода Якоби не
требует предварительных преобразований матриц жесткостных и
инерционных коэффициентов, что упрощает применение стандартных программ. Наличие вариантов программ с удвоенной
точностью вычислений позволяет рассчитывать собственные значения и собственные векторы с необходимой точностью для динамических моделей со многими степенями свободы.
Подпрограммы HSBG и ATEIG [12] позволяют рассчитывать все
собственные значения для действительных несимметричных матриц
в задачах (44) и (48). Первая из этих подпрограмм осуществляла приведение вещественной матрицы к верхней форме Хессенберга, вторая
– нахождение всех собственных значений матрицы в верхней форме
Хессенберга. Однако указанные подпрограммы не позволяют рассчитывать необходимые собственные векторы исследуемых систем. Для
получения необходимых собственных векторов подпрограммы HSBG
и ATEIG должны быть дополнены подпрограммами:
– масштабирования исходной матрицы;
– запоминания информации о преобразованиях уравновешенной матрицы в верхнюю форму Хессенберга;
– вычисления необходимых собственных векторов матрицы в
верхней форме Хессенберга методом обратной итерации;
– преобразования вычисленных собственных векторов матрицы в верхней форме Хессенберга в собственные векторы уравновешенной матрицы;
33
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
– преобразования собственных векторов уравновешенной матрицы в собственные векторы исходной матрицы.
Подпрограмма масштабирования предназначена для уменьшения ошибок, возникающих при работе QR-алгоритма и обусловленных округлением или способом окончания процесса
итераций.
Более подробно с вычислением собственных векторов действительной несимметричной матрицы можно ознакомиться в [15].
Подпрограммы HSBG и ATEIG не имеют вариантов двойной
точности и, следовательно, при расчетах собственных значений с
симметричными матрицами A и C в задаче (44) обеспечивают
меньшую точность, чем при использовании подпрограмм NROOT
и EIGEN в варианте с удвоенной точностью.
Среди зарубежных изданий в первую очередь следует отметить
руководства [17, 18] пакета EISPACK, содержащего один из лучших наборов программ на ФОРТРАНе, предназначенных для решения задач на собственные значения.
Кроме того, высококачественные программы для вычисления
собственных значений и собственных векторов включены в пакет
IMSL, поставляемый вместе с профессиональной версией Fortran
PowerStation 4.0 и Compag Visual Fortran. Ранние версии пакета
IMSL содержали подпрограммы, представлявшие собой прямой
перевод с АЛГОЛа на ФОРТРАН некоторых из программ справочника [15]. Однако более поздние варианты программ пакета
IMSL разработаны с учетом последних достижений в области
матричных вычислений, произошедших с момента выхода пакета
EISPACK. В частности, для решения ряда спектральных задач в
пакете IMSL предлагаются программы, составленные на основе
PWK-алгоритма.
При подборе математических программ из пакетов EISPACK и
IMSL для решения задач на собственные значения в соответствии
с табл. 1 возможны затруднения, обусловленные терминологической путаницей. Это связано с тем, что QL-алгоритм, являющийся
модификацией QR-алгоритма, в ряде литературных источников
называется QR-алгоритмом. В частности, в документации на английском языке, поставляемой вместе с пакетом Visual Fortran 6.5,
и ее переводе на русский язык [2] QR-алгоритм называется PWKалгоритмом. Поэтому для удобства пользователей приводятся
табл. 2 и 3, в которых предлагается список подпрограмм из пакетов EISPACK и IMSL, предназначенных для решения спектраль34
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ных задач в зависимости от типа задачи на собственные значения,
типа матриц и объема необходимых вычислений.
В пакете EISPACK для решения большинства типов задач на
собственные значения предусмотрено последовательное применение нескольких подпрограмм, выполняющих отдельные законченные этапы вычислений. Некоторые из этих подпрограмм осуществляют приведение матриц к более простой форме − трехдиагональной форме или верхней форме Хессенберга, другие определяют собственные значения и (или) собственные векторы, третьи
восстанавливают собственные векторы промежуточных или исходных матриц и т. д. Решение спектральных задач с использованием пакета IMSL несколько проще, поскольку для этого необходим вызов только одной подпрограммы.
В последние годы широкое распространение получают интерактивные системы для выполнения инженерных и научных расчетов, например системы MATLAB и MathCAD. В системе
MATLAB 5.х имеется обширная библиотека математических
функций [11]. Применяющееся программное обеспечение позволяет решать стандартную и обобщенную, полную и частичную проблему собственных значений для матриц общего вида. Указанные
задачи решаются с помощью математических функций eig и eigs.
Первая из них позволяет решать полную проблему собственных
значений, вторая − частичную. Обе программы разработаны на
основе QR- и QZ-алгоритмов и метода Арнольди. QZ-алгоритм
представляет собой расширение QR-алгоритма для решения обобщенной проблемы собственных значений (43) с несимметричными
матрицами [19]. Метод Арнольди предназначен для решения спектральных задач с несимметричными разреженными матрицами,
размерность которых на порядки превышает размерность динамических моделей силовых передач КМ с сосредоточенными параметрами. Универсальные функции пакета MATLAB могут быть
использованы для расчета собственного спектра силовых передач.
В более поздних версиях системы MATLAB, например системы
MATLAB 6.5, как это следует из текста встроенной справки, используются стандартные программы пакета LAPACK [16], учитывающие тип матриц.
Аналогичные программы имеются в системе MathCAD (табл. 4),
однако сведения о методах, заложенных в основу этих программ, отсутствуют.
35
Программы для решения задач на собственные значения
Обобщенная проблема собственных
Стандартная проблема собственных значений
Объем необзначений (C −λA) x = 0
и собственных векторов Dx = λx
ходимых
Симметричные матриСимметричСимметричные
вычислений
Матрицы
Матрицы
Симметричцы (А – положительно
ные ленточные трехдиагональобщего вида
общего вида
ные матрицы
определенная)
матрицы
ные матрицы
BALANC,
Все собственQZHES,
REDUC,
ELMHES,
TRED2,
BANDR,
ные значения
QZIT,
TRED2,
ELTRAN,
IMTQL2
и собственные
QZVAL,
TQL2,
TQL2
TQL2
HQR2,
векторы
QZVEC
REBAK
BALBAK
QZHES,
REDUC,
BALANC,
Все собственTRED1,
BANDR,
IMTQL1
QZIT,
TRED1,
ELMHES,
ные значения
TQL1
TQLRAT
QZVAL
TQLRAT
HQR
Часть собстBALANC,
REDUC, TRED1,
TRED1,
BANDR,
венных значеELMHES,
BISECT,
BISECT,
BISECT,
ний и собст–
HQR, INVIT,
BISECT,
TINVIT, TRBAK1,
TINVIT,
TINVIT
венных вектоELMBAK,
BANDV
REBAK
TRBAK1
ров *
BALBAK
Часть cобстREDUC,
TRED1,
BANDR,
–
–
BISECT
венных значеTRED1,
BISECT
BISECT
ний
BISECT
*Для матриц общего вида вычислению m собственных векторов предшествует расчет всего частотного спектра.
Программы пакета EISPACK, рекомендуемые для решения спектральных задач при исследовании
динамических процессов в силовых передачах колесных машин
Таблица 2
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
GVCRG
GVLRG
–
–
Все собственные значения
и собственные
векторы
Все собственные значения
Часть собственных значений
и собственных
векторов
Часть собственных значений
Объем необходимых вычислений
–
–
GVLSP
GVCSP
–
–
EVLRG
EVCRG
EVBSF
EVFSF
EVLSF
EVCSF
EVBSB
EVFSB
EVLSB
EVCSB
EVBSB
EVFSB
EVLSB
EVCSB
Программы для решения задач на собственные значения
Обобщенная проблема собстСтандартная проблема собственных значений
венных значений (C −λA) x = 0
и собственных векторов Dx = λx
Симметричные
Симметричные
Симметричные
Матрицы
Матрицы
Симметричные
матрицы (А –
трехдиаголенточные
общего вида
положительно
общего вида
матрицы
нальные
матрицы
определенная)
матрицы
Программы пакета IMSL, рекомендуемые для решения спектральных задач при исследовании
динамических процессов в силовых передачах колесных машин
Таблица 3
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Входные параметры
Квадратная матрица A
Квадратная матрица A
Квадратная матрица A и заданное собственное значение λ ,
полученное с помощью функции
eigenvals
Назначение
Нахождение полного собственного
спектра в стандартной задаче на
собственные значения Ax = λx
Нахождение всех собственных
векторов в стандартной задаче на
собственные значения Ax = λx
Нахождение только одного собственного вектора для заданного
собственного значения λ в стандартной задаче на собственные
значения Ax = λx
Функция
eigenvals ( A)
eigenvecs( A)
eigenvec( A, λ)
функции eigenvals
щью функции eigenvals
Функция возвращает вектор, представляющий собой собственный вектор*,
соответствующий i-му собственному
значению λ i , полученному с помощью
Функция возвращает вектор, координатами которого являются собственные значения
Функция возвращает квадратную матрицу
такой же размерности, что и у входной
матрицы A; i-й столбец возвращаемой
матрицы представляет собой собственный
вектор*, соответствующий i-му собственному значению λ i , полученному с помо-
Выходные параметры
Функции пакета MathCAD, предназначенные для решения задач на собственные значения
Таблица 4
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Квадратные матрицы A и C
Квадратные матрицы A и C
Нахождение всех собственных
векторов в обобщенной задаче
на собственные значения
(C − λA) x = 0
genvals (C , A)
genvecs (C , A)
ченному с помощью функции genvals
Функция возвращает квадратную матрицу
такой же размерности, что и у входной
матрицы A или C ; i-й столбец возвращаемой матрицы представляет собой
собственный вектор*, соответствующий
i-му собственному значению λ i , полу-
Функция возвращает вектор, координатами которого являются собственные
значения
Выходные параметры
* Все вычисленные собственные векторы нормированы: x т x = 1, где x т − транспонированный собственный вектор x.
Входные параметры
Назначение
Нахождение полного собственного
спектра в обобщенной задаче на
собственные значения
(C − λA) x = 0
Функция
Окончание табл. 4
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Таким образом, при решении задач на собственные значения
имеется выбор: разрабатывать соответствующую программу на
языке ФОРТРАН с использованием стандартных подпрограмм пакетов EISPACK и IMSL или программировать в интерактивных системах MATLAB и MathCAD. При любом из указанных подходов
обеспечивается доступ к библиотекам стандартных математических
подпрограмм, позволяющих решать самые разнообразные математические задачи. Главным преимуществом интерактивных систем
является простота применения и контроля за результатами промежуточных вычислений, наглядность математических действий и
возможность символьных вычислений. Поэтому эти системы и, в
первую очередь, система MathCAD могут быть рекомендованы для
решения учебных и сравнительно простых научных и инженерных
задач, для которых время вычислений не имеет большого значения.
Разрабатывать программы на языке ФОРТРАН, по всей видимости,
наиболее целесообразно при решении сложных или специализированных задач, когда необходим быстрый счет и, кроме того, алгоритм их выполнения требует сложных логических переходов.
5. ПРИМЕРЫ РЕШЕНИЯ СПЕКТРАЛЬНЫХ ЗАДАЧ
С ИСПОЛЬЗОВАНИЕМ ПАКЕТА MathCAD
Пример 1. Для динамической системы с тремя степенями свободы с параметрами: J1 = 1, J 2 = 2, J 3 = 3, кг ⋅ м 2 , k12 = 1,
k23 = 2, Н ⋅ м, найти все собственные частоты и формы колебаний
консервативной системы тремя способами в результате решения:
– обобщенной алгебраической проблемы (43);
– стандартной алгебраической проблемы (46) с симметричной
матрицей;
– стандартной алгебраической проблемы (48) с несимметричной матрицей.
Решение.
Способ 1. Решение обобщенной алгебраической проблемы (43)
1. Формирование исходных матриц:
 1 0 0


A :=  0 2 0 


 0 0 3 
40
 1 −1
0 


C := −1 3 − 2


 0 − 2
2 
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
2. Решение обобщенной алгебраической проблемы собственных значений (43):
а) вычисление собственных значений:
2.295


λ CA = 0.871


 0 
λ CA := genvals (C, A)
б) вычисление собственных векторов:
mz: = genvecs (C, A)
 0.581 0.916 0.577 


mz = −0.753 0.118 0.577 


 0.308 − 0.384 0.577
3. Расчет относительных координат собственных векторов (относительно амплитуды колебаний первой массы):
0
z0 :=
mz
mz 0,0
1
z1:=
mz
mz 0,1
2
z2 :=
mz
mz 0,2
4. Вычисление круговых частот колебаний:
ωCA := λ CA
1.515 


ωCA = 0.933


 0 
5. Результаты решения.
Круговая частота свободных
колебаний, с−1
Соответствующая форма
свободных колебаний
ωCA0 = 1.515
 1 


z0 = −1.295


 0.53 
41
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ωCA1 = 0.933
 1 


z1 =  0.129 


−0.419
ωCA2 = 0
 1 
 
z2 =  1 
 
 1 
Способ 2. Решение стандартной алгебраической проблемы (46)
с симметричной матрицей
1. Формирование исходных матриц:
 1 0 0


A :=  0 2 0 


 0 0 3 
 1

C := −1

 0
−1
0 

3 − 2 

−2
2
2. Разложение матрицы A в виде A = LLТ (45) – разложение
Холецкого:
L: = cholesky(A)
1
0


L =  0 1.414

 0
0
0 

0

1.732
3. Формирование симметричной матрицы F в (46):
F: = (L)
−1
−1 T
⋅ C ⋅ (L) 



1 − 0.707


F =  − 0.707
1.5

0 − 0.816

0 

− 0.816 

0.667 
4. Решение стандартной алгебраической проблемы собственных значений (46):
а) вычисление собственных значений симметричной матрицы F:
42
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
0.871


λ F =  0 


2.295
λ F := eigenvals (F)
б) вычисление собственных векторов симметричной матрицы F:
MVF : = eigenvecs(F)
−0.801 0.408

MVF = −0.146 0.577

 0.581 0.707
0.439 

− 0.803

0.403 
5. Восстановление собственных векторов исходной системы в соответствии с (47):
Т −1
( )
mx: = L
MVF
−0.801 0.408 0.439 


mx = −0.103 0.408 − 0.568


 0.336 0.408 0.233 
6. Расчет относительных координат собственных векторов (относительно амплитуды колебаний первой массы):
0
x0 :=
mx
mx 0,0
1
x1:=
mx
mx 0,1
2
x2 :=
mx
mx 0,2
7. Вычисление круговых частот колебаний:
ωF := λ F
 0.933 


ωF = 1.369 ×10−8 


 1.515 
8. Результаты решения.
43
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Круговая частота свободных
колебаний, с−1
Соответствующая форма
свободных колебаний
 1 


x0 =  0.129 


−0.419
ωF0 = 0.933
 1 
 
x1 =  1 
 
 1 
ωF1 = 1.369 ×10−8
 1 


x2 = −1.295


 0.53 
ωF2 = 1.515
Способ 3. Решение стандартной алгебраической проблемы
(48) с несимметричной матрицей
1. Формирование исходных матриц:
 1 0 0


A :=  0 2 0 


 0 0 3 
 0 0 0


B :=  0 0 0


 0 0 0
 1 −1
0 


C := −1 3 − 2


 0 − 2
2 
2. Задание нулевой О и единичной I матриц:
 0 0 0


O :=  0 0 0


 0 0 0
 1 0 0


I :=  0 1 0


 0 0 1
3. Формирование матрицы D в (48):
а) cлияние матриц по горизонтали:
44
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 0 0 0 1 0 0


D1 =  0 0 0 0 1 0


 0 0 0 0 0 1
D1: = augment (O, I)
(
)
D2: = augment −A−1 ⋅ C, − A−1 ⋅ B

1
0
0 0 0
 −1

D2 = 0.5 −1.5
1
0 0 0


 0 0.667 − 0.667 0 0 0 
б) cлияние матриц D1 и D2 по вертикали:
D: = stack ( D1, D2)
 0
0
0

 0
0
0

 0
0
0
D = 
1
0
−1

1
0.5 −1.5

 0 0.667 − 0.667
1 0 0 

0 1 0

0 0 1 

0 0 0 

0 0 0 

0 0 0
4. Решение стандартной алгебраической проблемы собственных значений (48):
а) вычисление собственных значений:
λ D := eigenvals ( D)
 1.515i 


−1.515i



0 
λ D = 

 0.933i 


−0.933i


0 

45
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
б) вычисление собственных векторов:
MX D : = eigenvecs(D)
5. Расчет относительных координат собственных векторов (относительно первых координат собственных векторов):
i
i := 0.. 5
MX D :=
i
MX D
MX D0,i
Матрица собственных векторов:

1
1

−1.295 −1.295


0.53
 0.53
MX D =
 1.515i −1.515i

−1.962i 1.962i

 0.803i −0.803i

1
1
1
1
0.129
0.129
1
−0.419
−0.419
0
0.933i
−0.933i
0
0.12i
−0.12i
0
−0.391i 0.391i



1


1


−15 
− 2.803×10 


0

−15 

−2.39×10

1
6. Расчет круговых собственных частот:
ωDi := Im (λ Di )
7. Формирование векторов, содержащих коэффициенты форм
колебаний:
x Di
8. Результаты решения.
46
MX 

D0,i 


:= MX D1,i 


MX D 


2,i 
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Круговая частота свободных
колебаний, с−1
ωD0 = 1.515
ωD1 = 1.515
ωD2 = 0
ωD3 = 0.933
ωD4 = 0.933
ωD5 = 0
Соответствующая форма
свободных колебаний
x D0
 1 


= −1.295


 0.53 
x D1
 1 


= −1.295


 0.53 
x D2
 1 
 
=  1 
 
 1 
x D3
 1 


=  0.129 


−0.419
x D4
 1 


=  0.129 


−0.419
x D5
 1 
 
=  1 
 
 1 
Представленный пример показывает, что для выбранной динамической системы все круговые собственные частоты и соответствующие формы колебаний, найденные тремя разными способами,
одинаковы.
Очевидно, что с точки зрения минимизации затрат на программирование в среде MathCAD наименее трудоемким является
первый способ решения, а наиболее трудоемким − третий.
47
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Результаты решения третьим способом позволяют наглядно
подтвердить ряд теоретических положений, изложенных выше.
Во-первых, для системы без демпфирования действительная часть
собственных значений равна нулю, что соответствует незатухающим колебаниям. Во-вторых, последние n координат произвольного i-го собственного вектора X D i , являющегося i-м столбцом
матрицы MX D , получаются из первых n координат того же вектора путем их умножения на i-е собственное значение λ i , причем
первые n координат вектора X Di характеризуют соответствующую форму колебаний. В-третьих, для системы без демпфирования относительные фазовые сдвиги ϕrs равны нулю, и потому относительные формы колебаний являются действительными.
Пример 2. В динамическую систему из примера 1 добавим
демпфирование: f12 = 0,01 , f 23 = 0,04 , Н⋅м⋅с. Частоты и формы
колебаний неконсервативной системы найдем двумя способами в
результате решения:
– стандартной алгебраической проблемы (48);
– обобщенной алгебраической проблемы (49).
Решение.
Способ 1. Решение стандартной алгебраической проблемы (48)
1. Формирование исходных матриц:
 1 0 0


A :=  0 2 0 


 0 0 3 
 0.01 − 0.01
0 


B := − 0.01
0.05 − 0.04



0 − 0.04 0.04 


 1 − 1 0
C := −1 3 − 2


 0 − 2 2 
2. Задание нулевой O и единичной I матриц:
 0 0 0


O :=  0 0 0


 0 0 0
 1 0 0


I :=  0 1 0


 0 0 1
3. Формирование матрицы D в (48):
а) слияние матриц по горизонтали:
48
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
D1: = augment (O, I)
(
)
D2: = augment −A−1 ⋅ C, − A−1 ⋅ B
б) слияние матриц D1 и D2 по вертикали:
D: = stack ( D1, D2)
 0
0
0

 0
0
0

 0
0
0
D = 
1
0
−1

1
0.5 − 1.5

 0 0.667 − 0.667
0 

0
1
0 

0
0
1 

− 0.01
0.01
0 

0.005 − 0.025
0.02 

0
0.13 − 0.013
1
0
4. Решение стандартной алгебраической проблемы собственных значений (48):
а) вычисление собственных значений:
λ D := eigenvals ( D)
б) вычисление собственных векторов:
MX D : = eigenvecs(D)
5. Расчет относительных координат собственных векторов (относительно первых координат собственных векторов):
i := 0…5
MX D
i
i
MX D
:=
MX D0,i
6. Расчет круговых собственных частот:
ωDi := Im (λ Di )
49
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
7. Формирование векторов, содержащих коэффициенты форм
колебаний:
x Di
MX 

D0,i 


:= MX D1,i 


MX D 


2,i 
8. Результаты решения.
Круговая частота свободных Соответствующая
свободных колебаний
колебаний, с−1
ωD0 = 1.515
ωD1 = 1.515
ωD2 = 0.933
ωD3 = 0.933
−9
ωD4 = 7.449×10
50
x D0


1


= −1.295 − 0.019i


 0.53 + 0.013i 
x D1


1


= −1.295 + 0.019i


 0.53 − 0.013i 
x D2


1


=  0.129 − 0.004i 


−0.419 + 0.002i
x D3


1


=  0.129 + 0.004i 


−0.419 − 0.002i
x D4
 1 
 
=  1 
 
 1 
форма
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 1 
 
ωD5 = −7.449 ×10
x D5 =  1 
 
 1 
Способ 2. Решение обобщенной алгебраической проблемы (49)
1. Формирование исходных матриц:
−9
 1 0 0

 1 −1
0 
0 
 0.01 − 0.01





A :=  0 2 0  B := − 0.01 0.05 − 0.04 C := −1 3 − 2






 0 − 0.04 0.04 
 0 − 2
 0 0 3 
2 
2. Задание нулевой O и единичной I матриц:
 0 0 0


O :=  0 0 0


 0 0 0
 1 0 0


I :=  0 1 0


 0 0 1
3. Формирование матрицы P в (49):
а) слияние матриц по горизонтали:
P1: = augment (O, А )
P2: = augment (−C, O)
б) слияние матриц P1 и P2 по вертикали:
 0
0

 0
0

 0
0
P: = stack ( P1, P2) P= 
1
−1

 1 − 3

2
 0
0
0
1
0
0
2
0
0
0
0
0
0
2
−2
0
0
0
0
0

0

3

0

0

0
4. Формирование матрицы S в (49):
51
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а) слияние матриц по горизонтали:
S1: = augment ( A, O)
S2: = augment (B, A )
б) слияние матриц S1 и S2 по вертикали:

1
0
0


0
2
0


0
0
3
S: = stack ( S1, S2) S = 
0
 0.01 − 0.01

 −0.01 0.05 − 0.04

0 − 0.04
0.04

0
0
0
0
0
0
1
0
0
2
0
0
0

0

0

0

0

3
5. Решение обобщенной алгебраической проблемы собственных значений (49):
а) вычисление собственных значений:
λ PS := genvals ( P, S)
б) вычисление собственных векторов:
MX PS : = genvecs(P, S)
6. Расчет относительных координат собственных векторов (относительно первых координат собственных векторов):
i := 0.. 5
i
MX PS :=
7. Расчет круговых собственных частот:
ωPSi := Im (λ PSi )
52
i
MX PS
MX PS0,i
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
8. Формирование векторов, содержащих коэффициенты форм
колебаний:
x PSi
MX


PS0,i 



:= MX PS1,i 


MX PS 


2,i 
9. Результаты решения.
Круговая частота свободных
колебаний, с−1
Соответствующая форма
свободных колебаний
x PS0


1


= −1.295 − 0.019i


 0.53 + 0.013i 
ωPS1 = 1.515
x PS1


1


= −1.295 + 0.019i


 0.53 − 0.013i 
ωPS2 = 0
 1 


x PS2 =  0.999


 1 
ωPS0 = 1.515
ωPS3 = 0.933
ωPS4 = 0.933
x PS3


1


=  0.129 + 0.02i 


−0.419 − 0.013i
x PS4


1



=  0.129 − 0.02i 


−0.419 + 0.013i
53
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ωPS5 = 0
x PS5
 1 


=  0.999


 1 
Результаты расчета показывают, что круговые собственные
частоты и формы колебаний, соответствующие круговой частоте
ω = 1,515, найденные двумя разными способами, как это и должно
быть, одинаковы. Однако формы колебаний для круговой частоты
ω = 0,933 несколько различаются, что обусловлено ошибками округления.
Введение демпфирования в систему привело к появлению ненулевых фазовых сдвигов ϕrs , из-за чего относительные формы
колебаний приобрели комплексный вид.
Пример 3. В динамическую систему из примера 1 введем
демпфирование, пропорциональное жесткости упругих участков,
соединяющих массы
fij = lkij ,
где l = 0,01 − коэффициент пропорциональности, с−1.
Частόты и формы колебаний неконсервативной системы с пропорциональным демпфированием найдем в результате решения
стандартной алгебраической проблемы (48).
Решение.
1. Формирование исходных матриц:
 1 0 0


A:= 0 2 0


 0 0 3 
B := 0.01⋅ C
 1

C:=−1

 0

 0.01
B=− 0.01


0
−1
3
−2
−0.01
0 

0.03 −0.02 

−0.02
0.02
2. Задание нулевой O и единичной I матриц:
54
0

− 2 

2
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 0 0 0


O :=  0 0 0


 0 0 0
 1 0 0


I :=  0 1 0


 0 0 1
3. Формирование матрицы D в (48):
а) слияние матриц по горизонтали:
D1: = augment (O, I)
(
)
D2: = augment −A−1 ⋅ C, − A−1 ⋅ B
б) слияние матриц D1 и D2 по вертикали:
D: = stack ( D1, D2)
4. Решение стандартной алгебраической проблемы собственных значений (48):
а) вычисление собственных значений:
λ D := eigenvals ( D)
б) вычисление собственных векторов:
MX D : = eigenvecs(D)
5. Расчет относительных координат собственных векторов (относительно амплитуды колебаний первой массы):
i := 0.. 5
i
MX D :=
i
MX D
MX D0,i
6. Расчет круговых собственных частот:
ωDi := Im (λ Di )
55
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
7. Формирование векторов, содержащих коэффициенты форм
колебаний:
x Di
MX 

D0,i 


:= MX D1,i 


MX D 


2,i 
8. Результаты решения.
Круговая частота свободных
колебаний, с−1
ωD0 = 1.515
ωD1 = 1.515
ωD2 = 0
ωD3 = 0.933
ωD4 = 0.933
56
Соответствующая форма
свободных колебаний
x D0
 1 


= −1.295


 0.53 
x D1
 1 


= −1.295


 0.53 
x D2
 1 
 
=  1 
 
 1 
x D3
 1 


=  0.129 


−0.419
x D4
 1 


=  0.129 


−0.419
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ωD5 = 0
x D5
 1 
 
=  1 
 
 1 
Следует отметить, что для силовых передач КМ не характерно
пропорциональное демпфирование, задаваемое в виде B = rA + tC ,
где A, C – матрицы инерционных и жесткостных коэффициентов; r , t − коэффициенты пропорциональности. Тем не менее,
данный пример приводится для подтверждения известного теоретического положения: формы колебаний неконсервативных
систем с пропорциональным демпфированием имеют действительный вид [3].
57
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
СПИСОК ЛИТЕРАТУРЫ
1. Аринчев С.В. Теория колебаний неконсервативных систем. М.: Изд-во
МГТУ им. Н.Э. Баумана, 2002. 464 с.
2. Бартеньев О.В. Фортран для профессионалов. Математическая
библиотека IMSL: Ч. 1. М.: ДИАЛОГ-МИФИ, 2000. 448 с.
3. Бидерман В.Л. Теория механических колебаний. М.: Высш. шк.,
1980. 408 с.
4. Вибрации в технике: Справ.: В 6 т. / Ред. совет: В.Н. Челомей
(пред.) и др. М.: Машиностроение, 1978. Т. 1: Колебания линейных систем / Под ред. В.В. Болотина. 352 с.
5. Гантмахер Ф.Р. Лекции по аналитической механике. М.: Наука,
1966. 300 с.
6. Икрамов Х.Д. Несимметричная проблема собственных значений.
Численные методы. М: Наука, 1991. 240 с.
7. Курош А.Г. Курс высшей алгебры. М: Наука, 1971. 432 с.
8. Мандельштам Л.И. Лекции по колебаниям. М.: Наука, 1972. 470 с.
9. Парлетт Б. Симметричная проблема собственных значений: Пер.
с англ. М.: Мир, 1983. 383 с.
10. Полунгян А.А., Фоминых А.Б. Методические указания к выполнению домашнего задания по курсу «Динамика колесных машин» / Под
ред. А.А. Полунгяна. М.: Изд-во МГТУ им. Н.Э. Баумана, 1993. 44 с.
11. Потемкин В.Г. Система инженерных и научных расчетов
MATLAB 5.x: В 2 т. М.: ДИАЛОГ-МИФИ, 1999. Т. 1. 366 с.; Т. 2. 304 с.
12. Сборник научных подпрограмм на Фортране / Пер. с англ.; Под.
ред. С.Я. Виленкина. Мн.: Статистика, 1974. Вып. 2. 224 с.
13. Терских В.П. Метод цепных дробей: В 2 т. М.–Л.: Машгиз, 1955.
Т. 1. 376 с.
14. Уилкинсон Дж. Алгебраическая проблема собственных значений:
Пер. с англ. М.: Наука, 1970. 564 с.
15. Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке
АЛГОЛ. Линейная алгебра: Пер. с англ. М.:Машиностроение, 1976. 389 с.
16. LAPACK User’s Guide. Anderson E. et al. Third Edition. Philadelphia: SIAM, 1999.
58
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
17. Matrix Eigensystem Routines − EISPACK Guide. V. 6. Smith B.T.
et al. Berlin; ⋅ Heidelberg New York: Springer-Verlag. 1974.
18. Matrix Eigensystem Routines − EISPACK Guide Extension. V. 51 /
Garbow B.S. et al. Berlin; Heidelberg; New York: Springer-Verlag. 1977.
19. Moler C.B., Stewart G.W. An Algoritm for Generalized Matrix Eigenvalue Problems // SIAM J. Nume. Anal. 1973. № 10. P. 241–256.
59
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ОГЛАВЛЕНИЕ
1. Уравнение свободных колебаний динамической системы
с конечным числом степеней свободы ..................................................
2. Эквивалентность нахождения собственных частот и решения
задачи на собственные значения ............................................................
2.1. Линейные консервативные системы ............................................
2.2. Линейные неконсервативные системы ........................................
3. Особенности расчета собственных значений и собственных
векторов ...................................................................................................
4. Программные средства для решения задач на собственные
значения ...................................................................................................
5. Примеры решения спектральных задач с использованием
пакета MathCAD ......................................................................................
Список литературы .....................................................................................
60
3
9
9
16
25
32
40
58
1/--страниц
Пожаловаться на содержимое документа