ПРИКЛАДНАЯ МАТЕМАТИКА И МЕХАНИКА
Том 78. Вып. 5, 2014
УДК 531.36:534.1
© 2014 г. Ж. Б. Бакиров, В. Ф. Михайлов
АНАЛИЗ НЕЛИНЕЙНЫХ СЛУЧАЙНЫХ КОЛЕБАНИЙ МЕТОДОМ УСРЕДНЕНИЯ
Рассматривается решение квазиконсервативной нелинейной колебательной системы, правые части которой пропорциональны малому параметру. Путем перехода к "медленным" переменным и сочетанием метода стохастического усреднения и теории марковских процессов получены основные соотношения для решения задачи. На базе быстрого преобразования Фурье разработан эффективный численный алгоритм, позволяющий получить плотность распределения выходных параметров и амплитуды колебаний. Рассмотрено приложение теории к решению уравнения Дуффин-га—Ван-дер-Поля при аддитивном и мультипликативном случайном воздействии.
Задачи о нелинейных случайных колебаниях имеют точное решение только для небольшого количества частных случаев. Обычно приходится довольствоваться приближенными методами, причем преобладают методы, заимствованные из нелинейной механики: метод линеаризации, метод малого параметра, стохастический аналог метода усреднения Крылова—Боголюбова. Наиболее полно теория нелинейных случайных колебаний изложена в известных монографиях
[1-3].
Для дифференциальных уравнений, правые части которых пропорциональны малому параметру и представляют собой периодические функции времени, справедлив принцип усреднения Крылова-Боголюбова, при помощи которого в сочетании с теорией марковских процессов можно получить эффективные решения нелинейных задач стохастической динамики.
Метод стохастического усреднения для нелинейных систем под воздействием белого шума был предложен вначале П.С. Ланда и Р. Л. Стратоновичем, развит Р.З. Хасьминским, позднее модифицирован [4]. Далее метод был распространен на широкополосные случайные воздействия М.Ф. Диментбергом [2, 5]. Позже [6] получено решение задачи для узкополосного воздействия, которое реализовано численно.
Метод стохастического усреднения используется для определения отклика нелинейных систем с одной степенью свободы в случаях, когда решения в аналитическом виде получить невозможно. Основные соотношения выведены в общем виде и представлены в форме, позволяющей разработать эффективную численную схему, использующей быстрое преобразование Фурье.
Рассмотрим нелинейную систему под действием мультипликативных и аддитивных случайных воздействий
m
x + f (x) = Vg (x, x) + V^ Z gk (x, x)n k (t) (1)
к=1
где ц — малый параметр, f (x) — нелинейная восстанавливающая сила, g(x, x) и gk(x, x) — периодические функции времени с периодом T (a) = 2л/ю(а), цк(t) — стационарные центрированные случайные процессы с матрицей корреляционных функций Kke(r).
Для малых ц движение системы близко к периодическому и решение может быть представлено в виде
x(t) = a(t) cos 9(t) + b, ф^) = \|/(t) + 0(t) (2)
где a(t) и 0(t) — "медленно" меняющиеся функции.
Введем новую переменную — энергию E(t), определенную как первый интеграл порождающего уравнения (ц = 0):
E = X2 /2 + V (х); V (х) = {f (z)dz (3)
о
где V(x) — потенциальная энергия системы, причем V(a + b) = E. Для порождающего уравнения a, 0 и b постоянные. Тогда, дифференцируя равенство (2) по времени, при учете последнего выражения получим
X(t) = -aP(a, ф) sin ф(0 (4)
Усредняя мгновенную частоту колебаний
d^= ^ [V(a + b) - V(a cosф + b)]1/2 = в(a,ф) (5)
dt a sin ф
можно записать
2п
ф(0 = (ü(a)t + 9; ю(а) = С0(а) = 1 f в(а, ф^ф (6)
п J
0
Подставляя выражения (2) и (4) в выражение (1), после преобразований получаем систему уравнений в стандартной форме
a = Ц qi(a, ф) + £ (a, ф)Пк
k=1 m
ф = yq2(a, ф) + 2k (a, ф)п
k=i (7)
m
k=1
здесь
a1k(a, ф) q1(a, ф) 7„. . . 7 a + b
1 = ч =-af P(a, ф)sin ф, f = f -—-gk (a, ф) g(a, ф) 1 + h
& 2k (a, ф) = q2(a, ф) = cos ф + h ^ = db
a1k q1 sin ф da
Поскольку функции q1(a, ф), aik(a, ф) периодические, для фиксированного а они могут быть разложены в ряд Фурье
да
q1(a, ф) = q1,0 + X (qí«C0S«9 + qí,si^)
n=1 (9)
да
Vik(a, ф) = <3¡ko + X (ст^созиф + ст^тиф)
n=1
Для автоматизации расчетов рекомендуется использовать быстрое преобразование Фурье. Тогда индексы r и i будут означать действительные и мнимые части разложения Фурье.
Для аналитического решения система уравнений в "стандартной форме" ряды (7) заменяются уравнениями Ито, являющимися результатом предельного перехода от некоторых стохастических разностных схем с независимыми нормальными прираще-
ниями [1, 2]. Решение уравнений Ито — предельный диффузионный марковский процесс, коэффициенты сноса и диффузии которого являются коэффициентами уравнения Ито и определяются через функции в системе (7). Вычисление коэффициентов сноса и диффузии марковского процесса зависит и от выбора разностной схемы. Здесь использована схема Стратоновича.
При малых значениях параметра ^ в правой части системы (7) справедлив метод стохастического усреднения, поэтому при вычислении коэффициентов предельного марковского процесса дополнительно проведено усреденение по времени. Для амплитуды уравнение Ито имеет вид
á = u(a) + a(a)í¡(t)
(10)
где £(г) — белый шум с единичной интенсивностью, и(а) и а(а) — усредненные коэффициенты сноса и диффузии, определяемые по формулам
l ^ m m
u(a) = J qx +1 j^X
2 k=1l =1 ^ da
d a
1k
da
1k
Ы + d0
Ы+T I Kkl(T)dT I
I с
(11)
u2(a) = J j H(aik|t aiit+x)Kki(T)dT ,
\-x k=11=i I
Индекс t означает, что берется значение функций в момент времени t.
Используя для ф(0 соотношение (6) и подставляя ряды (9) в соотношения (11), получим
т т
и(а) = дю + 0^1; оЯь, (0) +
k=11=1
m m
+zz
k=11=1
ПZ{(üík«)'a1ln + (a1kJG\ln + n(a1knOr2Jn -0Ík„0í2i„)}Skl(n^(a))
2 n=1
(12)
°2(a) = ZZ2no1ko01loSkl(0) +ZZ
k=11=1 k=11=1
nZ(c1kn^'\ln + ^1kn^1ln)Skl (n®(a)) . n=1
Штрихом обозначена производная по а; кроме того, учтено, что
да да
í Kkl(T)sin raxdx = 0, Skl(ю) = — [ Kkl(T) cos raxdx J 2 n J
—да —да
Стохастическому уравнению Ито (10) соответствует уравнение Колмогорова 2
dP = ~d [u(a)P] +-dT [a 2(a)P]
dt da
2da
где Р(а, г) — плотность распределения амплитуды. Стационарное решение этого уравнения имеет вид
P(a) = CQ(a); C =
\
-1
J Q(a)da
V o
-2/
, Q(a) = o (a) exp
2 [ ds X(s)
Постоянная С определена из условия нормировки.
да
Из уравнения (3) следует, что полная энергия Е, когда перемещение максимально, является функцией одной случайной переменной амплитуды. Из формулы преобразования плотности вероятности с использованием равенства (3) можно определить плотность распределения полной энергии
P(E) = P(a)/f (a + b); a + b = V _1(E) (14)
Было показано [2], что совместную плотность распределения перемещения и его производной можно выразить через плотность распределения полной энергии:
P(x, x) = P(E)/ T(E) (15)
Период изменения энергии T(E) получается из выражения T(a) = 2n/ro(a) с использованием соотношений (14) и (6). Теперь из равенства (16) можно определить одномерную плотность распределения перемещения
да
P(x) = J P(x, x)dx (16)
—да
Последовательность численного расчета. Исходя из вышеизложенного, предлагается следующий порядок расчета плотности распределения выходных параметров:
1) принимается множество значений a от нуля до заданного максимального значения на малом интервале (что обусловлено точностью);
2) выбираются значения ф на равных интервалах от 0 до 2п;
3) для каждого значения a из соотношения (6) с использованием выражения для ß(a, ф)определяем ra(a) и T (a);
4) для каждого значения a, используя быстрое преобразование Фурье определяем
?io, üii«, Ü12«, üiin и ci2n для n = 1,2,..., m (значение m должно выбираться в зависимости от отсечения частоты из функции спектральной плотности воздействия);
5) численным дифференцированием определяем производные (ü[kn)' и (a'ikn)';
6) по формулам (12) определяем коэффициенты сноса u(a) и диффузии a2(a);
7) по формуле (13) численным интегрированием вычисляем P (a);
8) для каждого a вычисляем E = V(a + b), f(a) и по формуле (14) вычисляем P(E);
9) из равенства (6) заменой a на E определяем a>(E) и T(E) = 2п/a>(E);
10) для разных сочетаний xи x по формуле (3) определяем Е и по формуле (15) вычисляем P(x, x);
11) численным интегрированием вычисляем P(x) по формуле (16). Пример. Рассмотрим нелинейное уравнение Дуффинга-Ван-дер-Поля
x + (ßi + ß2 x 2)x + asx + yx3 = xfi(t) + f2(i) (17)
где ю0, ß:, ß2 и у — постоянные, f1(() и f2(() — стационарные и эргодические процессы с нулевым математическим ожиданием и рациональной спектральной плотностью. Из уравнения (1) следует
f (x) = Юоx + yx\ g(x, x) = -(ßi + ß2x2)x
gi(x, -x)ni = xfl(t), g2(x, x)ll2 = f2(t)
x 2 4
V (x) = J f (u)du = «0 + Yr
Для функции V(x) из условия E = V(a + b) = V(-a + b) получено значение b = 0, и следовательно, h = 0. Подставив выражение для V(x) в равенство (6), получим
Р(А. Ф) - 1
( 2 4
2 a л 2 \ Ya 4 ч
ю0—(1 - cos ф) + J— (1 - cos ф)
1/2
a sin ф
Далее из соотношений (8) имеем
?1(a, Ф) = —г— (в1 + в 2a 2cos^)p2(a, ф>т2ф f (a)
Ф) = ü12(a, Ф) = --f- P(a, ф)sln ф cos ф f (a)
f (a) = ®0 + Ya3
Расчеты для уравнения (17) проведем, приняв для обоих случайных процессов одинаковый вид спектральной плотности:
£(ю) _ 2аст2_б2 + а2_
п (ю2 - 92 - а2)2 + 4а V
Меняя значение a, можно изменить форму спектральной плотности от узкополосной до широкополосной. Приняв
a = 0.35 с-1 , 9 = 2.45 с-1
получим узкополосное воздействие при Smax/a2 = 0.453 и ю = 2.42 с-1 . Выбирая подходящее значение дисперсии а2, можно интенсивность воздействия сделать достаточно малой.
При расчетах принято
Р1 = 0.01с-1, р2 = 0.02 м-2с-1, у = 1м-2с-2, с2 = 7.3 • 10-4, ю0 = 2.5 с-1
Результаты расчетов с применением ПК MATLAB приведены на фиг. 1 и 2.
На фиг. 1 слева показана плотность распределения амплитуды колебаний, полученная после реализации шага 7 вышеуказанного алгоритма расчета. Распределение од-номодально и симметрично относительно а = 1.4 м.
На фиг. 2 показана совместная плотность распределения перемещения и скорости. График получен после реализации шага 10 алгоритма. Он представляет собой седловидную поверхность в пространстве (впадина на кривой). Сеч
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.