Дифференциальные уравнения

Система тел и пружин

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

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

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

Данная система схематически изображена на рисунке 1. Она состоит из двух тел, массы которых равны m₁ и m₂, и трех пружин с коэффициентами жесткости k₁, k₂, k₃. Смещение тел от положений равновесия определяется координатами x₁ и x₂.


Рис.1		Рис.2

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

кинетическая и потенциальная энергии тел с пружинами

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

Составим дифференциальные уравнения Лагранжа:

Найдем частные производные:

В результате получаем следующую систему уравнений, описывающую движение тел:

Эту систему можно представить в матричной форме:

Будем искать решение X(t) в форме гармонических колебаний, т.е. в виде

где B₁, B₂ − амплитуды колебаний тел, ω − собственные частоты колебаний, подлежащие определению.

Подставляя пробные функции x₁(t), x₂(t) в матричное уравнение, получаем характеристическое уравнение для определения собственных частот колебаний:

Решая это биквадратное уравнение, находим собственные частоты колебаний. Вычислим сначала дискриминант:

Тогда квадраты собственных частот будут описываться формулой

Далее, чтобы избежать громоздких формул, рассмотрим более простой случай, когда жесткость всех пружин одинакова: k₁ = k₂ = k₃ = k. Кроме того, введем отношение масс тел: μ = m₂/m₁. Тогда формула для квадрата частот колебаний принимает такой вид:

Полученное соотношение описывает 2 собственных частоты: ω₁ (со знаком "плюс") и ω₂ (со знаком "минус"). Зависимости частот ω₁, ω₂ от отношения масс μ показаны на рисунке 2. В случае равных масс (μ = 1) собственные частоты колебаний будут описываться следующими компактными формулами:

Заметим, что частоты ω₁, ω₂ всегда являются действительными числами. Это следует из общих физических соображений. Действительно, при мнимой частоте возникла бы утечка энергии, что противоречило бы условию сохранения энергии в системе. Данный факт, однако, можно установить и чисто математическим путем. В самом деле, вопрос возникает лишь для частоты ω₂. Условие неотрицательности для ω₂² выглядит так:

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

что всегда выполняется.

Найдем собственный вектор H₁ = (H₁₁, V₂₁)^T (верхний индекс ^T обозначает операцию транспонирования), соответствующий частоте ω₁. Он определяется из матрично-векторного уравнения

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

В последней системе оба уравнения являются линейно зависимыми (поскольку определитель матрицы K равен нулю при ω² = ω₁²). Поэтому, координаты собственного вектора H₁ можно выразить, например, из первого уравнения. Пусть H₁₁ = 1. Тогда

Таким образом, вектор H₁ имеет следующие координаты:

Аналогично можно определить собственный вектор H₂ = (H₁₂, V₂₂)^T, соответствующий частоте ω₂. В этом случае для H₂ имеем уравнение

В развернутом виде оно записывается как

Полагая H₁₂ = 1, из первого уравнения найдем координату H₂₂:

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

После того, как собственные частоты ω₁, ω₂ и собственные векторы H₁, H₂ найдены, можно записать общее решение системы. Учтем, что каждый из собственных векторов соответствует квадрату собственной частоты, т.е. двум значениям частот с противоположными знаками. Вектору H₁ соответствуют две частоты ± ω₁, а вектору H₂ − частоты ± ω₂. В результате общее комплексное решение представляется в виде суммы четырех слагаемых:

где C₁,...C₄ − постоянные (в данном случае комплексные) числа. Пусть данные числа записываются в виде

Чтобы величины x₁(t), x₂(t) оставались действительными при любом значении t, необходимо, чтобы выполнялись соотношения

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

Выражения в квадратных скобках можно упростить по формуле Эйлера:

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

Далее удобно ввести фазовые углы φ₁, φ₂ и воспользоваться тригонометрической формулой

В результате общее решение будет записываться в следующем виде:

