Напомним основные формулы:

, где

x(0) = x₀; .

Установим правило для знака силы:

Гравитация:

Упругие силы (пружина):

Сопротивление воздуха, трение:

Создадим анимацию движения системы «масса – пружина». Для этого будет использовано решение ДУ для гармонического осциллятора y(t), задающего координаты точки, в которой находится движущаяся пружина, как функцию времени. Можно даже уточнить: нам нужно знать положение крайней точки пружины, той, где находится присоединенная масса. Начнем с решения ДУ.

Параметры: масса, постоянная упругости пружины, гравитационное ускорение.

Работаем в единицах СИ: g измеряется в м/с², m в кг, k в Н/м.

Сначала решаем уравнение в общем виде, а затем получаем то, что нам надо:

• with(plots):
• g:='g': m:='m': k:='k':

На массу действуют силы:

• сила тяжести (в форме F_y = –mg), не зависящая от положения на поверхности Земли;
• сила пружины, задаваемая законом Гука, которая противодействует смещению из положения равновесия: F_H = –k (у – у0).

Если y > y0, ускорение действует в отрицательном направлении оси y, и наоборот.

Обратите внимание, что пружины могут быть по-разному сделаны. Обычные пружины плотно намотаны, и их можно только растягивать. Кроме них есть пружины, которые применяются в подвесках автомобилей/мотоциклов/велосипедов: их можно и растягивать, и сжимать. При наличии силы тяжести такая пружина с прикрепленной массой имеет точку равновесия, из которой ее может как сжать, так и (дополнительно) растянуть, т. е. для нее точку равновесия больше нельзя рассматривать как длину пружины. Теперь точка равновесия определяется как точка, в которой гравитационное притяжение к поверхности Земли уравновешивается упругой силой пружины, возникающей в результате ее растяжения от начальной точки (точки равновесия), как противодействие движению массы в гравитационном поле. Упрощение: предполагаем, что пружина может только растягиваться; при полном сжатии (исходное состояние) она имеет нулевую длину. Обозначим высоту, на которой находится система, как y1, а исходное положение, в котором учтен закон Гука, – как y0 = 0. Получаем, что под действием силы тяжести положение равновесия, в котором находится система «масса – пружина»:

• y0:=0; y1:='y1':
• EQL:=-m*g-k*(y1-y0)=0;
• y1:=solve(EQL,y1);

Смещение массы отсчитывается от точки y = 0 (это – точка равновесия без присоединенной массы). Положительное значение y направлено вверх. Тогда сила гравитации F_y = –mg. Уравнение Ньютона:

• NE:=m*diff(y(t),t$2)=-m*g-k*(y(t)-y0);
• sol:=dsolve(NE,y(t));

Понятно, что решение в общем виде получить несложно вручную. Его можно использовать для проверки полученного в Maple результата. Две константы интегрирования определяются начальными условиями. Типичным условием является выход системы из состояния покоя на любой высоте. Можно получить общее решение, затем найти выражение для скорости путем дифференцирования и определить _C1 и _C2. Другой способ – задать начальные условия так, что dsolve предложит особое решение. Пойдем первым путем (произвольная высота и нулевая скорость).

• IC:=y(0)=-0.2,D(y)(0)=0;
• sol:=dsolve({NE,IC},y(t));

Введем константы для графика:

• g:=10; m:=1; k:=10;
• Y_s:=unapply(rhs(sol),t):

Рисуем отрезок длиной 10 с вместе с положением равновесия:

• y1;
• plot([Y_s(t),y1],t=0..10);

Получаем простое гармоническое движение относительно положения равновесия. Чтобы нарисовать анимацию, нужен вид сбоку на систему «масса – пружина». Для этого нужна последовательность отрезков прямых. Подбираем радиус и количество витков пружины:

• w:=0.2: nT:=10:

Масса находится в точке y(t). Делим эту длину на 2nT сегмента и генерируем список точек слева и справа.

• DL:=proc(t) local i,L,d; global w,nT;
• d:=evalf(Y_s(t)/(2*nT+1)); L:=[[0.,0.],[w,d*0.5]]:
• for i from 1 to nT do: L:=[op(L),[-w,d*(2*i-0.5)],[w,d*(2*i+0.5)]]; od:
  L:=[op(L),[0.,evalf(Y_s(t))],[w/4,evalf(Y_s(t))],[w/4,evalf(Y_s(t))-w/8],[-w/4,evalf(Y_s(t))-w/8],[-w/4,evalf(Y_s(t))],[0.,evalf(Y_s(t))]]; end:
• PL:=seq(plot(DL(j_t/10)),j_t=1..32):
• display(PL,insequence=true,axes=boxed);

Есть много способов сделать эту задачу более интересной:

1. введем трение, что сделает задачу реалистичной;
2. рассмотрим свободный гармонический осциллятор.

Начнем с первой задачи. Добавим член, пропорциональный скорости dy/dt, причем добавляемая сила противодействует движению. Тогда добавка должна быть для положительной скорости – отрицательной, для отрицательной скорости (вниз) – положительной. Получаем, что постоянная трения умножается на отрицательную величину скорости:

