Параметрический анализ устойчивости линейного осциллятора с мультипликативным шумом

В настоящее время во многих областях знаний большой интерес вызывают задачи, связанные с понятием колебаний или осцилляций различных динамических систем. Эти задачи возникают при наблюдении эволюции систем с течением времени. Например, эволюции простейших колебаний — колебания маятников, пружин с грузом, в радиотехнике это колебательный электрический контур — связанные конденсатор и катушка. Кроме того, поскольку временную эволюцию любой величины с помощью преобразования Фурье можно представить в виде суммы периодических членов, в настоящее время большое внимание уделяется изучению периодических явлений. С их рассмотрением тесно связано понятие аттрактора (притягивающего предельного множества в фазовом пространстве траекторий), которому соответствуют устойчивые периодические движения, и понятие устойчивого состояния равновесия — точки стабилизации фазовых траекторий. Известно, что и тому, и другому сопутствует общее сжатие (когда некоторая область в фазовом пространстве переходит с течением времени внутрь себя) и локальная устойчивость фазовых траекторий решения систем ОДУ, описывающих поведение соответствующей динамической системы. Выйти из этих состояний равновесия система может только в результате некоторого толчка извне, например, с появлением некоторого действующего шума достаточной интенсивности.

Появление неточности в измерениях сигналов различных приборов так же, как и влияние различных шумов, неизбежно влечет за собой применение в исследованиях такого мощного аппарата как стохастические дифференциальные уравнения. А вместе с тем возникают как вопросы устойчивости решений, так и вопросы качественного исследования влияния различных параметров системы на устойчивость решения СДУ.

Данный раздел посвящен теоретическому анализу устойчивости тривиального решения линейного СДУ второго порядка в смысле Ито с мультипликативным шумом. Это СДУ, описывающее осциллятор с линейным трением и флуктуирующей частотой [100], возникает практически во всех областях физики. Также это СДУ является линеаризацией двухмерной нелинейной системы СДУ в общем случае.

где Л, си, а — вещественные параметры, n (t) — процесс белого шума, т. е. гауссовский процесс, для которого Линейный осциллятор с мультипликативным шумом описывается стохастическим дифференциальным уравнением в смысле Ито второго порядка вида.

Если ввести новые переменные.

и добавить начальные данные тДО) = уо, Х2(0) = yi, то получим следующую задачу Коши для линейной системы СДУ в смысле Ито:

которую также можно записать в виде.

где X (t) = (a-'i (t),…, х_Пт (t))^T — п_х-мерный случайный процесс; f (X(t)) определенная на интервале [О, Т] измеримая п^-мерная векторфункция, S (X (t)) — определенная на интервале [О, Т] измеримая матричная функция размера п_х х n_w. W (t) — п_№-мерный стандартный винеровский процесс, А и В — матрицы размеров п_х х п_х и п_х х n_w соответственно.

или.

Нетрудно показать, что для задачи (2) СДУ в смысле Ито и в смысле Стратоновича совпадают.

Система (1.36) удовлетворяет условиям теоремы существования и единственности решения системы СДУ [79] (например, при К = max (l-t- 2А², 2u)⁴+ct²)). Следовательно, решение системы (1.36) существует, единственно и является марковским процессом.

Обозначим через M (t) = … т_П:г)'^! вектор первых моментов решения X (t) системы (1.38), а через T (i) — матрицу вторых моментов с элементами 7— 1 ,…, п_х. Тогда уравнения на первый и второй моменты случайного вектора X (t) будут иметь следующий вид [119]:

Так как система (1.40) является линейной системой первого порядка с постоянными коэффициентами, то ее решением, очевидно, будет.

Зная первый и второй моменты, можно выписать систему ОДУ для дисперсий компонент решения системы (1.38). Для этого введем матрицу

Следовательно,.

Окончательно получаем.

В силу того, что матрицы Г (?) и D (t) являются симметричными, то системы (1.41) и (1.43) для матриц (1.39) можно записать в виде.

Исследуем вопрос об устойчивости тривиального решения системы.

(1.36) по вероятности, в среднем и среднем квадратическом.

Определение 1.1 [145]. Решение X (t) = 0 системы уравнений (1.38) называется устойчивым по вероятности (при t > to), если для любых е > 0, 5 > 0 найдется г > 0:

Определение 1.2 [145]. Решение X (t) =0 системы уравнений (1.38) устойчиво в среднем, если для любых е > 0 существует г > 0:

Определение устойчивости нулевого решения в среднем квадратическом приведено в параграфе 2.4 учебного пособия [31].