где действительные постоянные A₁, A₂, φ₁, φ₂ зависят от начальных смещений и начальных скоростей тел, а собственные частоты ω₁, ω₂ и собственные векторы H₁, H₂ определяются соотношениями

Закон изменения скоростей тел находится путем дифференцирования общего решения:

Отсюда видно, что если в начальный момент t = 0 скорости тел равны нулю, то фазовые углы φ₁, φ₂ также равны нулю. Далее будем рассматривать именно такой случай. Общее решение представляет собой сумму двух гармоник с частотами ω₁, ω₂:

Вычислим постоянные A₁, A₂ в зависимости от начальных смещений. Пусть

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

Эту алгебраическую систему можно решить по формулам Крамера:

Итак, при начальных условиях

получаем следующую формулу общего решения:

где постоянные A₁, A₂ равны

а собственные векторы и собственные частоты выражаются через отношение масс μ, массу второго тела m₂ и коэффициент жесткости пружин k по приведенным выше формулам.

Полученные выражения значительно упрощаются, когда массы обоих тел одинаковы. Полагая μ = 1, получаем следующие формулы (при тех же самых начальных условиях):

Следовательно, в случае одинаковых масс и одинаковых коэффициентов жесткости законы движения тел определяются соотношениями:

закон движения тел, соединенных пружинами

В начале web-страницы представлена анимация, демонстрирующая характер колебаний тел, соединенных пружинами, при различных параметрах системы μ, k и начальных смещениях x₁₀, x₂₀. В модели принято значение m₂ = 2 (кг). Коэффициент жесткости k измеряется в Н/м. Смещения грузов показаны в сантиметрах в масштабе 1 см = 1 пиксель (график имеет масштаб 1 см = 3 пикселя).

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

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

Двойной маятник

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

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

Далее мы займемся построением математической модели двойного маятника в виде системы нелинейных дифференциальных уравнений. Начнем с вывода уравнений Лагранжа.

Уравнения Лагранжа

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

. Используя указанные переменные, построим лагранжиан двойного маятника и запишем дифференциальные уравнения Лагранжа.

Упрощенная модель двойного маятника показана на рисунке 1. Будем считать стержни невесомыми. Их длины равны l₁ и l₂. Массы точечных тел (они представлены шарами конечного радиуса) составляют m₁ и m₂. В точках подвеса трение отсутствует.


Рис.1		Рис.2

Введем систему координат Oxy, начало которой совпадает с точкой подвеса. Координаты маятников определяются следующими соотношениями:

Кинетическая и потенциальная энергия маятников (соответственно T и V) выражаются формулами

кинетическая и потенциальная энергии маятников

Тогда лагранжиан записывается в виде

Учтем, что

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

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

Теперь мы можем составить уравнения Лагранжа (иногда их называют уравнениями Эйлера-Лагранжа):

Входящие в уравнения частные производные выражаются следующими формулами:

Следовательно, первое уравнение Лагранжа записывается как

Сокращая на l₁ ≠ 0, получаем:

Аналогично выведем второе дифференциальное уравнение:

После сокращения на m₂l₁ ≠ 0 уравнение принимает такой вид:

Таким образом, нелинейная система двух дифференциальных уравнений Лагранжа записывается как

система нелинейных дифференциальных уравнений Лагранжа для двойного маятника

Малые колебания двойного маятника

Если считать углы α₁(t), α₂(t) малыми, то колебания маятников вблизи нулевого положения равновесия можно описать линейной системой уравнений. Чтобы получить такую систему, вернемся назад к исходному лагранжиану системы:

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

Здесь мы учли, что слагаемое с cos(α₁ − α₂) содержит произведение малых величин

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

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

квадратичный лагранжиан двойного маятника

Выведем дифференциальные уравнения Лагранжа для данного лагранжиана. Они записываются в таком виде:

Найдем частные производные:

Получаем систему двух дифференциальных уравнений Лагранжа:

или

Данную систему уравнений можно записать в компактной матричной форме. Введем матрицы

Тогда система дифференциальных уравнений представляется в виде

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

