Придется почитать (например: Пиппард А. Физика колебаний. М.: Высшая школа, 1985) и посчитать. Возможно, понадобятся знание Maple и С#.

Рассмотрим управляемый гармонический осциллятор с затуханием, пружина которого ослабевает при увеличении амплитуды – в соответствии с уравнением движения:

(3.6.1)

где возвращающая сила F(x)

(3.6.2)

в которое включены и линейное затухание (–2γx), и синусоидальная вынуждающая сила Acos ωt. Масса частицы задана: m = 1.

Нарисуйте F(x) и убедитесь, что ее внешний вид вполне приличный.

Нарисуйте F(x) и потенциальную энергию осциллятора U(x) в интервале x = –3 и x = 3. Эти величины связаны известным соотношением:

(3.6.3)

Постройте график потенциальной энергии вблизи минимума и оцените, насколько он параболичен (если это так, то на частицу в минимуме действует почти такая же сила, как и линейная возвращающая сила в простом гармоническом осцилляторе).

Чтобы в этом удостовериться, разложите F(x) в ряд Тейлора до седьмого порядка и докажите, что для малых смещений х собственная частота колебаний есть ω₀ = 1.

Подсказка: сравнивайте разложение с уравнением простого гармонического осциллятора.

Если частица приобретает большую энергию и поднимается выше, то попадает в область потенциальной энергии, существенно отличающуюся от параболической. В этой части траектории действующая на частицу возвращающая сила меньше и отличается от простой линейной F = –x (см. график F(x)), что означает меньшую частоту колебаний частицы. Такой нелинейный эффект, в результате которого собственная частота становится функцией амплитуды колебаний (такое наблюдается также у маятника – было раньше), имеет странные последствия для резонансной кривой осциллятора.

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

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

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

Представим два разных сценария (рис. 3.6.1).

• (i) Пусть вынуждающая сила относительно мала и мало затухание. Начнем с малой вынуждающей частоты ω и будем ее повышать по направлению к резонансу. Вначале амплитуда мала и совершаются простые гармонические колебания. При приближении к резонансу рост амплитуды сдвигает резонансную частоту вниз. Следует ожидать, что резонанс появится на более низкой частоте ω, чем та, что была рассчитана раньше – ω₀ = 1 (см. нижнюю кривую на рис. 3.6.1).



Рис. 3.6.1. Эскиз кривой нелинейного резонанса

• (ii) Пусть опять вынуждающая сила A мала и мало затухание. Начнем с больших частот и понижаем ω по направлению к резонансу. Вначале амплитуда мала и колебания – гармонические, при приближении к резонансу амплитуда возрастает, резонансная частота понижается, уходя от нас (а не к нам). Продолжаем понижать ω, амплитуда растет. Продолжая так, мы в конце концов окажемся на частоте ниже ω₀ = 1, и даже ниже частоты, при которой был резонанс в п. i (см. верхнюю кривую на рис. 3.6.1). Однако, поскольку понижение резонансной частоты продолжается, амплитуда остается все еще высокой.

Но это означает, что при некоторых значениях вынуждающей частоты возможны две амплитуды колебаний, что показано пунктиром на рис. 3.6.1. Этот эскиз принципиально отличается от обычных резонансных кривых гармонического осциллятора.

Цель работы – рассчитать верхнюю и нижнюю ветви колебаний, изображенных на рис. 3.6.1, и проанализировать условия возникновения таких колебаний в зависимости от параметров задачи.

Следует проанализировать (например, с помощью Maple), как может себя вести резонансная кривая. Анализ – качественный.

Первое приближение: пусть x настолько мало, что в разложении F(x) в ряд Тейлора можно оставить первые два ненулевых члена. Тогда сила в (3.6.2):

Второе приближение: пока обойдемся без затухания в (3.6.1), и тогда уравнение движения примет вид:

(3.6.4)

Мы интересуемся управляемым откликом системы.

Поскольку вынуждающая сила ∝ cosωt, то разумно ожидать, что в отклике x(t) тоже содержатся члены вроде cosωt. Но если x(t) ∝ cosωt, то из-за кубического члена в уравнении движения его решение, возможно, содержит член вроде cos³ωt. A так как cos³ωt = 3/4 cosωt + 1/4cos3ωt, то в x(t) следует ожидать cos3ωt. Эти два слагаемых будем учитывать для приближенного предсказания амплитуды таких установившихся колебаний. Конкретный вид приближения для x(t):