Очевидно, что устойчивость тривиального решения системы (1.36) будет зависеть от соотношений параметров Л, ш, а. Нас будет интересовать, при каких условиях на параметры тривиальное решение системы.

(1.36) будет устойчивым по вероятности, в среднем и среднем квадратическом.

I. Устойчивость по вероятности. Для нахождения ограничений на параметры, при которых тривиальное решение системы (1.36) устойчиво по вероятности, воспользуемся следующей теоремой.

Теорема 1.1 [145]. Пусть {t > 0} х U = U, C^(U) — класс функций дважды, непрерывно дифференцируемых по X? U и один

раз по t всюду, кроме, быть может, множества X = 0. Пусть в области U, содержащей прямую X (t) = 0, существует непрерывная положительно определенная в смысле Ляпунова функция V (t, X (t))? С2(Ui), удовлетворяющая при X (t) ф 0 условию

Тогда тривиальное решение системы (1.38) устойчиво по вероятности. А если выполнено условие.

при |Х| > е, то нулевое решение асимптотически устойчиво по вероятности.

Здесь ai (t, X) — (AX (t))i — г-й элемент вектора AX (t), Bij (t, X) — элемент матрицы [BX (t)][BX (t)]^J.

В работе [134] для данной задачи найдена положительно определенная функция V и получены следующие условия устойчивости по вероятности тривиального решения системы (1.36):

Нетрудно показать, что тривиальное решение при данных параметрах будет даже асимптотически устойчиво.

II. Устойчивость в среднем. Устойчивость в среднем тривиального решения системы (1.38) эквивалентна устойчивости, но Ляпунову тривиального решения системы ОДУ на первый момент (1.40).

Для определения условий на параметры А, ш, <�т, при которых тривиальное решение системы (1.40) устойчиво, воспользуемся следующей теоремой.

Теорема 1.2 (Ляпунова) [145|. Рассмотрим систему

где, А = {a^fc} - постоянная матрица, ipi — бесконечно малые первого порядка. Тогда, если все собственные значения матрицы, А имеют отрицательные вещественные части, нулевое решение рассматриваемой системы асимптотически устойчиво-, если же хоть одно собственное значение имеет положительную вещественную часть, то нулевое решение неустойчиво.

Система (1−40) для матрицы (1.39) имеет вид.

т. е. является линейной с постоянной матрицей. Вычислим собственные значения матрицы А:

откуда следует, что Поэтому:

• а) при |А| > 2 И все собственные значения матрицы А вещественны;
• б) при | А | < 2 Н собственные значения матрицы А являются комплексно сопряженными.

В обоих случаях Дет, 0.

Следовательно, система (1.36) устойчива в среднем, если выполнено условие.

III. Устойчивость в среднем квадратическом. Аналогично случаю устойчивости в среднем устойчивость в среднем квадратическом тривиального решения системы (1.38) эквивалентна устойчивости нулевого решения системы ОДУ, в данном случае системы (1.41).

Воспользуемся третьей теоремой. Для этого оценим действительные части собственных значений матрицы А, т. е. корней уравнения.

В общем случае данное уравнение имеет один действительный и два комплексно сопряженных корня. Обозначим их как: то = Xq, Tj = Х + iY, Т‘2 = Х — гУ]. По теореме Виета имеем.

Откуда получаем, что Дет* <�т²/(2ш²). Таким образом, доказана следующая теорема.

Теорема 1.3. Тривиальное решение системы (1.36):

1) устойчиво по вероятности и в среднем квадратическом при

2) устойчиво в среднем при А > 0.

Теперь исследуем вопрос численного решения осциллирующей системы (1.36) для различных режимов устойчивости различными численными методами.

Введем равномерную сетку на временном интервале [fo, X¹] с шагом h > 0:

Обозначим через Х_п+ численное решение в (п + 1)-м узле; - п," — мерный случайный вектор с независимыми между собой стандартными нормальными компонентами.

Проведем сравнение четырех методов по точности вычисления функционалов от решения исследуемого СДУ для различных значений параметров модели, определяющих разные режимы устойчивости. Также исследуем вопрос применимости экстраполяции по Ричардсону для повышения точности оценки функционалов.