• b:='b': m:='m': g:='g': k:='k':
• NE:=m*diff(y(t),t$2)=-m*g-k*(y(t)-y0)-b*diff(y(t),t);

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

• sol:=dsolve(NE,y(t));

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

Зададим те же самые начальные условия:

• sol:=dsolve({NE,IC},y(t));

В дополнение к ранее выбранным параметрам (b измеряется в Н · с/м) выберем небольшой коэффициент затухания:

• g:=10; m:=1; k:=10; b:=1/5;
• Y_s:=unapply(rhs(sol),t):
• plot([Y_s(t),y1],t=0..20);

Получаем колебательное решение:

• PL:=seq(plot(DL(j_t/10)),j_t=1..320):
• display(PL,insequence=true,axes=boxed);

Подумаем о реальных приложениях.

Демпфирование в системе подвески автомобиля создают амортизаторы, которые устроены по-разному. На самом деле интересно не конкретное устройство и тип движения (например, который показан выше), а необходимость поглощения энергии колебаний за одно-два колебания. Из аналитического решения с неопределенными константами видно, когда изменяется характер решения. Оно происходит, когда (b² – 4km) = 0 и в экспоненте нет комплексных аргументов (чтобы избежать нулевого знаменателя, приходится добавлять очень малое число):

• g:=10; m:=1; k:=10; b:=sqrt(4*k)+10^(-Digits);
• Y_s:=unapply(rhs(sol),t):
• plot([Y_s(t),y1],t=0..10);

Это называется критическим затухающим решением.

Вопросы: что произойдет, если демпфирование еще больше увеличится? Желательно ли это делать в приложениях?

Удвоим демпфирование:

• g:=10; m:=1; k:=10; b:=2*sqrt(4*k)+10^(-Digits);
• Y_s:=unapply(rhs(sol),t):
• plot([Y_s(t),y1],t=0..10);

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

Есть интересное представление, которое позволяет получить дальнейшее понимание. На так называемой фазово-пространственной диаграмме график зависимости скорости от положения со временем играет роль параметра. Эта кривая (фазовая траектория), позволяет определить, как энергия рассеивается в системе. Для случая без трения получается замкнутая кривая. Решения с различными начальными условиями характеризуются похожими фазовыми траекториями, размер которых является мерой количества полной энергии, сохраняющейся и преобразующейся между формами кинетической и потенциальной энергии. Для случая с демпфированием получаются спирали, поскольку энергия удаляется из системы.

Начнем с подкритического демпфирования:

• g:=10; m:=1; k:=10; b:=1/2;
• Y_s:=unapply(rhs(sol),t):
• V_s:=D(Y_s);
• plot([Y_s(t),V_s(t),t=0..20]);

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

• Etot:=t->1/2*m*V_s(t)^2+m*g*Y_s(t)+1/2*k*Y_s(t)^2;
• plot([Etot(t)-Etot(0),V_s(t)],t=0..5);

Рассмотрите самостоятельно диаграммы для случаев критического и закритического затухания.

Теперь обратимся к случаю управляемого гармонического осциллятора. Выбираем гармоническую движущую силу с малой амплитудой, чтобы вывести гармонический осциллятор из равновесия. Круговая частота ν характеризует движущую силу:

• b:='b': m:='m': g:='g': k:='k': nu:='nu':
• NE:=m*diff(y(t),t$2)=-m*g-k*(y(t)-y0)-b*diff(y(t),t)+(1/10)*sin(nu*t);

Для недемпфированного случая внутренняя частота осциллятора:

• omega:=sqrt(k/m);

Исследуем происходящее, когда ν изменяется в широком диапазоне, и особенно то, что происходит, когда нуль приближается к ω (резонансу).

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

• g:=10; m:=1; k:=10; b:=1/2;
• omega;
• y1;
• nu:=3.1;
• sol:=dsolve({NE,y(0)=y1,D(y)(0)=0},y(t)):

На длинный результат не стоит смотреть:

• Y_s:=unapply(rhs(sol),t):
• plot([Y_s(t),y1],t=0..20);

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

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

• V_s:=D(Y_s):
• plot([Y_s(t),V_s(t),t=0..20]);

Появление предельного цикла говорит нам о том, что наступает момент, когда движущая сила сбрасывает в систему «масса – пружина» столько энергии, сколько энергии теряется из-за трения.

Рассмотрите случаи меньшего трения, когда можно наблюдать движение с большей амплитудой.

Изучите различный набор частот движущей силы.

В конце рассмотрения управляемого осциллятора с демпфированием выясним зависимость конечной амплитуды от частоты движущей силы:

• b:='b': m:='m': g:='g': k:='k': nu:='nu':
• NE:=m*diff(y(t),t$2)=-m*g-k*(y(t)-y0)-b*diff(y(t),t)+(1/10)*sin(nu*t);
• y1;
• sol:=dsolve({NE,y(0)=y1,D(y)(0)=0},y(t)):
• Y_s:=unapply(rhs(sol),t);