(3.6.5)

Что надо сделать:

(a) Используя Maple, подставьте выражение x(t) (3.6.5) в (3.6.4). Превратите степени косинуса в члены вида cosnωt, n = 1, 3, 5, ..., и затем сгруппируйте члены, пропорциональные cosωt и cos3ωt. Все высшие гармоники игнорируйте в надежде, что они малы.

В результате получится выражение вида . Поскольку cosωt и cos3ωt – независимые функции времени, то их коэффициенты (обозначены как (···)) не одновременно равны 0, что дает уравнение для оценки коэффициентов c₁ и c₃ через A и ω.

(b) Пусть c₁ мало, а c₃ еще меньше. Получим из компоненты уравнения с cos3ωt простое приближенное выражение для c₃ в зависимости от c₁ и ω. При этом надо быть уверенным, что сохраняются только члены, линейные по c₃ (убраны все члены со степенями c₃). При рассмотрении нужно придумать, что делать с членами, содержащими и , сохраняя и убирая , потому что по предположению c₃ << c₁. Получится

(3.6.6)

При каких условиях c₃ не столь мало, как мы предполагаем? В особенности какие ω полностью разрушают это предположение? Это особенное значение ω еще появится, но позже.

(c) Теперь используем это приближение для c₃ в уравнении для cosωt, чтобы получить приближенное кубическое уравнение для c₁ через ω и A, исключая члены высшего порядка по c₁. К удивлению, получится, что c₃ не играет роли в приближенном уравнении для c₁:

(3.6.7)

Решите его, получите три корня при A = 0.136 и подставьте формулы для этих трех значений c₁ в переменные s₁, s₂ и s₃.

(d) Постройте графики s₁, s₂ и s₃:

• (i) постройте |Re(c₁)| при ω = 0..2 (все три корня на одном графике. Строится график модуля |Re(c₁)|, потому что в c₁ нам не надо различать положительные и отрицательные значения);
• (ii) постройте Im(c₁) при ω = 0..2 (все три корня на одном графике).

Вторая часть важна, потому что если c₁ имеет комплексную часть, то этот корень не имеет физического смысла и соответствующую часть кривой из графика Re части в п. (i) не надо учитывать.

(f) Теперь известно, какие корни физичны, они видны на графике |c₁|, и отметим, что для вынуждающей частоты ω ниже 0.6 есть три действительных корня c₁, но выше 0.6 есть только один.

Эффект, при котором из трех решений получается одно, делает такой осциллятор интересным и подсказывает, как дополнить рис. 3.6.1. Чтобы узнать причину, рассмотрим три сценария, начав с некоторой вынуждающей частоты ω и медленно изменяя ее с помощью уже полученных кривых c₁.

• (i) Начнем с большой вынуждающей частоты ω и станем медленно понижать ее. Выше 0.6 есть единственная амплитуда c₁, которая должна повышаться и при понижении частоты (так предсказывает график). При достижении 0.6 появляются три возможных решения, но если мы меняем частоту медленно, то, возможно, останемся на верхней ветке и продолжится рост c₁.
• (ii) Пусть стартуем с низкой частоты и медленно повышаем ее до 0.6. Амплитуда c₁ будет медленно расти вдоль нижней ветки, но вскорости она начнет расти быстро по причине эффекта «раннего резонанса», обсужденного в сценарии (i). При переходе через 0.6 нижняя ветка прекращает существование, и осциллятор вынужден перепрыгнуть в верхнюю ветку, после чего амплитуда уменьшается с ростом частоты ω.
• (iii) Теперь о промежуточной ветке амплитудной зависимости для частот между ω = 0 и ω = 0.6. Впоследствии окажется, что эта ветка нестабильна, и осциллятор с этой ветки быстро перескочит либо на нижнюю, либо на верхнюю. При добавлении затухания окажется, что две верхние ветви объединяются, а не продолжаются до ω = 0, как на ранее полученных кривых. Неожиданностью является то, что кривые нелинейного резонанса похожи на кривые обычного гармонического осциллятора, которые выпячиваются, как шляпа ведьмы. Для получения таких кривых придется написать код Maple.