Метод Эйлера [200] является часто используемым методом, и поэтому хотелось исследовать возможность решения данной задачи с его помощью. В работе [42] рассмотрен полунеявный метод Эйлера, который является модификацией метода Эйлера для уравнений второго порядка. Он обладает лучшими свойствами устойчивости численного решения, чем метод Эйлера. Обобщенный метод типа Розенброка (с параметром а = ½) является асимптотически несмещенным численным методом решения СДУ с любым шагом интегрирования. Алгоритм переменного шага на основе обобщенных вложенных трехстадийных методов Рунге Кутты Фельберга (параграф 2.2 учебного пособия [31]) был специально построен для решения осциллирующих систем.

Рассмотрим несколько методов, предназначенных для статистического моделирования решений систем СДУ (1.37). Приведем вид этих методов для решаемой задачи.

1. Метод Эйлера [200] для решений систем СДУ в смысле Ито (1.37):

Для упрощения записи будем опускать знак суммирования по индексу, который встречается в слагаемом дважды.

Для СДУ (1.36) метод Эйлера покомпонентно имеет следующий вид:

1а. Метод Мильштейна [112] для систем СДУ в смысле Ито (137) с одним шумом:

Для СДУ (1.36) метод Мильштейна совпадает с методом Эйлера. Поэтому в дальнейшем этот метод рассматриваться не будет.