Интуитивно понятно, что при подкритическом затухании внутренние колебания затухают в масштабе времени, определяемом постоянной затухания b. В приведенном выше выражении остается первый член, изменение времени которого просто связано с неоднородным членом в уравнении Ньютона, т. е. только с периодической движущей силой.

Предположим, что находимся в асимптотическом режиме, который достигается, когда exp(–bt/2m) пренебрежимо мало по сравнению с первым членом. В решении выбираем первый член. Для b = 1/2 (сильное, но подкритическое демпфирование) должно подойти значение t_a = 20. Из графика видно, что достигнут предельный цикл, т. е. управляемое гармоническое движение с постоянной амплитудой:

• Y_as:=unapply(op(1,Y_s(t)),t);
• plot([subs(g=10,m=1,k=10,b=1/2,nu=31/10,Y_as(t)),
  subs(g=10,m=1,k=10,y1)+1/20*sin(31/10*t)],t=20..40,
  color=[red,blue]);

График показывает, что связанное с предельным циклом гармоническое движение центрировано в точке равновесия y1 и сдвинуто по фазе по сравнению с обычным sin(ν*t). Сдвиг фазы достигается сочетанием синусов и косинусов, присутствующих в решении.

Получим зависимость амплитуды от параметров. Сначала выделим положение равновесия:

• assume(nu,real);
• A:=combine(Y_as(t)-y1,trig);
• A0:=simplify(subs(b=0,A));

В случае без трения асимптотическая амплитуда «взрывается» на резонансе (учтем, что ω² = k/m). Это означает, что предельный цикл не достигается, но амплитуда растет линейно со временем.

Получить этот результат самостоятельно из анализа решения уравнения Ньютона.

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

Вначале посмотрим на асимптотическое решение, в котором ν выбирается на частоте резонанса для гармонического осциллятора без трения (т. е. нет затухания). Для малых констант трения b этот ответ должен быть разумным:

• simplify(algsubs(nu=sqrt(k/m),A));

Пусть вблизи резонанса амплитуда задается по члену с косинусом (в случае слабого затухания два соседних синусных члена вблизи резонанса практически равны нулю):

• A1:=simplify(nu*b*cos(nu*t)/denom(A));

Амплитуда колебаний:

• A1a:=simplify(nu*b/denom(A));

Зададим k и m при резонансе для

• A1as:=subs(m=1,k=10,A1a);
• P0:=plot([sqrt(10),t,t=0..1.5],color=black):
• P1:=plot(subs(b=1/5,A1as),nu=3..3.3,color=red):
• P2:=plot(subs(b=1/15,A1as),nu=3..3.3,color=blue):
• P3:=plot(subs(b=1/45,A1as),nu=3..3.3,color=green):
• display(P0,P1,P2,P3,axes=boxed,labels=[frequency,
  Amplitude]);

График показывает, что демпфированный осциллятор способен поглощать энергию от источника следующим образом:

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

Обсудим колебания системы, расположенной на горизонтальной плоскости цепочки последовательно соединенных пружин с массой. Пружина 1 закреплена с одного конца (скажем, слева), к ее противоположному концу прикреплена масса m1, к ней прикреплена еще одна пружина 2 с массой m2. Внешняя сила действует справа. Рассматривается только горизонтальное движение. Сила тяжести уравновешена движением в воздухе с некоторым трением о поверхность. Коэффициент трения о поверхность равен b. Трение в самих пружинах игнорируется. Упругие константы пружин – k1 и k2, координаты масс (смещения из положений равновесия пружин 1 и 2 соответственно) – x1(t) и x2(t). На массу m1 действует сила Гука от пружины 1, равная –k1*x1(t), и от пружины 2, равная –k2*(x1(t) – x2(t)). Сила для пружины 1 очевидна, а последняя запись проверяется так:

1. если положение равновесия для пружины 2, x2(t), совпадает с положением равновесия пружины 1, т. е. для x1(t) нет силы;
2. если x2(t) больше, чем x1(t), вторая пружина стремится уравновесить эффект первой пружины, а если оно меньше, она будет тянуться в том же направлении, что и пружина 1.

Запишем два уравнения Ньютона для этих пружин с массой.

Для первой пружины – три вклада:

(в обычной записи: ).

Первое слагаемое – движение первой пружины относительно фиксированного конца, второе – движение второй пружины относительно первой, причем при х2 > х1 пружина движется вправо, третье слагаемое – трение первой пружины.

Так как вторая масса с одного конца не закреплена, то для нее распределение сил проще, однако сила, действующая на нее от пружины 2, зависит от обеих координат (аналогично, но противоположно по знаку, как для силы, действующей на массу m1, т. е. при х2 > х1 пружина движется влево).

Для второй пружины с массой:

(в обычной записи: ).

Для пружины 1 в Maple записываем:

• restart; with(plots):
• m1:='m1': m2:='m2': k1:='k1': k2:='k2': b:='b':
• x1(t):='x1(t)': x2(t):='x2(t)':
• NE1:=m1*diff(x1(t),t$2)=-k1*x1(t)-k2*(x1(t)-x2(t))-b*diff(x1(t),t);

Для пружины 2 в Maple записываем:

• NE2:=m2*diff(x2(t),t$2)=-k2*(x2(t)-x1(t))-b*diff(x2(t),t);