Важным нелинейным эффектом является то, что для частот ниже 0.6 нельзя указать, какой будет вынужденная амплитуда осциллятора, без знания деталей изменения вынуждающей частоты. Этот эффект назван гистерезисом (в переводе с греческого – отставание, запаздывание). Такой термин уместен, поскольку состояние осциллятора не определяется однозначно текущим значением ω, но и тем, какой она была: система помнит, по какому пути изменялась ω.

Расчеты (Maple)

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

Прямой путь получить отклик осциллятора – задать начальные условия x(0), v(0) и запустить расчет до получения устойчивого состояния. К сожалению, это долго и не позволяет учесть нестабильную среднюю ветвь резонансной кривой.

Лучше сделать так: учесть, что отклик осциллятора на вынуждающую силу устойчив в том, что x(t) и v(t) обе должны быть периодичны с одинаковым периодом, таким, как у вынуждающей силы T = 2π/ω. Для каждой вынуждающей частоты ω должен быть набор начальных условий x(0), v(0), создающий конечное устойчивое состояние.

Искать это особое состояние можно, задав начальные значения x(0) и v(0), решив затем ДУ в области от 0 до T и проверив, соответствуют ли x(T) и v(T) их начальным значениям. Если нет – то продолжаем подгонку x(0), v(0), пока не совпадут. Если повезет, решение ДУ по короткому интервалу [0, T] окажется тем, что мы ищем, а это ускорит расчеты по сравнению с обычным поиском устойчивого состояния. Но чтобы быть удачливым, надо по-умному искать начальные условия, для чего запустим задачу об оптимизации (поиск минимума).

Понадобится пошаговая инструкция для действий.

Попытаемся написать код для получения рисунка вроде рис. 3.6.2. Но это будут подсказки. Остальное делайте сами.

Рис. 3.6.2. Резонансная кривая нелинейного осциллятора

Подсказки

Вам придется использовать процедуры и циклы в Maple:

1. Чтобы рассчитать отклик осциллятора при прохождении частоты ω через резонанс, нужен цикл, который изменяет значения ω с ω₁ до ω₂ за N шагов:

(a) глобальные переменные: частота ω, вынуждающая амплитуда A, коэффициент затухания и S (ее назначение станет понятно позже). Их начальные значения: A = 0.136 и γ = 0.2;

(b) задается диапазон частот ω₁, ω₂ и число точек в нем N. Расчет шага h = (ω₂ – ω₁)/N, задание wscan-массива величин ω.

2. Организация цикла по значениям ω и использование команды поиска минимума функции для того, чтобы, стартовав с начальных x(0) и v(0), перейти к периодическим x(t) и v(t) с периодом T = 2π/ω – таким же, как период вынуждающей силы. Удобно поиск минимума применить к скалярной функции .

Возможно, вы предложите свой вариант решения.

Для отладки используйте значения A = 0.136, γ = 0.2, ω₁ = 1.8, ω₂ = 2, N = 5 (т. е. расчеты будут выполняться намного выше резонанса ω = 1). При таких частотах амплитуда должна быть мала и x(t) должно быть на 180° смещено по фазе от вынуждающей силы, поэтому мы ожидаем что-то вроде x(0) = –A/ω₂ и v(0) = 0 для получения периодического решения. Добейтесь от кода выдачи приемлемого результата. При численной проверке с γ = 0.2 и ω = 2 правильное начальное значение есть (x₀, v₀) = (–0.0423, 0.0225).

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

Важно: везде в вычислениях используйте A = 0.136.

(a) Задайте γ = 0.015 и сделайте скан по частотам от ω₁ = .1 до ω₂ = .4 при N = 80. Теперь окажется очень плохо, если частота маленькая, а x(0) = –A/ω². Возьмите вместо этого x(0) = 0.15. Значение N = 80 должно быть достаточно велико, чтобы можно было видеть предсказанное резонансное поведение при ω = 1/3, а также другой малый резонанс при ω = 1/5. Он получается из-за того, что при анализе мы не учитывали высшие гармоники, а этот вклад – от cos5ωt (кстати, вкладами от cos7ωt, cos9ωt и т. д. тоже пренебрегали). Можете ли увидеть что-то на ω = 1/7? А на ω = 1/9? Почему их труднее увидеть и почему резонансы немного сдвинуты от этих значений? (Увеличьте графики, чтобы увидеть.)