2. Полунеявный метод Эйлера [42| для решений систем СДУ в смысле Ито (1.37). Это модификация метода Эйлера для систем второго порядка вида (1.36):

Отличие от метода Эйлера состоит в вычислении первой компоненты решения.

• 3. Обобщенный метод типа Розенброка (метод (2) из параграфа 1.1) для решения систем СДУ в смысле Стратоновича (1.37).
• 4. Алгоритм переменного шага (на основе обобщенных вложенных трехстадийных методов Рунге — Кутты — Фельберга) (см. параграф 2.2 из пособия [31]). Этот метод был разработан специально для осциллирующих систем.

Все перечисленные выше методы для произвольной системы СДУ в смысле Ито имеют первый порядок слабой и среднеквадратической сходимости. Для среднеквадратической сходимости это означает, что разложение численного решения Х_п+ в ряд Тейлора в окрестности точки t_n отличается от выражения.

на величину р: (р) = 0(/i^Pl), (р²) = 0(/*²), где р > 1.5 [112].

Для того, чтобы численный метод имел второй порядок среднеквадратической сходимости, достаточно, чтобы разложение численного решения в ряд Тейлора в окрестности точки t_n отличалось от.

на величину р': ((//)) = 0(h²), ((р')²} = 0(h³). Следовательно, для колебательной задачи (1.36) ввиду того, что в разложении решения в ряд Тейлора отсутствует член, содержащий стохастический интеграл У Wj₂ (т) dwj₁ (г), порядок среднеквадратической сходимости методов равен двум.

5. Экстраполяция Ричардсона [144]. В теории обыкновенных дифференциальных уравнений для повышения точности вычислений используется экстраполяция по Ричардсону [144, 212]. Для повышения порядка слабой сходимости численных методов решения СДУ можно также попытаться ее использовать.

Пусть J^h(X_T), J^h/²(X_T) — некоторые функционалы от решения СДУ, полученные численным методом порядка р слабой сходимости с шагами h и h/2 соответственно. Тогда оценка.

имеет порядок р 4- 1 слабой сходимости при /г —" 0.

В работе Д. Талэ (D. Talay) [221] данная процедура называется экстраполяцией Ромберга согласно работе В. Грата (W. В. Gragg) [182] (также см. работу Ромберга [213]).

Зная точное решение функционала J (Xt) и полученные оценки с шагами h и h/2, можно оценить порядок метода по формуле.

В заключение пункта отметим, что на каждом временнбм шаге метод Эйлера и полуиеявный метод Эйлера имеют равное число арифметических операций. Метод типа Розенброка для систем с линейным сносом не требует обращения матрицы на каждом шаге. Поэтому однократное обращение матрицы фактически не увеличивает время вычисления по сравнению с методом Эйлера. Алгоритм переменного шага на каждом шаге требует вычисления вектора сноса в промежуточных точках. Однако это незначительно увеличивает время счета, так как метод позволяет значительно увеличивать шаг.

Все тосты рассматривались на интервале t € [0,50]. Начальные данные Х (0) задавались как гауссовский случайный вектор е параметрами.

Вычисление пары нормальных случайных величин? и у с нулевым математическим ожиданием и единичной дисперсией осуществлялось по формулам.

где «1, о.2 ~ равномерные случайные числа на [0,1].

Тест 1. Тривиальное решение задачи (1.36) устойчиво в среднем и неустойчиво по вероятности и в среднем квадратическом. В этом случае дисперсия будет со временем неограниченно возрастать, в то время как математическое ожидание будет стремиться к нулю. Параметры принимали следующие значения: Л = 0.1, а = 0.5, ш² = 1.

Тест 2. Тривиальное решение задачи (1.36) устойчиво в среднем, по вероятности и в среднем квадратическом. В этом случае со временем дисперсия и математическое ожидание будут стремиться к нулю. Параметры принимали следующие значения: Л = 0.1. гг = 0.1. и)² = 1.

Тест 3. Тривиальное решение задачи (1.36) неустойчиво в среднем, по вероятности и в среднем квадратическом. В этом случае со временем дисперсия и математическое ожидание будут неограниченно возрастать. Параметры принимали следующие значения: А = —0.1, а = 0.1, ш² = 1.

Тест 4. Тривиальное решение задачи (1.36) неустойчиво по вероятности и в среднем квадратическом. В этом случае со временем дисперсия будет неограниченно возрастать, а математическое ожидание характеризуется установившимися периодическими колебаниями. Параметры принимали следующие значения: А = 0, о = 0.1, ш² = 1.

Чтобы численные методы решения ОДУ «отлавливали» амплитуду и фазу численного решения колебательного ОДУ, необходимо брать не меньше 16 узлов на период. На рассматриваемом интервале [0, 50] математическое ожидание первой компоненты решения задачи (1.36) имеет 8 периодов. Следовательно, количество узлов не должно быть меньше 128, таким образом, шаг интегрирования не должен превышать 0.38.

Поэтому при численном моделировании СДУ начальный шаг для исследований был выбран равным 0.1 и моделировалось N = 10⁴ траекторий (см. табл. 1.6, 1.7, рис. 1.19−1.23). Оказалось, что все методы, кроме метода Эйлера, с таким шагом дают реальную картину поведения решения для различных режимов устойчивости.

Использование метода Эйлера с шагом h = 0.1 дает совершенно неверные результаты (см. табл. 1.6). В связи с этим был увеличен объем выборки до N = 10⁶ (см. табл. 1.8−1.11) и проводилось последовательное уменьшение шага, интегрирования в 2 раза с использованием для метода Эйлера экстраполяции, но Ричардсону (см. табл. 1.12).

Для уравнивания статистической и детерминированной погрешностей целесообразно полагать N = 0(h~^2p) (или К^2р C/VN).

Так как все рассматриваемые в статье методы имеют первый порядок слабой сходимости, то для объема выборки N = 10⁶ не рекомендуется брать шаг h < 0.001 (аналогично для N = 10⁴ при h < 0.01 будет преобладать статистическая погрешность).

На рис. 1.19−1.22 приведены графики дисперсий первой компоненты: точные значения (жирная линия) и оценки, полученные численными методами (полунеявный метод Эйлера — точки; метод Розенброка — сплошная линия). Так как метод Эйлера с шагами h = 0.1 и h = 0.05 считает очень плохо, то мы приводим для него оценки, полученные с шагом h = 0.025 (сплошная линия с жирными точками). Для других методов приводятся оценки дисперсии, полученные с шагом h = 0.1. На рисунках не приведены графики для алгоритма переменного шага, так как нас интересовали оценки фзчжционалов в последний момент времени, и внутри интервала мы не требовали точного попадания в узлы сетки. Математическое ожидание имеет вид либо затухающих, либо установившихся, либо «разбачтттывающихся» колебаний. Так как все методы, кроме метода Эйлера, с шагом h = 0.1 во всех тестах оценивают математическое ожидание очень хороню и графики оценок сливаются с точным решением, то мы эти графики не приводим.

Рис. 1.19. Дисперсия первой компоненты. Тест 1.

Чтобы представить поведение математического ожидания решения, мы приводим его только на рис. 1.23 (жирная линия) для теста 4. Здесь также приводятся графики оценок, которые получены методом Эйлера с последовательным делением шага пополам для N = 1()⁴ траекторий (сплошная линия — h = 0.1; точки — h — 0.05; линия с жирными точками -/i = 0.025).

В табл. 1.6, 1.7 представлены точные значения математического ожидания и дисперсии первой компоненты решения задачи (1.36) при t = 50 для тестов 1 4 и их оценки, полученные с помощью рассматриваемых методов для N = 10⁴ траекторий и h = 0.1. Для оценок математического ожидания также приведены доверительные интервалы, полученные с учетом замечания 2. Так как метод Эйлера (1.50) с шагом h = 0.1 очень плохо считает, то в таблицах также приведены оценки, полученные им при h = 0.05 и h = 0.025. Для алгоритма переменного шага расчеты проводились с точностью е = 10″³. Точные значения вычисля;

Рис. 1.20. Дисперсия первой компоненты. Тест 2.

емых функционалов выделены жирным шрифтом.

В табл. 1.8 представлены выходные данные алгоритма переменного шага при моделировании тестов 1−4 на интервале [0,50] для N = 10⁴ траекторий. В табл. 1.9−1.12 можно найти оценки, аналогичные приведенным в табл. 1.6, 1.7, но для большего числа траекторий N = 10⁶. Доверительные интервалы для оценок математического ожидания мы не приводим потому, что их легко получить из доверительных интервалов табл. 1.6, 1.7. Так как оценка дисперсии почти не изменилась, а число траекторий увеличилось в 100 раз, то границы доверительных интервачтов уменьшатся в 10 раз.

Увеличение числа траекторий улучшило оценки, полученные методом Розенброка и полунеявным методом Эйлера, так как уменьшилась статистическая погрешность. У метода Эйлера точность либо не изменилась, либо стала даже хуже (тест 1). Это говорит о том, что для метода Эйлера статистическая погрешность незначительна по сравнению с погрешностью метода. Здесь также приведены результаты численных расчетов для N = 10⁶ траекторий с последовательным делением шага пополам h = 0.1, h = 0.05, h = 0,025, h = 0.0125, h = 0.625,.

Рис. 1.21. Дисперсия первой компоненты. Тест 3.

h = 0.3 125: табл. 1.9 — тест 1, табл. 1.10 — тест 2, табл. 1.11 — тест 3, табл. 1.12 — тест 4. В таблицах жирным шрифтом выделены оценки, полученные методом (1.51) и (2) с шагом h = 0.1, чтобы отметить, что даже с таким большим шагом получены хорошие оценки. У метода Эйлера жирным шрифтом выделены лучшие оценки, чтобы отметить, что даже с мелким шагом h = 1/320 его оценки менее точны.

На основе данных, приведенных в табл. 1.9−1.12, вычислена экстраполяция по Ричардсону для метода Эйлера и приведена в табл. 1.13, 1.14. Точные значения выделены жирным шрифтом.

В данном разделе исследован вопрос устойчивости тривиального решения линейной системы СДУ с мультипликативным шумом в смысле Ито. Получены следующие условия на параметры модели: при Л > 0 решение устойчиво в среднем.

На данной задаче для различных режимов устойчивости проведено сравнение четырех численных методов решения СДУ: метода Эйлера, полунеявного метода Эйлера, метода Розенброка, алгоритма переменного шага на основе обобщенных вложенных трехстадийных методов.

Рис. 1.22. Дисперсия первой компоненты. Тест 4.

Рунге Кутты. Результаты экспериментов позволяют сделать следующие выводы:^[1][2][3]

Рис. 1.23. Математическое ожидание первой компоненты. Тест 5.

4. Экстраполяция по Ричардсону улучшает точность метода Эйлера при мелком шаге (начиная с 1/80). Полученные оценки функционалов с шагами h = 1/160 и h = 1/320 сравнимы с оценками, полученными другими методами с шагом h = 1/10.

Таблица 1.7. Численные результаты (N = 10⁴).

Таблица 1.8. Размер шага в алгоритме переменного шага (N = 10⁴).

Таблица 1.10. Численные результаты. Тест 2 (N = 10⁶).

Таблица 1.12. Численные результаты. Тест 4 {N = 10⁶).

Таблица 1.13. Экстраполяция по Ричардсону для метода Эйлера.

• [1] Для численного решения колебательных задач из рассмотренныхметодов предпочтительно использовать метод Розенброка, алгоритм переменного шага, иолунеявный метод Эйлера. Эти методыдают реальную картину устойчивости исследуемой системы дажепри относительно большом шаге интегрирования h = 0.1.
• [2] Применение метода Эйлера для решения колебательных задач возможно только при очень мелком шаге. При относительно больших шагах метод Эйлера дает неверную картину устойчивостирешения: тривиальное решение задачи (1.36) всегда получаетсянеустойчивым в среднем.
• [3] Для получения численных результатов методом Эйлера, сравнимых по точности с тремя другими рассмотренными методами, необходимо брать шаг интегрирования в 50−100 раз меньший. Ни наодном тесте метод Эйлера с шагом h = 1/320 не достиг точностидругих методов с шагом h = 1/10.