Размерности величин возьмем в системе СИ: массы в килограммах, расстояния в метрах, время в секундах и т. п.

Начальные условия зададим исходя из того, что для двух связанных ОДУ второго порядка нужны четыре условия. Вначале считаем, что обе пружины находятся в равновесии и обе массы (m1 и m2) движутся с одинаковой скоростью.

Пусть m1 = m2, и изменяем k1 и k2.

Посмотрим, можем ли мы установить их так, чтобы m2 следовала за массой m1?

В лабораторной системе отсчета положение массы m2 определяется как x1(t) + x2(t).

Проверим наличие ограничений:

1. Пусть k1 >> k2, т. е. пружина 1 намного жестче, чем пружина 2. Движение начинается со сдвига пружины 1 из равновесия. В этом случае уравнение NE1 соответствует гармоническому осциллятору для x1 с небольшим возмущением. Тогда имеет смысл найти приближенное решение для x1(t) без связи с x2(t); в уравнении для массы m2 можно подумать о члене, включающем x1(t) как заранее определенную движущую силу. Движение для m2 в лабораторной системе будет сложным.
2. Изменим в п. 1: пружина 2 запускается вдали от равновесия, в то время как пружина 1 запускается из состояния покоя в положении равновесия x1(0) = 0. В этом случае движение в x1 остается малым; масса m2 ведет себя как возмущенный гармонический осциллятор.
3. Пусть теперь наоборот: k2 >> k1. Независимо от того, начинается ли движение с пружины 1 или с пружины 2, оно соответствует колебаниям около x1(t) = x2(t), результатом чего будут затухающие почти гармонические колебания массы m2 в лабораторной системе координат.
4. Пусть теперь k1 и k2 сопоставимы. Движение осложняется непрерывным обменом энергией между пружинами. Это результат того, что оба генератора находятся в резонансе.

• m1:=1; m2:=1; k1:=4; k2:=5; b:=1/10;
• IC:=x1(0)=0,D(x1)(0)=1/2,x2(0)=0,D(x2)(0)=0;

Это начальное условие подразумевает, что обе массы движутся в унисон при t = 0, т. е. нет относительного смещения и относительной скорости между массами m2 и m1.

• sol:=dsolve({NE1,NE2,IC},{x1(t),x2(t)},method=laplace);
• assign(sol);

Теперь, когда решение назначено, можно записать:

• x1(t);

Но не

• x1(1);
• simplify(%);

Картинки получаются из зависимостей x1(t) и x2(t). Назначение решения имеет недостаток: перед повторным использованием придется повторно назначить зависимые переменные. Именно для этого в начале, до определения уравнения Ньютона, пришлось записать: x1(t): = 'x1(t)'. Это же следует сделать перед повторным использованием назначения с новым выбором параметров (в том же сеансе работы).

• X1:=unapply(Re(x1(t)),t): X2:=unapply(Re(x2(t)),t):
• plot([X1(t),X2(t),X1(t)+X2(t)],t=0..40,color=[red,
  blue,green]);

Из графических решений для обеих масс следует, что:

1. Масса m2 действительно следует за m1. Смысл координаты x2(t) в том, что вторая пружинная масса совершает колебания вокруг положения своего равновесия.
2. Решения подразумевают, что анимация запускается для пружины, которая может смещаться в обе стороны от положения равновесия (т. е. существует предварительное растяжение), но, что более важно, две массы могут переключать свои положения (они могут проходить друг через друга).
3. Для многих начальных условий движения массы m2 в лабораторной системе отсчета – это искаженные затухающие гармонические колебания.

Следующий шаг решения – это создание анимации. Расположим систему горизонтально, чтобы подчеркнуть, что гравитация не учитывается (есть только движение в воздухе с трением о поверхность). Рисуем массы в виде тонких вертикальных плиток, широкие пружины и число их витков:

• w:=0.2: nT:=10:

Для данного положения массы y(t) делим эту длину на 2nT участка и создаем список чередующихся точек слева и справа:

• DL:=proc(t) local i,L,d1,d2,le1,le2,t0,x10; global w,nT; t0:=evalf(t);
• d1:=evalf(X1(t0)/(2*nT+1)); le1:=1; le2:=1; L:=[[0.,0.],[d1*0.5,w]]:
• for i from 1 to nT do: L:=[op(L),[d1*(2*i-0.5),-w],[d1*(2*i+0.5),w]]; od: L:=[op(L),[X1(t0),0.],[X1(t0)-w/80,1.1*w],[X1(t0)+w/80,1.1*w],[X1(t0)+w/80,-1.1*w],[X1(t0)-w/80,-1.1*w],[X1(t0),0.]]; d2:=evalf(X2(t0)/(2*nT+1)); x10:=X1(t0); L:=[op(L),[x10+d1*0.5,w]]: for i from 1 to nT do: L:=[op(L),[x10+d2*(2*i-0.5),-w],[x10+d2*(2*i+0.5),w]]; od: L:=[op(L),[x10+X2(t0),0.],[x10+X2(t0)-w/80,1.1*w],[x10+X2(t0)+w/80,1.1*w],[x10+X2(t0)+w/80,-1.1*w],[x10+X2(t0)-w/80,-1.1*w],[x10+X2(t0),0.]]; end:
• PL:=seq(plot(DL(j_t/10)),j_t=1..320):
• display(PL,insequence=true,axes=boxed);