где H₁, H₂ − собственные векторы, ω − действительная частота. Значения нормальных частот ω_{1, 2} определяются из решения характеристического уравнения

Выведем общие формулы для циклических частот ω_{1, 2} в случае произвольных масс m₁, m₂ и длин l₁, l₂:

Мы получили биквадратное уравнение для частот ω. Вычислим дискриминант:

Таким образом, квадраты нормальных частот ω_{1, 2} равны

или

нормальные частоты колебаний двойного маятника

Данное выражение является несколько громоздким. Поэтому далее рассмотрим случай, когда длины стержней обоих маятников равны: l₁ = l₂ = l. Тогда нормальные частоты будут определяться более компактной формулой

Как видно, собственные частоты ω_{1, 2} зависят лишь от отношения масс μ = m₂/m₁. Зависимости частот ω₁, ω₂ от параметра μ (при условии g/l = 1) показаны выше на рисунке 2. В частности, при равных массах m₁ = m₂ = m, т.е. при μ = 1, собственные частоты равны

собственные частоты колебаний двойного маятника при равных масссах

Теперь, после того как собственные частоты ω_{1, 2} известны, для описания нормальных мод нужно определить еще собственные векторы H_{1, 2}. Они находятся из решения векторно-матричного уравнения

Пусть собственный вектор H₁ = (H₁₁, H₂₁)^T (верхний индекс ^T означает операцию транспонирования) соответствует нормальной частоте ω₁. Тогда получаем следующее уравнение для определения H₁:

Координаты собственного вектора H₁ удовлетворяют уравнению

Таким образом, собственный вектор H₁ равен

Аналогичным образом найдем координаты второго собственного вектора H₂ = (H₁₂, H₂₂)^T :

Следовательно, собственный вектор H₂ имеет такие координаты:

Общее решение матричного уравнения записывается в виде

где постоянные C₁, C₂, φ₁, φ₂ зависят от начальных положений и скоростей маятников.

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

В этом случае начальные фазы равны нулю: φ₁ = φ₂ = 0. Определим постоянные C₁ и C₂:

Тогда закон колебаний маятников выражается формулами

где циклические частоты ω_{1, 2} определяются соотношением

Здесь углы α₁(t), α₂(t) выражаются в радианах, а время t в секундах. На рисунках 3-5 приведены графики малых колебаний маятников для трех значений μ: μ₁ = 0.2, μ₂ = 1, μ₃ = 5, при условии l = l₁ = l₂ = 0.25 м, g = 9.8 м/c². Углы отклонения маятников для удобства приведены в градусах. Из графиков видно, что в системе происходят биения, при которых энергия циклически переходит от одного маятника к другому. Когда один маятник почти останавливается, другой раскачивается с максимальной амплитудой. Через некоторое время маятники "меняются ролями" и так далее. Колебания с большей частотой ω₁ модулируются более низкочастотными колебаниями с частотой ω₂. Это особенно хорошо заметно на рисунке 5 при большом значении μ (μ₃ = 5), когда разница между частотами ω₁ и ω₂ велика.


Рис.3		Рис.4


Рис.5		Рис.6

Итак, малые колебания двойного маятника имеют периодический характер и описываются суммой двух гармоник с частотами ω₁, ω₂, зависящими от параметров системы. Характерным свойством малых колебаний двойного маятника является эффект биений.

Преобразование Лежандра и уравнения Гамильтона

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

Приведенные выше уравнения Лагранжа являются дифференциальными уравнениями второго порядка. Их удобнее преобразовать в форму канонических уравнений Гамильтона. В результате вместо 2 уравнений второго порядка мы получим систему 4 дифференциальных уравнений первого порядка.

В гамильтоновой механике состояние системы определяется обобщенными координатами и обобщенными импульсами. В нашем случае в качестве обобщенных координат мы снова, как и в уравнениях Лагранжа, будем использовать углы α₁, α₂. Вместо обобщенных скоростей

