Название: Моделирование и анализ нелинейной динамики - Лабораторный практикум (Хиценко В.Е.)

Жанр: Технические

Просмотров: 1279


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

 

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

 

                                   Теоретическое описание

 

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

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

Состояние динамической системы в момент t характеризуется вектором  x(t) = [ x1(t), x2(t), … , xn(t)]. Элементами вектора могут быть, например, концентрации веществ в химической реакции, численности популяций в системе их взаимодействия, положение, скорость и высшие производные координат материального тела и т.п. Пространство векторов состояний M называется фазовым пространством. Рассмотрим в М автономную динамическую систему, описываемую дифференциальным уравнением

                                         .                                                                                   (3.1)     Решения   x(t) , полученные при всевозможных начальных условиях  x(0) ÎМ , образуют в М множество непересекающихся фазовых траекторий (фазовый портрет). Оператор эволюции ( здесь вектор-функция   f (x) = [f1(x), f2(x), … , fn(x)] ) формирует вектор скорости, с которым фазовая траектория выходит из состояния  x . Так что f (x) определяет поле скоростей конкретной системы в фазовом пространстве.

Механизм (3.1), однозначно определяющий динамику системы, доказывает существование отображения  Фt  фазового пространства М в себя

 

                              x(t) = Фtx(0) ÎМ ,                                                                               (3.2)

 

называемого потоком.  Этот поток переносит некоторую область фазового пространства W0 Ì  М в область  Wt = Фt W0  Ì  М  за время t .

Область W0  при этом в общем случае деформируется, но если  ее объем, равный  V0= , в среднем не меняется,  то система (3.1) называется консервативной или сохраняющей объем. Это эквивалентно отсутствию трения. Такие системы описывают движение планет, частиц в разреженном газе, волн в плазме.

 

Рис.3.1. Поток  Фt  в фазовом пространстве.

 

Наличие трения, рассеяния энергии свойствено реальным системам,  которые по этой причине принято называть диссипативными. Этот термин применяют и в немеханических системах. При этом происходит уменьшению фазового объема области Wt с течением времени, так что в пределе этот объем равен нулю:

                                  

                                    .

 

Относительное изменение фазового объема называется дивергенцией и равно

 

                              . 

 

Фазовый поток при отрицательной дивергенции называется сжимающим.  

Траектории втягиваются в некоторую зону В Ì М  нулевого объема, значит, размерность этой зоны, называемой аттрактором, меньше размерности пространства n . Именно наличие внутреннего трения (в самом широком смысле) влечет за собой существование аттрактора, т.е. асимптотического предела решений, предела на который не оказывает прямого влияния начальное условие – исходная точка.

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

Аттрактором диссипативной системы может быть точка в М. Это будет устойчивая неподвижная точка отображения Фt (устойчивый фокус, узел). Аттрактором может оказаться замкнутая кривая (предельный цикл), поверхность тора размерностью n-1. Это соответствует периодическим и квазипериодическим режимам. Эти аттракторы могут быть описаны алгебраическими уравнениями и являются дифференцируемыми многообразиями в М.

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

 Именно условие n>2 необходимо, чтобы исключить пересечение и самопересечение траекторий в этом аттракторе и тем обеспечить единственность бесконечной траектории, начинающейся в конкретной начальной точке x(0). 

 

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

 

Понятно, что на фазовой плоскости это невозможно.

Фазовые траектории диссипативной системы, начинающиеся  в некоторой области притяжения странного аттрактора (бассейн аттрактора), сближаясь между собой, вполне предсказуемо втягиваются в него и, оставаясь там навсегда, флуктуируют внутри аттрактора неожиданным образом. Две рядом идущие, близкие траектории вдруг экспоненциально быстро разбегаются к границам аттрактора, затем рано или поздно, хотя и неясно когда, вновь оказываются близкими и вновь расходятся. Проявляется уже отмеченная чувствительность к точности задания начальных условий, что принципиально ограничивает прогноз поведения системы. Даже имея точную модель реального явления, мы все равно не  сможем абсолютно безошибочно измерить и установить стартовые значения состояний для компьютерной имитации поведения с целью прогноза.  Таким образом, при вполне детерминированном механизме в виде (3.1) мы сталкиваемся с непредсказуемостью, свойственной хаосу и случайности. А.Пуанкаре в 1908 г., говоря об особенностях динамики неинтегрируемых систем, заметил, что “совершенно ничтожная причина вызывает значительное действие, которое невозможно было предусмотреть”.

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

На рис.3.3 показан аттрактор и решения системы дифференциальных уравнений

                                                                                                                (3. 3)

 

 

Рис.3.3. Странный аттрактор и соответствующие ему хаотические решения системы (3.3).

 