Интересный график получается при параметрическом построении положения m2 в зависимости от положения m1: если выбрать начальное положение обеих масс в равновесии, график пойдет из начала координат:

• plot([X1(t),X2(t),t=0..10]);

Продолжим расчеты и рассмотрим новую систему, которую легче настроить на воздухе. Три пружины (например, с одинаковой пружинной постоянной k) удерживают две (например, равные) массы m1 и m2 между собой. Пружины предварительно растянуты, так что движение масс может их сжать или растянуть.

Масса m1 прикреплена с одной стороны к пружине 1 (k1), масса m2 – к пружине 2 (k2), а обе массы соединены пружиной k3.

Лучше всего выбрать систему координат, в которой положения равновесия для пружин 1 и 2 представляются одинаковыми координатами с противоположными друг другу знаками. Это приводит уравнения движения в симметричную форму. Сжатие средней пружины измеряется как х1 + х2.

Для такого выбора координат действующая на m1 сила состоит из –k1*x1 и –k3*(x1 + x2). Точно так же для m2 имеем –k2*x2 и –k3*(x1 + x2):

• restart; with(plots):
• m1:='m1': m2:='m2': k1:='k1': k2:='k2': k3:='k3': b:='b':
• x1(t):='x1(t)': x2(t):='x2(t)':
• NE1:=m1*diff(x1(t),t$2)=-k1*x1(t)-k3*(x1(t)+x2(t))-b*diff(x1(t),t);
• NE2:=m2*diff(x2(t),t$2)=-k2*x2(t)-k3*(x1(t)+x2(t))-b*diff(x2(t),t);

Это помогает выбрать слабую связь между двумя генераторами, т. е. слабую пружину между массами m1 и m2:

• m1:=1; m2:=1; k1:=1; k2:=1; k3:=1/10; b:=1/50;
• IC:=x1(0)=0,D(x1)(0)=1,x2(0)=0,D(x2)(0)=0;
• sol:=dsolve({NE1,NE2,IC},{x1(t),x2(t)},method=laplace);
• assign(sol);
• X1:=unapply(Re(x1(t)),t): X2:=unapply(Re(x2(t)),t):
• plot([X1(t),X2(t)],t=0..60,color=[red,blue]);

Если посмотреть на различие в решениях, то окажется, что оно удовлетворяет простому гармоническому движению с затуханием:

• plot(X1(t)-X2(t),t=0..60,color=blue);

Наблюдаются так называемые биения: массы колеблются попеременно. Задав противоположные знаки для координат x1 и x2, получаем очевидное движение массы m2 в направлении, противоположном колебаниям m1.

Теперь надо переписать программу анимации. Если явление биений было легко понять, основываясь на отклонениях масс m1 и m2 от положений равновесия, то для графика движения всей системы понадобится перейти в лабораторную систему отсчета.

Рисунок делаем так: две внешние пружины исключим, а для начальной длины средней пружины выберем произвольное значение (2*Le1):

• w:=0.2: nT:=10:

Для данного положения массы y(t) делим эту длину на 2nT сегмента и генерируем список чередующихся слева и справа точек:

• DL:=proc(t) local i,L,d1,d2,t0,x10,x20,Le1; global w,nT; t0:=evalf(t);
• Le1:=1; d1:=evalf((2*Le1-X1(t0)-X2(t0))/(2*nT+1)); x10:=-Le1+X1(t0); x20:=Le1-X2(t0); L:=[[x10,0.],[x10-w/80,1.1*w],[x10+w/80,1.1*w],[x10+w/80,-1.1*w],[x10-w/80,-1.1*w],[x10-w/80,1.1*w],[x10,0.],[x10+d1*0.5,w]]:
• for i from 1 to nT do: L:=[op(L),[x10+d1*(2*i-0.5),-w],[x10+d1*(2*i+0.5),w]]; od: L:=[op(L),[x20,0.],[x20-w/80,1.1*w],[x20+w/80,1.1*w],[x20+w/80,-1.1*w],[x20-w/80,-1.1*w],[x20,0.]]; end:
• PL:=seq(plot(DL(2*j_t/10)),j_t=0..350):
• display(PL,insequence=true,axes=boxed);

Чтобы картинка выглядела реалистичной, количество фрагментов должно быть велико, тогда конечное решение окажется почти тем же, что было в начале (для достаточно малого трения это достижимо).

Следует попробовать начальные условия, которые соответствуют двум основным модам:

1. Массы движутся совместно (без сжатия пружины, соединяющей m1 и m2). Этого можно достигнуть путем выбора противоположно направленных и равных по величине скоростей и при t = 0 (поскольку x1 и x2 имеют противоположный знак).
2. Массы движутся навстречу друг другу.