(b) Теперь возьмите ω₁ = 0.2, ω₂ = 2, γ = 0.2 и N = 40. Увидите, что при таком затухании гистерезиса нет, хотя намечается легкий изгиб на резонансных кривых, подсказывающий, что может произойти.

Теперь – главное событие.

Задайте γ = 0.15 и прокрутите цикл по пяти разным наборам ω₁ и ω₂.

(a) Use ω₁ = 2, ω₂ = 0.7, чтобы отследить верхнюю кривую при уменьшении частоты от низких амплитуд при высокой частоте. Задайте N = 30. Удостоверьтесь, что на экран выдаются ω, (x₀, v₀) и S, – чтобы их можно было отследить. Результаты (массивы wscan и Amp – частота и амплитуда) запоминайте в файле, чтобы впоследствии можно было построить свою версию рис. 3.6.2. При поиске минимума учитывайте особенности применения команды поиска минимума, внимательно читая справку.

Лучше всего, конечно, угадывать хорошие стартовые значения (x0, v0).

(b) Для завершения верхней кривой возьмите ω₁ = 0.7, ω₂ = 0.55 и N = 30. Для старта найдите (x₀, v₀) из первого скана, где было ω = 0.7. Запоминайте массивы wscan и Amp (с другим именем).

(c) Для нижней кривой задайте ω₁ = 0.2 и ω₂ = 0.7, чтобы отследить ее до точки, где она сворачивает и прыгает на верхнюю ветвь. N = 30. Всегда ли скачок вверх получается вблизи значения, предсказанного графиком, сделанным ранее? Случается ли скачок при ω больше, чем последнее хорошее значение на верхней кривой в (ii)? Если да, то вы достигли гистерезиса. Поздравляю. Запоминайте массивы значений частоты и амплитуды (с другим именем).

(d) Для средней ветви задайте ω₁ = 0.6, ω₂ = 0.55, (x₀, v₀) = (0.35, 0.31), N = 20. Получите левую половину средней ветви.

(e) Теперь задайте ω₁ = 0.6, ω₂ = 0.7, (x₀, v₀) = (0.35, 0.31), N = 20. Получите правую часть средней ветви, и это последний скан.

(f) Создайте полную картину резонанса, для чего объедините в одном графике все пять сканов. Возможно, для этого понадобится еще один скрипт.

(g) Модифицируйте код, чтобы поисследовать стабильность колебательных состояний, найденных ранее. Например, замените интервал [0, T] на [0, 100*T], чтобы увидеть поведение в большом интервале времен.

• (i) Начните с ω₁ = 0.3, а ω₂ равно чему-то другому, N = 2. Не надо заботиться об ω₂, так как нам нужна только первая точка, ω = ω₁. Получившиеся графики x(t) и v(t) проанализируйте на предмет стабильности (или не-) колебательного состояния. Перейдите на следующее значение.
• (ii) Теперь начните скан с ω₁ = 1, посмотрите на (не)стабильность состояния на верхней кривой.
• (iii) В конце ω₁ = 0.6, (x₀, v₀) = (0.35, 0.31), и решите – стабильно или нет это состояние. Вы должны найти, что оно нестабильно, но нестабильность насыщается установкой в состояние либо верхней, либо нижней кривой, каждая из которых стабильна. При перескоке на верхнюю кривую измените начальные условия в строке скрипта, задав ответ, получаемый при решении ДУ с коэффициентом 0.999999, и запустите снова. При прыжке вниз измените коэффициент на 1.000001. Получите, что направление прыжка – вверх или вниз – очень чувствительно к начальным условиям. Поэтому при экспериментальном исследовании подобных систем мы не можем точно знать, где средняя ветвь.

Итак, вы нашли кривую, подобную изображенной на рис. 3.6.2, и определили, что верхняя и нижняя ветви стабильны, а средняя нестабильна. При малом затухании есть малые подрезонансы на ω =ω₀/3, ω = ω₀/5 и т. д.