CC BY
53
Вестник СамГУ. 2015. № 6(128) 141
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 621.373.12+519.622
В.В. Зайцев, А.В. Карлов, Ар.В. Карлов1
О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТОМСОНОВСКИХ АВТОКОЛЕБАТЕЛЬНЫХ СИСТЕМ
Предложен алгоритм численного интегрирования задачи Коши для уравнений движения автоколебательных систем томсоновского типа. Алгоритм основан на использовании отсчетов импульсной характеристики линейной резонансной системы в качестве дискретизирующей последовательности при переходе к дискретному времени в интегральной форме уравнений движения. Даны оценки погрешности численных решений. Обсуждается трансформация конечно-разностного вычислительного алгоритма в объект нелинейной динамики в дискретном времени. Предложен вариант дискретного отображения осциллятора Ван дер Поля.
Ключевые слова: автоколебательная система, интегральное уравнение Вольтерра, импульсная характеристика резонатора, конечно-разностный алгоритм, нелинейная динамика в дискретном времени, дискретное отображение осциллятора Ван дер Поля.
1. Предварительные сведения
В соответствии с общепринятой классификацией к автоколебательным системам (АКС) томсоновского типа относятся системы с малым энергообменом за период автоколебаний [1]. Такие системы, как правило, содержат высокодобротный колебательный контур и слабо нелинейный активный элемент. В радиотехнических устройствах они применяются в качестве источников квазигармонических колебаний.
Приближенный анализ автоколебаний в томсоновских системах базируется на различных вариантах асимптотических методов теории нелинейных колебаний [2; 3]. Для численного моделирования вполне применимы стандартные методы интегрирования задачи Коши для систем обыкновенных дифференциальных урав-
х© Зайцев В.В., Карлов А.В., Карлов Ар.В., 2015
Зайцев Валерий Васильевич (zaitsev@ssu.samara.ru), кафедра радиофизики и полупроводниковой микро- и наноэлектроники, Самарский государственный университет, 443011, Российская Федерация, г. Самара, ул. Акад. Павлова, 1.
Карлов Алексей Владимирович (aleksey.karlov@gmail.com), ОАО "РКЦ "Прогресс", 443009, Российская Федерация, г. Самара, ул. Земеца, 18.
Карлов Артем Владимирович (Ar.Karlov@gmail.com), кафедра радиофизики и полупроводниковой микро- и наноэлектроники, Самарский государственный университет, 443011, Российская Федерация, г. Самара, ул. Акад. Павлова, 1.
нений [4]. Рассматриваются также численные решения интегральных уравнений Вольтерра второго рода [5].
В настоящей работе предлагается метод численного моделирования томсонов-сих АКС, основанный на использовании осциллирующего отклика их резонаторов на импульсное воздействие.
2. Осцилляторный метод моделирования
В качестве базовой дифференциальной модели томсоновской АКС рассмотрим уравнение
d?u "0 du , 2 с ( du\
W + Qdt + " u = ^[U'A),
^ + ~ + "lu = "lf u^ , (2.1)
где ш0 и Q — собственная частота и добротность резонансного контура, /(о) — нелинейная функция, описывающая кольцо обратной связи, содержащей активный элемент.
Альтернативной (2.1) формой описания АКС является ее интегральная модель в виде уравнения Вольтерра второго рода [5]
i(t)= [ f (u(r),u(r)) h(t - r)dr + U(t). (2.2)
J 0
г* /0
Здесь функция времени и(I) описывает свободные колебания в резонаторе.
Ядро уравнения (2.2) Н(1 — т) —импульсная характеристика резонатора, удовлетворяющая уравнению
d?h "0 dh п, , di+Q -ж+"0h="ls(t) (2.3)
с дельта-функцией в правой части и нулевыми начальными условиями. Отметим, что свободные колебания U(t) также подчиняются уравнению (2.3), но зависят от начальных условий. Решением (2.3) при t ^ 0 является функция
, . . "l f"0\ h(t) = — exp -— t sm("1t). "1 V 2Q J
Частота осцилляций "i в высокодобротном контуре с точностью до величин второго порядка малости по параметру ц = 1/Q совпадает с "0:
Ш1 = Ш0\/1 — = ш0 + 0(р2).
Потому в дальнейшем, считая эти частоты одинаковыми, будем использовать импульсную характеристику
Н(г) = Ш0 ехр ^ — ^ вт(иог). (2.4)
Легко установить, что отсчеты характеристики (2.4) на равномерной сетке Ьп = = п&Ь при п ^ 2 удовлетворяют разностному уравнению
Г>2 [Н(гп)} = Н(гп) — 2а ссв(ш0^г)Н(гп-1) + а2Н(Ьп-2) = 0, (2.5)
где
Ч - i 4
и начальным условиям h(to) = h(0) = 0, h(ti) = h(At) = ашо sin(woAt).
Последовательность отсчетов h(tn) применим при дискретизации времени в интегральном уравнении (2.2). Для этого создадим дискретизирующую последовательность
hd(t) = AtJ2h(tm)S(t — tm). (2.6)
m=0
Импульсную характеристику h(t — т) в (2.2) заменим последовательностью (2.6). После интегрирования для значений решения на временной сетке tn получим уравнение дискретной свертки
n
u(tn) = At^2 f (u(tn — tm), u(tn — tm)) h(tm) + U(tn). m=0
Это уравнение имеет вторую форму
n
u(tn) = AtJ2f (u(tk), ü(tk)) h(tn — tk) + U(tn). (2.7)
k=o
Теперь обе части уравнения (2.7) подвергнем воздействию разностного оператора D2{o} (см. (2.5)). При n ^ 2 это приводит к следующему результату:
D2{u(tn)} = Ath(At)f (u(tn — At), U(tn — At)). (2.8)
Отметим, что в процессе преобразования правой части (2.8) использованы равенство h(0) =0, уравнение (2.5) и аналогичное ему уравнение D2{U(tn)} =0.
Следует обратить внимание на то, что равенство h(0) = 0 приводит к зависимости правой части (2.8) от момента времени tn — At = tn—i. А это, в свою очередь, позволяет получить явный алгоритм моделирования автоколебаний, если аппроксимировать производную u(tn—i) левыми разностями. Ограничимся трехточечной аппроксимацией второго порядка точности:
3un-i —4un-2+u
n — 3 /f-ч r.\
un-1 =-2At-. (2J)
Более высокий порядок точности представления производной использовать не следует, т. к. левая часть равенства (2.8), точная для квазигармонических колебаний, для колебаний, содержащих гармоники основной частоты, не может аппроксимировать левую часть исходного уравнения (2.1) с порядком выше второго.
Таким образом, алгоритм моделирования автоколебаний (нелинейных колебаний) второго порядка определяется рекуррентной формулой
un — 2а cos(^o At)un—i + a2un—2 = Ath(At)f (un—i ,un— 1), (2.10)
где n = 3, 4,..., а для значения ui и u2 в соответствии с интегральным уравнением (2.2) определяются по начальным условиям uo, uо с помощью выражений
ui = 1Ath(At)f (uo, uо) + Ui,
u2 = Ath(2At)f (uo, uo) + U2. Физическая интерпретация алгоритма (2.10) состоит в том, что он является уравнением движения дискретного осциллятора, возбуждаемого сигналом кольца обратной связи, содержащей нелинейный активный элемент. Учитывая это, в дальнейшем предложенный метод условно будем называть осцилляторным.
а
3. Пример результатов моделирования
В качестве тестовой рассмотрим задачу о расчете процесса установления автоколебаний в осцилляторе Ван дер Поля с равнением движения вида
З2и ш° ¿и 2 , ¿и
+ + и = ^ - (31)
где 7 — пареметр глубины положительной обратной связи в системе. Результаты расчетов по алгоритму (2.10) будем сравнивать с расчетами методами Рун-ге — Кутты 2-го (методом Хойна) и 4-го порядков. Последний имеет заведомо большую точность, чем (2.10), и его результаты рассматриваются в качестве эталонных. Расчеты проведены для осциллятора с параметрами Q = 10 и 7 = 0,3 (превышение порога генерации р = YQ = 3). Графики, представленные на рис. 1, получены интегрированием задачи Коши для уравнения (3.1) с начальными условиями и(0) = 0,13, и(0) =0 с шагом ш°Д; = 0, 586.
На рис. 1, а точками показаны отсчеты и(;п), вычисленные по методу (2.10), а ромбами — по стандартному методу Рунге — Кутты 4-го порядка. Временные зависимости амплитуд автоколебаний прведены на рис. 1, б: непрерывная линия — метод Рунге—Кутты 4-го порядка, пунктирная линия — осцилляторный метод (2.10), штрихпунктирная линия — метод Хойна. Амплитуды оценивались путем детектирования и низкочастотной фильтрации временного ряда:
Л(1п) = ЕБ{\и(гп)\],
где ГБ{о} —оператор дискретного фильтра Баттерворта четвертого порядка. При анализе разультатов фильтрации следует принимать во внимание динамические эффекты фильтра, в частности, наличие групповой задержки.
Амплитудные спектры автоколебаний, вычисленные по их реализациям, представлены на рис. 1, в: перерывная линия — метод Рунге — Кутты 4-го порядка, точки —метод (2.10). Частота анализа ш на графиках нормирована на частоту дискретизации: О = шД;/2п.
Сопоставление результатов, полученных различными методами, позволяет сделать вывод о том, что на практике нашел подтверждение заявленный второй порядок точности предложенного осцилляторного метода (2.10). При этом осцилля-торный метод имеет меньшую погрешность, чем формально эквивалентный ему метод Рунге — Кутты второго порядка. С ростом шага интегрирования (параметра ш°Д;) решение осцилляторным методом дает завышенные частоту и уровень гармоник автоколебаний, чем метод Рунге — Кутты четвертого порядка. В конечном итоге это приводит к неустойчивости решения при шагах интегрирования, не удовлетворяющих условию ш°Д; ^ п/3. Следует, однако, отметить, что в области п/3 ^ ш°Д; ^ п неустойчивость численного решения уравнения (3.1) осцилля-торным методом не всегда сопропождается неограниченным ростом членов временного ряда и(;п). Существуют локальные области параметра ш°Д;, в которых алгоритм (2.10) становится генератором хаотического временного ряда. В этом случае (2.10) следует рассматривать как самостоятельный объект нелинейной динамики дискретного времени — дискретный осциллятор (ДВ-осциллятор) Ван дер Поля (см. также [6]).
0.01
Рис. 1. Характеристики автоколебаний: шаг интегрирования и>оА1 = 0, 586
шпг
4. Трансформация вычислительного алгоритма в ДВ-систему
В динамике дискретных во времени нелинейных систем (нелинейных ДВ-систем) численный алгоритм вида (2.10) становится самостоятельным объектом исследований. Он задает динамическую систему, свойства которой могут существенно отличаться от свойств аналогового прототипа, моделируемого с помощью (2.10). Подробнее рассмотрим это на примере осциллятора (3.1).
Естественным временным масштабом в ДВ-системах является интервал дискретизации ДЬ. С ним связан частотный масштаб — частота дискретизации О = = 2п/ ДЬ, в единицах которой измеряются все частоты ДВ-системы. В частности, собственная частота дискретного резонатора О0 = ш0/О^, а параметр диссипации а определяется выражением
а = ехр ^—п
С учетом введенных обозначений ДВ-осциллятор Ван дер Поля задается разностным уравнением движения
ип — 2а ео8(2п0о)м„_1 + а2ип-2 = 7(1 — «П-1)(3ип-1 — 4м„-2 + «п-з) (4.1)
с начальными условиями ио,«1,«2. Уравнение содержит три назависимых параметра: Оо , 7 и Q, впрочем, вместо добротности прототипа Q независимым можно считать параметр диссипации а. При различных комбинациях их значений реализуются как регулярные, так и хаотические режимы автоколебаний.
Для примера на рис. 2 приведены усредненные амплитудные спектры автоколебаний осциллятора (4.1) с параметрами Оо = 0, 32, Q = 20. Пунктирная линия соответствует значению 7 = 0,1, сплошная — 7 = 0,147. Спектральные оценки проведены методом Бартлетта с 512-точечным дискретным преобразованием Фурье по отрезку реализации ип длины N = 65 536. Широкий спектр автоколебаний является одним из эвристических признаков [7] генерации динамического хаоса при 7 = 0,147.
Рис. 2. Амплитудные спектры ДВ-осциллятора
В целом спектр хаотических автоколебаний состоит из уширенной линии основной частоты и хаотического пьедестала, включающего уширенные спектральные линии подменных гармоник. Анализ методами теории аналитического сигнала динамики огибающей и фазы хаотических автоколебаний позволил выявить механизмы уширения линий и формирования пьедестала. Далее приведены некоторые результаты анализа. Отметим, что они полностью соответствуют классической теории формы спектральной линии автоколебаний [8].
Для реализации ДВ-автоколебаний выбрано представление вида
ип = Ап соъ(фп) = Ап СОв(2п0.аП + фп),
где огибающая Ап и фаза фп —хаотические процессы, а частота 0.а вычисляется путем усреднения полной фазы фп по всей длине реализации. В хаотическом спектре, показанном на рис. 2, средняя частота Па = 0,366.
Типичный график хаотической огибающей представлен на рис. 3.
Рис. 3. Хаотическая огибающая автоколебаний
То, что такой хаотический процесс не приводит к уширению спектральной линии автоколебаний, а лишь к появлению амплитуды пьедестала, подтверждается спектром автоколебаний
Уп = Ап сов(2пИап),
показанным на рис. 4 линией У.
Процесс фп представляет собой хаотическую диффузию фазы автоколебаний. Реализации процесса для различных наборов начальных значений ио,и\,П2 отображены на рис. 5 линиями 1-4.
Низкочастотная составляющая хаотической диффузии приводит к уширению спектральной линии, а высокочастотная — к образованию фазового пьедестала, как это показано на рис. 4 пунктирной линией и. Она соответствует спектру автоколебаний
ип = Асов(2пИап + фп)
с усредненной амплитудой (для рассматриваемой реализации А = 1, 75) и диффундирующей фазой.
Рис. 4. Компоненты амплитудного спектра хаотических автоколебаний
Рис. 5. Диффузия фазы хаотических автоколебаний
Таким образом, основной механизм хаотизации автоколебаний ДВ-осциллятора (4.1) — это хаотическая диффузия фазы, приводящая к уширению спектральной линии на основной частоте и ее гармониках.
Заключение
Описанный выше вычислительный метод предлагается использовать в численных экспериментах по исследованию временных и частотных характеристик нелинейных резонансных систем. По сравнению с близким ему по точности методом Рунге — Кутты второго порядка он имеет большую вычислительную эффективность за счет однократного, а не двукратного вычисления нелинейных функций в уравнениях движения колебательных систем. Кроме того, метод позволяет прово-
дить синтез нелинейных динамических систем, функционирующих в дискретном времени.
Литература
[1] Основы теории колебаний / под ред. В.В Мигулина. М.: Наука, 1978. 392 с.
[2] Боголюбов Н.Н., Митропольский Ю.А. Асимптотические методы теории нелинейных колебаний. М.: Наука, 1974. 504 с.
[3] Найфе А. Введение в методы возмущений. М.: Мир, 1984. 536 с.
[4] Хайер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990. 512 с.
[5] Зайцев В.В., Зайцев О.В., Никулин В.В. Интегральные модели автоколебательных систем // Физика волновых процессов и радиотехнические системы. 2006. Т. 9. № 1. С. 53-57.
[6] ДВ-осцилляторы, порождаемые томсоновскими автоколебательными системами / В.В. Зайцев [и др.] // Физика волновых процессов и радиотехнические системы. 2008. Т. 11. № 4. С. 98-103.
[7] Многоликий хаос / Е.Ф. Мищенко [и др.]. М.: Физматлит, 2013. 432 с.
[8] Малахов А.Н. Флуктуации в автоколебательных системах. М.: Наука, 1968. 660 с.
References
[1] Fundamentals of vibration theory. V.V. Migulin (Ed.). M., Nauka, 1978, 392 p. [in Russian].
[2] Bogolyubov N.N., Mitropolsky Yu.A. Asymptotic methods in the theory of non-linear oscillations. M., Nauka, 1974, 504 p. [in Russian].
[3] Nayfe A. Introduction to perturbation methods. M., Mir, 1984, 536 p. [in Russian].
[4] Hairer E., Nersett S., Wanner G. Solving of ordinary differential equations. Nonrigid problems. M., Mir, 1990, 512 p. [in Russian].
[5] Zaitsev V.V., Zaitsev O.V., Nikulin V.V. Integral models of self-oscillating systems. Fizika volnovykh protsessov i radiotekhnicheskie sistemy [Physics of wave processes and radio systems], 2006, Vol. 9, no. 1, pp. 53-57 [in Russian].
[6] Zaitsev V.V., Zaitsev O.V., Karlov A.V., Karlov A.V.(junior). Discrete-time oscillators produced by Thomson autooscillator systems. Fizika volnovykh protsessov i radiotekhnicheskie sistemy [Physics of wave processes and radio systems], 2008, Vol. 11, no. 4, pp. 98-103 [in Russian].
[7] Mishchenko E.F. [et al.] Many face of chaos. M., Fizmatlit, 2013, 432 p. [in Russian].
[8] Malakhov A.N. Fluctuations in self-oscillating systems. M., Nauka, 1968, 660 p. [in Russian].
V.V. Zaitsev, A.V. Karlov, Ar.V. Karlov2
ABOUT NUMERICAL MODELLING OF THOMSON SELF-OSCILLATORY SYSTEMS
The algorithm of numerical integration of a task of Cauchy for the equations of the movement of self-oscillatory systems of Thomson type is offered. The algorithm is based on the use of samples of impulse response of linear resonant system as discretization sequences at the transition to the discrete time in the integral form of the equations of motion. Estimates of an error of numerical decisions are given. Transformation of finite difference algorithm in object of nonlinear dynamics in discrete time is discussed. Version of discrete mapping of Van der Pol oscillator is proposed.
Key words: self-oscillatory system, Volterra integral equation, impulse response of resonator, finite difference algorithm, nonlinear dynamics in discrete time, discrete mapping of Van der Pol oscillator.
Статья поступила в редакцию 28/V/2015. The article received 28/V/2015.
2 Zaitsev Valeriy Vasilievich (zaitsev@ssu.samara.ru), Department of Radiophysics and
Semiconductor Micro- and Nanoelectronics, Samara State University, 1, Acad. Pavlov Street, Samara, 443011, Russian Federation.
Karlov Alexey Vladimirovich (aleksey.karlov@gmail.com), Joint Stock Company Space Rocket Centre Progress, 18, Zemetsa Street, Samara, 443009, Russian Federation.
Karlov Artem Vladimirovich (Ar.Karlov@gmail.com), Department of Radiophysics and
Semiconductor Micro- and Nanoelectronics, Samara State University, 1, Acad. Pavlov Street, Samara, 443011, Russian Federation.