Эти две моды соответствуют движениям, описываемым одним уравнением Ньютона в координатах q1 = x1 + x2 и q2 = x1 – x2 соответственно. Можно проверить и убедиться, что два уравнения Ньютона в этих двух случаях расцепляются и что частоты колебаний различны.

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

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

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

• evalf(X1(t));

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

Анализ Фурье можно выполнить аналитически на простых функциях:

• with(inttrans);
• fourier(sin(t),t,f);
• fourier(sin(t)+cos(t/2),t,f);

У компонент косинуса амплитуды действительные, в то время как у синуса – мнимые, а частотный спектр симметричен.

Нужно преобразование, которое интегрируется только при положительных значениях аргумента (время), так как есть проблема с начальными значениями, а демпфирование даст неправильный результат для отрицательных времен. Очевидными кандидатами являются Фурье-косинус и синус. Синусное преобразование Фурье дает очень разумный результат:

• fouriersin(x1(t),t,f);
• plot(%^2,f=0.9..1.2);

Частотный спектр имеет конечную ширину относительно собственных частот.

Изменяя постоянную демпфирования, проверьте, что при слабом демпфировании на собственных частотах спектр приближается к спектру дельта-функции, и наоборот, что спектр расширяется, если увеличивается константа трения b.

Для преобразования косинуса Фурье находим, что квадрат амплитуды обращается в нуль на собственных частотах:

• fouriercos(x2(t),t,f);
• plot(%^2,f=0.9..1.2);

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

Объясните, почему для кажущейся симметричной системы двух одинаковых масс появляются две разные частоты. Брались две одинаковые пружины, к соответствующим сторонам прикреплялись массы, соединенные между собой слабой пружиной. Повторно проанализируйте случаи более сильной и более слабой связи между двумя массами. Очевидно, что для симметричной системы возникли две разные естественные моды: та, где массы не сжимают связывающую их пружину (предположительно ее энергия ниже), и та, где в дополнение к сжатию/растяжению собственных пружин прикрепленными к ним массами возможно сжатие/растяжение сцепляющей их пружины (высокочастотный режим). Именно из-за слабости связи частоты становятся почти вырожденными (равными), и возникают биения. Какое из двух движений имеет более высокую / более низкую частоту, показано приведенным ниже разделением уравнений. Для выбранных параметров преобразуем уравнения Ньютона:

• x1(t):='x1(t)': x2(t):='x2(t)':
• m1:=1; m2:=1; k1:=1; k2:=1; k3:=1/10; b:=1/50;
• NE1:=m1*diff(x1(t),t$2)=-k1*x1(t)-k3*(x1(t)+x2(t))-b*diff(x1(t),t);
• NE2:=m2*diff(x2(t),t$2)=-k2*x2(t)-k3*(x1(t)+x2(t))-b*diff(x2(t),t);
• NE1p:=subs(x1(t)=q1(t)-q2(t),x2(t)=q1(t)+q2(t),NE1);
• NE2p:=subs(x1(t)=q1(t)-q2(t),x2(t)=q1(t)+q2(t),NE2);
• NE1pp:=NE1p+NE2p;
• NE2pp:=NE1p-NE2p;

Развязали два ДУ и теперь ищем их решения:

• sol1:=dsolve(NE1pp,q1(t));
• evalf(rhs(sol1));
• sol2:=dsolve(NE2pp,q2(t));
• evalf(rhs(sol2));

Более низкая частота возникает в режиме, описываемом как q2(t) = 2*(x2(t) – x1(t)). В этом режиме массы движутся параллельно (напомним: в соответствии с соглашением x1 и x2 имеют противоположный знак, а x1/x2 – это мера смещения соответствующей массы от точки равновесия соединяющей пружины).

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

Математическим маятником называют точечную массу m, подвешенную к креплению на невесомой нити длиной l. В декартовой системе координат потребуются два связанных ДУ с ограничением (в обычной записи: ).

Рассмотрим колебания этой системы в поле тяжести. Угол отклонения маятника от вертикали обозначим . Для физического маятника в рассмотрение включают его момент инерции: I = ml². Движение маятника удобно представлять как вращательное. Получаются уравнения Ньютона:

где – угловое ускорение; N – перпендикулярная плоскости движения компонента момента вращения.

Вектор

направлен внутрь плоскости, т. е. N_z < 0. Если увеличивается, то направлены из плоскости (наружу).

Получается:

Гравитационный момент вращения действует как возвращающий.

Для малых амплитуд колебаний sin в правой части можно разложить в ряд Тейлора

и оставить только первое слагаемое.

Получим линеаризованное уравнение движения маятника:

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

.

Нелинейное решение сложнее, его удобнее рассматривать с помощью Maple.

Нам понадобятся все решения, в том числе – случай, когда масса проходит через верхнюю точку. Чтобы представить, как ведут себя разные решения как функция полной энергии, удобно воспользоваться фазовым пространством для переменных .

Кинетическая энергия:

.

Потенциальная энергия:

.

Найдем период математического маятника. Выражения для кинетической и потенциальной энергий записаны выше.

Кинетическая энергия выражается через угловую скорость и момент вращения точечной массы m, подвешенной на нити длиной l:

• restart;
• Tkin:=1/2*m*l^2*omega^2;