(в лагранжиане) введем теперь обобщенные импульсы p₁, p₂, связанные со скоростями формулами

или в краткой записи:

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

Предположим, что f(x) − гладкая выпуклая вниз функция (рисунок 6). Рассмотрим прямую y = px, проходящую через начало координат. Расстояние между прямой y = px и функцией y = f(x) вдоль оси Oy зависит от координаты x. При определенном значении x это расстояние будет максимальным. Ясно, что оно зависит от наклона прямой, т.е. от параметра p. Таким образом мы вводим новую функцию g(p):

Такое преобразование функции f(x) в сопряженную функцию g(p) называется преобразованием Лежандра. Заметим, что функция g(p) достигает максимального значения по переменной x когда p = df/dx. Действительно,

Зная зависимость p(x), можно найти обратную функцию x(p). Тогда преобразование Лежандра будет выражаться соотношением

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

преобразование Лежандра для перехода от лагранжиана к гамильтониану

В этом выражении L является лагранжианом, а функция H представляет собой гамильтониан системы, который зависит от обобщенных координат α₁, α₂ и обобщенных импульсов p₁, p₂.

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

Определим теперь конкретный вид уравнений Гамильтона для двойного маятника. Обобщенные импульсы p₁, p₂ выражаются через частные производные лагранжиана в виде

Решим эту систему уравнений и выразим угловые скорости через обобщенные координаты и импульсы. Воспользуемся формулами Крамера и вычислим соответствующие определители:

Отсюда получаем следующие выражения для угловых скоростей:

Эти формулы представляют собой первые 2 (из 4) дифференциальных уравнений Гамильтона. С учетом данных выражений гамильтониан можно записать в следующем виде:

Последнюю формулу можно представить как

Числитель N в этом выражении является весьма громоздким. Упростим его:

Следовательно, функция Гамильтона принимает такой вид:

Здесь первое слагаемое представляет собой обобщенную кинетическую энергию T, а два других слагаемых − потенциальную энергию V, т.е. гамильтониан H определяется как

где

Теперь мы можем составить еще два дифференциальных уравнения Гамильтона для обобщенных импульсов:

Вычислим отдельно частные производные обобщенной кинетической энергии:

где символами A₁ и A₂ обозначены выражения

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

Отсюда получаем уравнения Гамильтона в виде:

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

система уравнений Гамильтона для двойного маятника

где

Теперь можно приступить к численному анализу уравнений.

Численное моделирование хаотических колебаний

Наиболее распространенным методом численного решения дифференциальных уравнений является метод Рунге-Кутты 4-го или 5-го порядка точности. Различные вариации этого метода используются в большинстве математических пакетов (MatLab, Maple, Mathematica, Mathcad), как правило, с автоматическим контролем точности и адаптивным временным шагом.

Для моделирования движения двойного маятника мы также воспользуемся классическим методом Рунге-Кутты 4-го порядка точности. Предварительно несколько упростим дифференциальные уравнения, полагая, что длины маятников одинаковы: l₁ = l₂ = l. Введем также параметр μ, равный отношению массы второго маятника к массе первого: μ = m₂/m₁. Тогда система уравнений принимает следующий вид:

где

Данную систему можно переписать в векторной форме:

Вектор Z составлен из 4-х канонических переменных данной системы, а компоненты вектора f соответствуют правым частям дифференциальных уравнений.

Метод Рунге-Кутты предполагает на каждом шаге последовательное вычисление четырех промежуточных векторов:

Значение вектора Z_n+1 в следующем временном узле вычисляется по формуле

Суммарная ошибка данного алгоритма на конечном интервале имеет порядок O(τ⁴), т.е. точность вычислений возрастает в 16 раз при уменьшении временного шага τ в два раза.

Описанная модель реализована в анимации, приведенной в начале web-страницы. Для упрощения мы положили начальные углы отклонения маятников равными: α₁ = α₂ = α. Данное приложение наглядно демонстрирует хаотическую динамику двойного маятника при различных значениях параметров μ и α. Интересно,