Некоторую закономерность, детерминированность в странном аттракторе можно усмотреть, тогда как решения системы (3.3) кажутся совершенно непредсказуемыми.

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

Все известные странные аттракторы трехмерного пространства имеют размерность 2 < df < 3.  Можно предположить, что размерность странного аттрактора в пространстве произвольно высокой размерности n также будет лежать в диапазоне n-1< df < n. То есть, поток Фt хаотической системы, сжимая фазовый объем по пути к аттрактору, все  же не превращает этот объем в поверхность. Странный аттрактор это меньше чем объем, но больше чем поверхность, именно это позволяет фазовым траекториям при n  ³ 3 бесконечно избегать пересечений.

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

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

Другим надежным критерием хаотичности поведения системы типа (3.1) является положительность одного из показателей Ляпунова. Каждая из составляющих вектора расхождения двух близких траекторий в аттракторе Dx(t) = çx1(t)-x2(t)ç аппроксимируется экспонентой Dxi(t)= Dxi(0)explit, i=1,…,n. Показатели этих экспонент дают спектр показателей Ляпунова системы : l1, l2, … , ln . В условиях диссипативного сжатия в спектре должны быть в основном отрицательные показатели Ляпунова, но для хаотичности поведения траекторий максимальный показатель должен быть положительным, поддерживающим экспоненциально быстрый разбег двух соседних траекторий вдоль одной из осей фазового пространства. Известны несколько методов вычисления показателей Ляпунова, часть из которых сводится к вычислению верхнего предела вида     

li = ½Dxi(t)½].

Дивергенция системы (3.3) равна

 

                              

то есть, она флуктуирует вместе с  x2 (t). Интегрируя систему (3.3) совместно с

 

                                        

мы можем увидеть как меняется во времени первоначальный фазовый объем V(0).

В этой лабораторной работе мы рассмотрим особенности моделирования и анализа динамических систем с помощью программы ODE. Описание этой программы приведено в приложении.

 

                   Порядок проведения лабораторной работы

 

Запустите программу ODE и откройте файл VANDPOL1.ode. В этом файле записано уравнение осциллятора Ван дер Поля

             в нормальной форме. С ростом α от 0.1 до 10 форма автоколебаний искажается. Посмотрите фазовый портрет и решения этой системы для различных α.

 

Откройте файл TOR3.ode и посмотрите в режиме имитации движения как наматывается фазовая траектория на тор в трехмерном фазовом пространстве. Посмотрите эту динамику во временной развертке. Тоже для четырехмерного тора (файл TOR4.ode).

Откройте файл LOTVOL2.ode, где имитируется динамика  взаимодействия хищников и жертв. Биологические модели Лоттки-Вольтерра имеют следующий вид

                         .

Здесь x это численность «зайцев», α – их рождаемость, β – коэффициент гибельности встреч с «волками», y – численность «волков», λ -  их смертность, μ – коэффициент прироста «волков» от хорошего питания.  Изменяя  эти параметры и начальные условия, увеличивайте минимальное число зайцев в колебаниях их численности, которые возникают в системе Лоттки-Вольтерра. Добейтесь наилучшего результата. Получите зависимости нижних пределов популяций жертв и хищников от параметров модели. Оцените влияние β на период колебаний.

Вычислите дивергенцию этой системы. Получите процесс изменения дивергенции и фазового объема. Сделайте вывод о характере системы Лоттки-Вольтерра (консервативна или диссипативна).

Посмотрите в режиме имитации движения ход фазовой траектории по аттрактору Ресслера (файл REUSSL.ode). Посмотрите решения системы Ресслера. Проверьте диссипативность этой системы.

Проделайте аналогичный анализ  для системы Лоренца (файл LORENZ.ode). Исследуйте влияние параметра r, первоначально равного 28 на динамику системы Лоренца.   ( 0,5 < r < 13,93; 24,06; 28,0). Убедитесь в гиперчувствительности хаотических решений к точности задания начальных условий. Для этого, используя возможности программы ODE, покажите на одном графике решение x(t) и то же решение, но при x(0) , отличающемся на 0,1\%.

Проведите аналогичный анализ  системы, выданной преподавателем.

 

                                         Содержание отчета

 

Уравнение Ван дер Поля в обычной и нормальной форме. Семейство предельных циклов и колебаний этого осциллятора для различных значений параметра α.

Уравнения Лоттки-Вольтерра. Колебания численности популяций, предельные циклы. Зависимости нижних пределов численности и периода от параметров модели взаимодействия.

Выражение дивергенции системы Лоттки-Вольтерра. Зависимость ее от времени. Вывод о диссипативности системы.

Уравнения Ресслера. Аттрактор этой системы. Процесс изменения дивергенции системы Ресслера. Вывод о характере системы.