Начало отсчета потенциальной энергии выбираем в самой нижней точке (масса находится внизу):

• Upot:=m*g*l*(1-cos(phi));

Для построения графика следует задать значения параметров:

• l:=1;

• g:=10;

• m:=1;

• plots[contourplot](Tkin+Upot,phi=-3/2*Pi..3/2*Pi,omega=-10..10,
  contours=[1,5,9,13,17,21,25,29],axes=boxed,filled=true, coloring=[red,green]);

Показанная на рис. 4.8.1 фазовая диаграмма периодична с периодом 2π. Для малых угловых скоростей система похожа на гармонический осциллятор, а для больших скоростей (или энергий; т. е. для более далеких от центра контуров на графике) есть решения, описывающие простое вращение (кручение) маятника. В периодических фазовых диаграммах, например, каждой точке справа найдется точка слева с той же угловой скоростью. Движение происходит слева направо в верхней половине (для открытых кривых) и справа налево в нижней половине (посмотрите на знак угловой скорости). Единицы энергии (выбранные для контурных линий) указаны в СГС l, g, m.

Рис. 4.8.1

Существует критическая энергия, при которой происходит разделение колебательных решений (замкнутые контуры, т. е. в некотором абстрактном фазовом пространстве можно получить эллиптические или эллипсоподобные «траектории», которыми затем можно представлять колебания, для чего надо следить за знаком ω при прохождении через = 0, что соответствует нижнему положению маятника)

Точки = + π, ω = 0 или = –π, ω = 0 соответствуют неустойчивому равновесию (маятник покоится в верхней точке). Можно показать, что траектория в фазовом пространстве превращается в косинус-кривые для /2.

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

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

Определить траектории в фазовом пространстве было просто. Фактического интегрирования (t) нет. Как оказалось, полная энергия записывается через потенциальную энергию в точке поворота (масса маятника покоится), т. е. через начальный угол 0. Кроме того, вводится энергетический параметр k = sin(0/2). Для записи в общем виде выражения для потенциальной энергии зададим параметры: длина, ускорение свободного падения, масса:

• l:='l': g:='g': m:='m': Upot;

• Upot:=algsubs(cos(phi)=1-2*sin(phi/2)^2,Upot);

• Etot:=subs(phi=phi0,Upot);

Получается, что для данной траектории полная энергия зависит только от k = sin(/2). Следовательно, период рассчитывается как функция k. Чтобы вычислить период колебаний, используем выражение энергии, которое зависит от и ω = diff(, t), выделим ω:

• Tkin+Upot=Etot;

• solve(%,omega);
• solve(%,omega);

• DE:=diff(phi(t),t)=subs(phi=phi(t),%[1]);

• sol:=dsolve({DE,phi(0)=phi0},phi(t));

Для определения периода T можно поступить так: дифференциальное уравнение DE можно разделить так, чтобы получить выражение для dt:

• subs(diff(phi(t),t)=dphi/dt,DE);

• solve(%,dt);

Если правую часть интегрировать по от 0 до 0, левая часть даст нам четверть периода T (для полного периода придется продолжить через 0 до –0 и вернуться к 0). Вычисляемый интеграл оказывается эллиптическим:

• assume(k>0);

Его можно выразить в виде эллиптической функции:

• Int(1/sqrt((1-z^2)*(1-k^2*z^2)),z=0..1);

• value(%);

• omega[0]:=sqrt(g/l);

Удобно считать начальную фазу положительной и меньшей π:

• assume(phi0>0,phi0<Pi);
• T:=2/omega[0]*Int(1/sqrt(sin(phi0/2)^2-sin(phi/2)^2),phi=0..phi0);

• T:=value(T);

Для построения графиков упростим задачу, выбрав длину маятника и гравитационное ускорение в единицах СГС:

• g:=10; l:=10;

• T;

Теперь рассмотрим этот ответ как функцию начального угла 0 (получится не вполне то, что хочется, поэтому придется либо немного уточнить интервал, либо оставить все как есть и перейти к следующему шагу):

• plot(T,phi0=0..Pi);

В пределе начальной фазы, стремящейся к π (неустойчивое равновесие наверху), получается бесконечный период:

• limit(T,phi0=Pi,left);

Для малых начальных фазовых углов (в пределе малой амплитуды мы должны восстановить изохронность) получаем (рис. 4.8.2):

• plot(T,phi0=0..Pi/16);

Рис. 4.8.2

• taylor(T,phi0=0,5);

• T1:=algsubs(sin(phi0/2)=k,T);

• series(T1,k=0,8);

• T1t:=convert(%,polynom);

Это соответствует теории (рис. 4.8.3).

• plot({Re(T1),T1t},k=0..1);

Рис. 4.8.3

Возможен другой способ: ввести переменную k непосредственно в подынтегральную функцию:

• T2:=2/omega[0]*int(1/sqrt(k^2-sin(phi/2)^2),phi=0..2*arcsin(k));

• simplify(T2);

• simplify(T2,symbolic);

Maple неохотно выполняет отмену функции и вычисление обратной ей функции.

Теперь рассчитаем траекторию (t).

Вначале строим фазовую диаграмму (график производной по времени динамической переменной (фазовый угол ) от нее самой). Для этого нужно иметь выражение для энергии (поскольку механическая энергия сохраняется в математическом маятнике без трения, то достаточно нарисовать контурные линии постоянной энергии):

• l:='l': g:='g': m:='m':
• phi:='phi';

• DE;

Анализ:

• omega[0]:=sqrt(g/l);

• Etot;

• E_s:=Etot/(m*l^2*omega[0]^2);

Параметр масштабированной энергии позволяет легко различать различные типы движения:

• для E_s < 2 – колебания;
• для E_s > 2 – вращение;
• E_s = 2 соответствует асимптотическому решению, в котором масса маятника остается в точке неустойчивого равновесия.

Выразим масштабированную энергию через параметр k = sin (0/2):

• k0:=7/8;

• E_s:=2*k0^2;

• phi:='phi':

Записываем ДУ, в котором возможны положительные и отрицательные угловые скорости:

• DE:=(diff(phi(t),t))^2=omega[0]^2*(2*(cos(phi(t))-1+E_s));

Осциллирующие решения уравнения движения (E_s < 2) задаются так:

• theta:=(k,t)->2*arcsin(k*JacobiSN(omega[0]*t,k));

• g:=10; l:=1;

• subs(phi(t)=theta(k0,t),DE);

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

• simplify(%);
• evalf(subs(t=1,%));

Без дальнейшего анализа нельзя понять, является ли выражение для θ(t) решением дифференциального уравнения. Тем не менее можно проверить, что в пределах числовой точности дифференциальное уравнение удовлетворяется.

• plot([theta(0.5,t),theta(0.8,t),theta(0.9,t), theta(0.999,t)],t=0..2*Pi,color=[red,blue,green,black], view=[0..2*Pi,-Pi..Pi]);

Из рис. 4.8.4 видно, что для средних значений энергии решение по-прежнему выглядит как синусоида, хотя зависимость периода от энергии уже очевидна. Только для значений k, близких к единице, видны качественные изменения в поведении решения: в окрестности точки поворота качание маятника резко замедляется.

Рис. 4.8.4

Дополним это анимацией. Чтобы отобразить движение маятника, нужно вычислить координаты x и y. Пусть длина l = 1, и определим координаты так, чтобы высота была нулевой в предполагаемой точке маятника:

• x:=(k,t)->evalf(sin(theta(k,t))): y:=(k,t)->-evalf(cos(theta(k,t))):
• x(0.99,1),y(0.99,1);

Нарисуем линии из начала координат к массе маятника:

• PL:=(k,t)->plot([[0,0],[x(k,t),y(k,t)]],style=line, color=magenta,title=cat(cat("k= ",convert(evalf(k,3),string)), "time= ",convert(evalf(t,3),string))):
• PLs:=seq(PL(0.995,it/10),it=1..200):
• with(plots):
• display(PLs,insequence=true,axes=boxed, scaling=constrained);

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

Рис. 4.8.5

Покажем численное решение. Этот общий подход можно применить к любому дифференциальному уравнению:

• DE:=k->(diff(phi(t),t))^2=omega[0]^2*(2*(cos(phi(t))-1+2*k^2));

• omega[0]:=1;

• IC:=phi(0)=0;

• sol:=dsolve({DE(0.999),IC},phi(t),numeric);

Перепишем заново (если получится сбой в решении):

• omega:='omega': # it is not good to mix omega(t), and omega[0], which makes omega a table.
• DE:=diff(phi(t),t)=omega(t),diff(omega(t),t)=-(g/l)*sin(phi(t));

Начальным условием выбирается точка поворота, где вся энергия – потенциальная:

• phifun:=k->2*arcsin(k);

• phi0:=phifun(0.999);

• IC:=omega(0)=0,phi(0)=phi0;

• sol:=dsolve({DE,IC},[omega(t),phi(t)],numeric, output=listprocedure);

• theta:=subs(sol,phi(t));

В команде plot следует записать θ(t) в кавычках, так как выражение можно оценить, только задав численные значения t (рис. 4.8.6, 4.8.7):

• plot('theta(t)',t=0..2*Pi);

Рис. 4.8.6

• Omega:=subs(sol,omega(t));

• plot('Omega(t)',t=0..2*Pi);

Рис. 4.8.7

Очевидно, что существуют длинные отрезки времени, когда маятник движется с малой скоростью. Полезно посмотреть, как при уменьшении параметра k происходит возврат соотношения между угловой скоростью и углом к случаю гармонического осциллятора. Можно также посмотреть решения, в которых маятник вращается (рис. 4.8.8):

• theta:=(t,k)->2*JacobiAM(sqrt(g/l)*t,k);

• g:=10; l:=1;

• plot(Re(theta(t,0.6*I)),t=0..4);

Рис. 4.8.8

• evalf(theta(1,5.5*I));

Мнимые значения k получены из аналитического продолжения соотношения между энергией и начальным углом:

• E_s:=2.1;

• solve(1-E_s=cos(phi1),phi1);

• k:=evalf(sin(%));

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