Уравнения Лоренца. Решения уравнений. Аттрактор Лоренца. Анализ влияния параметра r. График, иллюстрирующий чувствительность к начальным условиям.

Результаты анализа индивидуальной системы.

 

                                 Контрольные вопросы

 

Определение динамической системы. Различные математические модели динамики.

Что такое теория динамических систем? Что такое поток непрерывной динамической системы?

О чем говорит дивергенция потока? Как ее вычислить?

Что называется аттрактором динамической системы?

В чем проявляется гиперчувствительность к начальным условиям?

Почему фазовые траектории автономных систем не могут пересекаться?

Как теоретически проверить устойчивость решений нелинейных систем дифференциальных уравнений?

Как перейти от дифференциальных уравнений к разностным при малом Δt? Как проверить диссипативность полученных разностных уравнений?

Как получить уравнение динамики изменения фазового объема?

Найдите стационарное решение уравнений Лотки-Вольтерра и проверьте его устойчивость.

То же для уравнений Лоренца.

Как вывести на экран изменение дивергенции и фазового объема?

Как найти стационарное решение системы нелинейных дифференциальных уравнений и проверить его устойчивость?

Как можно, используя программу ODE, получить все проекции трехмерного аттрактора на координатные плоскости?

 

                Приложение «Описание программы ODE»

 

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

Взаимодействие с программой:

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

- Интерактивная справка доступна, если нажать <F1>.

- При нажатии <F10> вызывается операционная система DOS, в то время как программа ODE остается в памяти и повторно загружается по команде "Выход".

ЛИНИЯ МЕНЮ в верхней части используется, чтобы выбрать различные средства:

- File name (имя файла): название загруженного файла; вызываются файлы с расширением ".ODE", файлы с другими расширениями игнорируются.

- Remark (примечание): Линия для информативного текста без влияния на программу.

- Load ODE (загрузка): Загружает файл " *. ODE " содержащий дифференциальные уравнения, параметры, и предыдущие решения.

- Save ODE (сохранить ODE): Создает файл " *. ODE " содержащий дифференциальные уравнения, параметры и предыдущие решения. Обновляет предыдущий файл или создает новый с использованием кнопки «TAB»

- Edit ODE (редактирование): Ввод дифференциальные уравнений, начальных состояний и параметров режима интегрирования (см. ниже).

Integrate (интегрировать): Начинает численное интегрирование (см. ниже).

Edit Representation (редактирование представления результатов): вид графиков и т.п. (см. ниже).

- Display (дисплей): Дает графическое представление решения или вывод

  данных в файле с расширением " *. TBL ", если активирована опция Tabl.

- Edit ODE (редактирование): активизирует центральную часть экрана. Нажимая <ESC>, можно выйти из режима редактирования.

Дифференциальные уравнения записаны как

                dxi/dt = fi (t, x1, .., xn), i = 1 .. n,

где xi – обозначает переменные.

- d (..) /dt: ввод названия переменной.

- f (....): ввод функции fi. (список разрешенных математических постоянных дан ниже).

Пример: дифференциальные уравнения

            dv/dt = - x - v

            dx/dt = v.

Требуют ввода

 1        v = -x-v

            2        x = v

Эти строки должны быть написаны в последовательном порядке без пустых строк. Строка 0 предопределена как dt/dt = 1 и не может быть изменена.

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

- Range of integration (характер интегрирования):

  - Start value (начальные состояния): ввод начальных значений переменных в том же самом порядке как в столбце d ( )/dt.

- Stop-difference (интервал времени интегрирования): если абсолютное значение разности между значением переменной t и начальной t0 превышает число, заданное здесь, интегрирование автоматически останавливается.

- Store at difference (сохранить в разности): Если абсолютное значение разности между значением переменной и последним сохраненным ее значением не превышает число, указанное здесь, то текущее значение не сохраняется. Если введен - "0", то сохраняются все значения в векторе решения. Т.е. это точность сохраненных результатов. Везде кроме режима «PoincarВ отображение» требуется, чтобы по крайней мере одна из величин, указанных здесь отличалась от ноля.

- PoincarВ map (Poincare отображение): из данных по непрерывной траектории решений сохраняются только точки пересечения с n-мерной (гипер)плоскостью в (n +1) – мерном пространстве решений. Эта плоскость дана в форме Гессе

                      "x" · "n" = "a",

где "x" – вектор решения, "n" - вектор, ортогональный к плоскости,

и "a" – постоянная. Первый номер в колонке «PoincarВ отображение» – это величина "a", затем идут составляющие вектора "n". Минимальное расстояние плоскости от начала координат  - "a"/ | "n" |.

  Если требуется использовать режим PoincarВ, то ВСЕ значения в столбце  «Сохранить в разности» должны быть нулевыми и "n" должен быть отличным от нулевого вектора.

- Accuracy/time unit (точность / единица времени): здесь, верхний предел допустимой ошибки дан относительно единицы времени. Переменные с назначенным допуском "0" не проверяются. Первая строка для переменной "t" не может быть введена и по крайней мере один из других входов должен быть отличен от ноля.

- Integrate (интегрируйте): Начинает интегрирование дифференциальных уравнений. Программа генерирует временный файл "ODE.TMP". Численное интегрирование может быть прервано при нажатии любой клавиши. Предыдущее интегрирование может быть продолжено. В таком случае необходимо выбрать один из следующих пунктов:

- Restart integration with start values (перезапустить интегрирование со старыми начальными условиями): интегрирование начинается снова с прежними начальными условиями.

- Restart integration with current values (перезапустить интегрирование с текущими величинами в качестве начальных).

- Continue previous integration (продолжить предыдущее интегрирование): вычисление решения продолжается и все вычисленные данные остаются сохраненными.

- Start new integration, hold old values (начать новое интегрирование с сохранением старых результатов (векторов решения)): интегрирование начинается снова, используя новые начальные значения, новые параметры или даже дифференциальные уравнения. При этом ранее вычисленные решения остаются сохраненными и могут быть сравнены с новыми на графике.

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

Edit Representation (редактирование представления результатов): активизирует нижнюю часть экрана. Нажимая <ESC>, возвращаемся к верхней строке экрана.

- <Representation mode> ( Вид демонстрации): Здесь мы выбираем содержание осей координат на экране. Вычисленные переменные решения, например,  x(t), y(t), ... могут быть показаны как функции времени, для этого оси Х присваиваем t, а оси ординат Y – значение соответствующей переменной (возможно выводить на ось координат функцию от нескольких переменных). Присваивая осям значения переменных x, y, z, получаем фазовые траектории системы.

  - "2" -мерный. " n / (n + k) ": одномерное отображение с запаздыванием на k итераций.

                         X -ось: " x (n) "; Y -ось: " x (n + k) ".

    Величина k может быть изменена (от 1 до 32); k = 1 предварительно установлено.

  -"2" -мерный: Плоское изображение.              

  -"3"-мерный параллельный: Трехмерное изображение, показанное в параллельном проектировании.

  -"3"-мерный центральный: Трехмерное изображение, показанное в центральном проектировании.

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

- Representation of solution (представление решения): формирование графика функций "X,Y,Z". Графическое представление или вывод данной последовательности векторов решения на принтер. Графически решение показано в двумерной или трехмерной декартовой системе координат (X, Y) или (X, Y, Z), которая спроектирована на экран. Координаты " X, Y, Z " могут быть определены как функции переменных.

- Range of representation (амплитуда представления): выбор интервалов переменных, показанных на экране. <Авто> опция – активизирует автоматический выбор интервалов, содержащих полное вычисленное решение.

- <Simulation of motion> (имитация движения): когда этот режим активен, можно наблюдать перемещение сферы по траектории. Скорость движения определяется заданной ниже записью «время = » как реальное мгновенное время. Возможно масштабирование времени, например, 0.1*t, ln (t) или любая функция. Это "время" никогда не выполняется в обратном направлении по кривой решения.

- <Trajectory> (траектория): опция для аппроксимации вычисленных точек решения. Используется сплайновый алгоритм.

- <Plot graphics> ( Чертят графику): опция для отображения решения с помощью графопостроителя. Это невозможно, когда режим «Моделирование движения» активен. Также возможен вывод на Epson HI-80 и HPGL-графопостроитель, а также на PostScript принтер.

- <Graphics Hardcopy> ( Графическая копия): опция для вывода решения на принтер. Это невозможно, когда режим «Моделирование движения» является активным.

- <Table file> ( Таблица): когда этот пункт активен, решение не выводится на экран. Вместо этого все векторы решения записываются в файл с расширением " *.TBL", что позволяет производить обработку результатов другой программой.

 

Список математических операций:

  ( )                      скобки

  + - * /              сложение, вычитание, умножение, деление

  ^                      возведение в степень

  sqr                   квадрат

  sqrt                  квадратный корень

  abs                  абсолютное значение

  sgn                 знак функции

  ln                   натуральный логарифм

  lg                   логарифм по основанию 10

  exp                экспоненциальная функция

  sin cos tan     тригонометрические функции

  asn acs atn      обратные тригонометрические функции

  sinh cosh tanh    гиперболические функции

  asinh acosh atanh   обратные гиперболические функции

  Hv             функция  Хевисайда

 round         округление к целому числу

 trunc           усечение к целому числу