Определение показателей Ляпунова на примере модели Селькова в присутствии внешней периодической силы Текст научной статьи по специальности «Математика»

УДК 519.8

ОПРЕДЕЛЕНИЕ ПОКАЗАТЕЛЕЙ ЛЯПУНОВА НА ПРИМЕРЕ МОДЕЛИ СЕЛЬКОВА В ПРИСУТСТВИИ ВНЕШНЕЙ ПЕРИОДИЧЕСКОЙ СИЛЫ

В работе обсуждаются вычислительные особенности расчёта показателей Ляпунова на примере системы Селькова с периодическим втоком, которая описывает гликолитические колебания. Детально описана программная реализация алгоритма Вольфа в среде MATLAB и выбор его параметров. Численное решение рассматриваемой системы иллюстрирует существование трёх возможных режимов: а) гармонические колебания, представляющие собой предельные циклы на фазовой плоскости, внутри областей захвата; б) хаотический режим со странными аттракторами на фазовой плоскости между языками Арнольда; в) квазипериодические колебания с двумерными торами в качестве фазовых портретов на границах областей захвата.

Ключевые слова: хаотические колебания, показатели Ляпунова, модель Селькова, синхронизация, языки Арнольда, гликолитическая реакция, MATLAB.

При исследовании нелинейных систем одной из важных задач является определение типа колебаний - периодического, квазипериодического, случайного, хаотического. Особенно сложно отличить квазипериодические колебания от хаотических и случайных, так как квазипериодические колебания часто имеют очень сложную форму, визуально слабо отличимую от «случайной». В настоящее время существуют различные критерии определения хаоса [Лоскутов, Михайлов 2007]. Простейшим методом является исследование спектра колебаний на основе анализа Фурье. Дискретность спектра идентифицирует периодические или квазипериодические колебания, в случае непрерывности спектра колебания являются либо хаотическими, либо случайными. В качестве альтернативы анализу Фурье также успешно применяется вейвлет-анализ динамических систем. Другой метод основан на применении отображений Пуанкаре, то есть сечений фазовой траектории при помощи секущей поверхности. Отображение Пуанкаре случайного процесса будет иметь вид облака, а для квазипериодических и хаотических решений - форму некоторой линии.

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

© 2013 А. Ю. Верисокин

аспирант каф. общей физики e-mail: ffalconn@mail. ru

Курский государственный университет

ВВЕДЕНИЕ

= Fi (^..^Xn), i = U.П

(l)

существуют n показателей Ляпунова, определяемых формулой

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

показателей Ляпунова состоит в том, что два решения, начальные значения которых расположены в некоторой окрестности радиуса £, за время T разойдутся в n -мерный эллипсоид по n главным полуосям и в момент времени t радиусы будут определяться значениями ee^, i = 1,...,n [Кузнецов 2006]. Знаки показателей Ляпунова полностью характеризуют тип колебаний решения динамической системы. Наличие положительного показателя является критерием хаотичности решения динамической системы.

Для большинства динамических систем расчёт показателей Ляпунова возможен только численно. В настоящее время существует несколько алгоритмов. Наиболее важно определение старшего (наибольшего) показателя Ляпунова, так как именно он описывает тип колебаний. Для его вычисления обычно используется алгоритм Бенеттина [Benettin et. al. 1980]. Так как компонент решения, отвечающий старшему показателю Ляпунова, является доминирующим по величине, то для вычисления младших показателей приходится использовать специальные методы, один из наиболее точных и надёжных из которых основан на ортогонализации Грама - Шмидта [Wolf et. al. 1985].

Предложенная в цитируемой работе программная реализация приведена на языке FORTRAN, который в настоящее время используется для математического моделирования достаточно редко, в отличие от программной среды MATLAB. Один из вариантов трансляции алгоритма Вольфа для MATLAB был предложен в рамках базирующегося на нём модуля для исследования динамических систем MATDS [MATDS 2004]. Однако этот модуль использует специальный графический интерфейс, не дающий полного доступа к параметрам алгоритма при задании системы исследуемых дифференциальных уравнений, а также включает в подпрограмму обращение к ряду внутренних переменных модуля.

Поэтому практическую пользу представляет реализация в MATLAB алгоритма Вольфа «в чистом виде», в форме, явным образом идентифицирующей контролирующие параметры. Для её тестирования проводится поиск показателей Ляпунова на примере математической модели Селькова фосфофруктокиназной фазы гликолитической реакции [Selkov 1968] с добавлением внешнего периодического возмущения и исследование типа решений при различной его интенсивности.

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

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

уравнений:

dx 2 dy 2 «ч

— = v - xy2, -f = xy2 - wy. (3)

dt dt

Здесь x - концентрация субстрата реакции (АТФ), у - концентрация продукта (АДФ), v - скорость втока субстрата, w - скорость оттока продукта, необратимо удаляемого из автокаталитической реакции.

Проведём исследование системы (3). Стационарные значения переменных: xs = w2 /v, ys = v/w. Пусть x = xs +g, y = ys +n, где n - отклонения от

стационарного состояния. Обозначим f (x, y ) = v - xy2, g(x, y ) = xy2 - wy. Тогда,

используя разложение в ряд Тейлора, получим

^ = /1 (X, у, )£+/ (X, у, )п = - ^2 £-2 щ,

Ш М

ШП = 8X (X,у,)£ + 8У (X,У,)п = £ + МП-

аг м

\у

Рис. 1. Бифуркационная диаграмма системы Селькова, у1 = ^wJW, v2 = w/2

Характеристическое уравнение имеет вид Л2 - (w - V2 / м>2 )Л + 2v3 / м>2 - V2 / w = 0,

дискриминант В = м*2 + IV2 / м - 8у3 / м*2 + V4 / м4, то есть точки бифуркации системы (3): у1 = , у2 = w /2. В зависимости от знака В бифуркационная диаграмма системы

Селькова примет вид, представленный на рис. 1. Точка бифуркации Хопфа системы (1): V = w 4^. В зависимости от параметров V и н в системе возможно возникновение двух периодических типов колебаний: гармонических и релаксационных.

Рассмотрим поведение системы при фиксированном значении параметра оттока = 2 . Для w = 2 точкой бифуркации Хопфа будет V « 2.83 , то есть при V < 2.83 стационарная точка является неустойчивым фокусом и решение системы (3) имеет вид предельного цикла. С уменьшением V стационарная точка меняет свой тип сначала на неустойчивый узел, а затем - на седло, в результате чего предельный цикл исчезает.

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

V = у0 + А 8т(2П / Т'), (4)

где V - значение втока вблизи точки бифуркации Хопфа (V = 2.78), А - амплитуда колебаний втока, а Г' - период. При значении А = 0 система имеет форму (3) и её решение представляет собой квазигармонические колебания с периодом Т0 = 3.3. Если А ^ 0, то система принимает вид

Жх . . _ ,ч 2 Су

— = у0 + А зт(2П / Т') - ху , — = ху - wy.

(5)

С сИ

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

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

гармонический осциллятор и осциллятор Селькова, также можно ожидать большого разнообразия динамических режимов в зависимости от параметров A и T’.

МЕТОДЫ

Аналитическое определение показателей Ляпунова для системы (5) представляется невозможным из-за отсутствия аналитического решения указанной системы нелинейных дифференциальных уравнений. Поэтому была использована MATLAB-реализация алгоритма Вольфа [Wolf et. al. 1985] (код программы приведён в приложении). В программном коде были частично использованы функции из модуля для исследования динамических систем MATDS [MATDS 2004], однако в данной реализации используются только переменные, непосредственно требуемые алгоритмом.

Метод основан на численном решении системы совместно с уравнениями в вариациях, которые описывают эволюцию бесконечно малого возмущения траектории. Другими словами, рассматриваются две траектории, удалённые друг от друга на малое расстояние R, которое через промежуток времени 6t достигает значения R. Максимальный показатель Ляпунова на каждом шаге находится по формуле Лпах = log21 R1/ R, | / dt. С увеличением времени вычисления расчётное значение

показателя Ляпунова приближается к реальному. Для реализации описанного алгоритма используется программная среда математического пакета MATLAB.

В программе система Селькова задаётся функцией Selkov, в которую в качестве входных параметров передаются параметры системы, начальные значения переменных и временная переменная. В теле функции заданы правые части системы (5) и якобиан системы, на основе которых строятся уравнения в вариациях. Поиск показателей Ляпунова производится функцией lyapunov, в которую передаются параметры исследуемой системы, начальные значения, время расчёта и параметры решателя ОДУ. В блоке 1 функции lyapunov происходит инициализация необходимых при вычислении переменных. Непосредственное определение показателей Ляпунова происходит в основном цикле. В блоке 2 производится численное решение исследуемой системы с уравнениями в вариациях, которые позволяют определить отклонения траекторий. В блоке 3 происходит процесс ортогонализации Грама - Шмидта решения системы, позволяющий выделить все показатели Ляпунова и исключить влияние доминирующего старшего показателя на младшие. После этого подсчитываются суммы ортогонализованных отклонений. Наконец, в блоке 4 находятся показатели Ляпунова перенормировкой полученных сумм по времени вычисления. Цикл продолжается большое число раз, что позволяет точно приблизить расчётными величинами реальные значения показателей Ляпунова.

РЕЗУЛЬТАТЫ ЧИСЛЕННОГО АНАЛИЗА

Систему (5), несмотря на то, что она состоит из двух ОДУ, необходимо считать трёхмерной, так как параметр втока v, в отличие от системы (3), не является независимым, а колеблется периодически во времени по закону (4). Это означает, что система (5) имеет три показателя Ляпунова, один из которых всегда равен нулю в связи с периодичностью v . В зависимости от значений оставшихся показателей фазовые траектории системы могут относиться к одному из трёх типов: показатели

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

х ю"4

а) 1 хЮ5

Рис. 2. Вычисление старшего показателя Ляпунова системы Селькова (5) при значениях параметра периодического втока А = 1, Т' = 8 и начальных условиях х0 = 1.5, у0 = 1.4: а) линейные координаты, б) логарифмические координаты, в) вычисления показателей Ляпунова при тех же параметрах и

начальных условиях х0 = 2, у0 = 2

На рисунках 2 а), б) показан пример вычисления старшего показателя Ляпунова для системы (5) при значениях параметров периодического втока А = 0.1, Т' = 8 и начальных условиях х0 = 1.5, у0 = 1.4. Шаг итераций брался равным значению

собственного периода системы Т', что соответствует нахождению значения показателя Ляпунова, определяющего расхождение траекторий за один период колебаний. Абсолютная и относительная погрешности вычислений составляют порядка 10-10. Из рисунка 2 видно, что показатель Ляпунова при указанном наборе параметров положителен (кривая с течением времени сходится к малому положительному значению).

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

а)

б) * хЮ4

в) 1

Рис. 3. Поиск показателей Ляпунова при шаге итераций: а) 0.01, б) 2; в) сравнение результатов вычисления показателей Ляпунова при Ыз1ер=0.01 и Ыз1ер=2

На участках квазипериодических колебаний величина отклонения траекторий остаётся практически постоянной, а так как Атах обратно пропорционально дХ, то

мгновенное значение показателя Ляпунова уменьшается по степенному закону. Такие интервалы степенной зависимости могут быть наглядно представлены в логарифмических координатах - график в этом случае принимает вид наклонной прямой (рис. 2 б).

Тем не менее дальнейшие вычисления показателей Ляпунова для времён, больших /»105, демонстрируют преобладание хаоса в системе: вычисляемый показатель Атах быстро стремится к постоянному положительному значению. Графики в логарифмических координатах, более наглядно показывающие малые отклонения, подтверждают это.

На нижнем графике на рисунке 2 в) показан пример вычисления старшего показателя Ляпунова Атах для тех же значений параметров втока, что и на рисунках

2 а), б), но при других начальных условиях. Показатель Атах сначала имеет

отрицательное значение, но затем пересекает ноль и стабилизируется в области положительных значений при том же значении, что и на рисунке 2 а), что согласуется с теорией: значения показателей Ляпунова в общем случае не зависят от начальных условий.

При нахождении старшего показателя Ляпунова системы (3), результаты которого представлены на рисунке 2, шаг итераций ЫБ1ер брался равным 8 (период колебаний втока). Данный параметр определяет, на каких промежутках времени будут вычисляться отклонения траекторий. Выбор Ь1Б1ер равным периоду втока обоснован тем, что при анализе типа решения системы с внешней силой важным представляется определение влияния силы на протяжении полного периода колебаний, то есть величины отклонения траекторий за целый период. Кроме того, такой выбор шага аналогичен выбору плоскости сечения Пуанкаре при одном и том же значении величины втока. На вставке рисунка 1 а) показана динамика изменения вычисляемых показателей Ляпунова во времени с шагом итерации, равным 8. При меньших значениях шага итераций численное значение сходится к тому же значению, что и на рисунке 2. На рисунке 3 представлены примеры поиска максимального показателя Ляпунова при ЫБ1ер=0.01 и ЫБ1ер=2. На рисунке 3 в) представлена сравнительная динамика изменения расчётных значений показателя Ляпунова в окрестности ^ « 500. Как видно из графиков, несмотря на более высокую мгновенную точность вычисления показателей Ляпунова при ЫБ1ер=0.01 по сравнению с ^1ер=2, выбор большего значения параметра Ь1Б1ер не влияет на точность расчётов, так как вычисление показателей Ляпунова проводится в среднем по фиксированному времени вычислений (в одинаковые значения времени значения показателя при разных ЫБ1ер равны - рис. 3 в)). Выбор маленького шага итераций существенно увеличивает время расчёта показателей Ляпунова. В связи с этим при вычислении показателей Ляпунова параметр ^1ер целесообразно брать достаточно большим (как на рис. 2), при этом необходимо, чтобы амплитуда возмущения траекторий была меньше линейных размеров неоднородностей фазового пространства и, очевидно, размеров самого аттрактора.

На рисунке 4 приведена карта распределения старшего показателя Ляпунова системы (5). При численном анализе амплитуда втока изменялась в пределах от 0 до 1.5 с шагом 0.01, период колебаний величины втока брался в пределах от 0 до 10 с шагом

0.1. Участки, на которых старший показатель Ляпунова отрицательный (холодные тона), - области захвата собственных колебаний системы периодическим втоком (языки Арнольда).

Вблизи значений периода Т'/ т = Т0/ п, где п, т Е • в системе происходит захват колебаний втоком и наблюдаются синхронные колебания порядка п: т (п колебаний втока соответствуют т колебаниям решения системы). Максимальная ширина области синхронизации наблюдается при Т' » Т0.

0123456789 10

Т

Рис. 4. Карта распределения старшего показателя Ляпунова системы Селькова (5)

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

Отдельного внимания заслуживают области, возникающие внутри языков Арнольда при больших значениях втока, в которых показатели Ляпунова являются положительными (для языка Арнольда порядка 1:1 такая область возникает в районе А = 0.4, для 1:2 - А = 0.6 , для 1:3 - А = 0.7 ). На этих участках большое значение втока меняет тип стационарной точки системы на устойчивый фокус (см. рис. 1) на продолжительные промежутки времени. Предельный цикл уже не успевает притянуть к себе фазовые траектории за промежутки времени, когда параметр втока меньше значения в точке бифуркации Хопфа. Решения выходят из области притяжения предельного цикла, то есть расходятся, поэтому на таких участках не наблюдается колебательной динамики и определение показателей Ляпунова в таких областях смысла не имеет.

ВЫВОДЫ

Расчет показателей Ляпунова, вычисляемых с помощью алгоритма, предложенного в [Wolf et. al. 1985], допускает эффективную реализацию в среде

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

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

1. Периодические колебания внутри языков Арнольда (Amax < 0). Имеет место как захват частот, так и захват фаз колебаний. Фазовые портреты колебаний в этих областях имеют форму устойчивых предельных циклов.

2. Хаотический режим между областями захвата колебаний (Amax > 0), соответствующий странным аттракторам на фазовой плоскости.

3. Квазипериодические колебания на границах доменов (Amax = 0). Фазовые портреты системы будут представлять собой устойчивые двумерные торы.

Работа выполнена при поддержке фонда некоммерческих программ Дмитрия Зимина «Династия».

Библиографический список

Benettin, G., Galgani L., Giorgilli A., Strelcyn J.-M. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems: Acomputing all of them. P. I: Theory. P. II: Numerical application // Meccanica, 1980. V. 15. P. 9-30.

MATDS, 2004. URL: http://www.mathworks.com/MATLABcentral/linkexchange/links/1253-matds-toolbox (дата обращения: 12.09.2012)

Selkov E.E. Self-oscillations in glycolysis. A simple kinetic model // Eur J. Biochem, 1968. V. 4. P. 79-86.

Verisokin A.Yu., Verveyko D.V. Non-Turing Mechanism of Self-Sustained Structures Formation // International Journal of Bifurcation and Chaos, 2013. V. 23. No. 2. 1350037.

Wolf A., Swift J., Swinney H., Vastano J. Determining Lyapunov exponents from a time series // Physica D., 1985. V. 16. P. 285-317.

Кузнецов С.П. Динамический хаос. М.: Физматлит, 2006. 355 с.

Лоскутов А.Ю., Михайлов А.С. Основы теории сложных систем. М. - Ижевск: Институт компьютерных исследований, 2007. 620 с.

ПРИЛОЖЕНИЕ: ПРОГРАММНЫЙ КОД

x0=2; y0=2; % начальные значения tend=5e5; % время расчёта

% параметры системы Селькова v=2.78; % среднее значения втока субстрата w=2; % скорость оттока продукта

A=0.1; % амплитуда колебаний втока T=8; % период колебаний втока

Ltstep=8; % временной шаг итераций вычисления показателей Ляпунова

[Texp,Lexp]=SelkLyapExp([x0 y0], [v w A T], tend,Ltstep); plot(Texp,Lexp(:,1),Texp,Lexp(:,2));

title('Dynamics of Lyapunov exponents'); xlabel('t'); ylabel('\lambda_1, \lambda_2');

FUNCTION [Texp,Lexp]=SelkLyapExp(InVal, Param, TimeInt, Ltstep);

% Вычисление показателей Ляпунова системы Селькова % dx/dt=v(t)-xyA2; dy/dt=-wy+xyA2

% Периодический вток v(t)=v+A*(sin(2*pi*t/T))

% Входные параметры:

% InVal - начальные значения: массив [x0 y0]

% Param - параметры системы Селькова: массив [v w A T]

% TimeInt - время расчёта

% Ltstep - временной шаг итераций вычисления показателей Ляпунова % Выходные параметры:

% Texp - переменная времени, на котором происходил расчёт показателей Ляпунова % Lexp - показатели Ляпунова для каждого значения

odeParam=[1e-10 1e-10 1e-5 1e-6]; % параметры солвера

[Texp,Lexp]=lyapunov(2,@Selkov,@ode45,TimeInt,Ltstep,InVal,odeParam,Param);

END

FUNCTION f=Selkov(t,X,Param);

v=Param(1); w=Param(2); A=Param(3); T1=Param(4); % параметры системы x=X(1); y=X(2); % начальная точка траектории решения системы

Y= [X(3), X(5); X(4), X(6);]; % квадратная матрица размерности Якобиана f=zeros(6,1);

% правая часть системы

f(1)=v+A*(sin(2*pi*t/T1))-x*yЛ2; f(2)=x*yЛ2-w*y;

Jac=[^2,-2*x*y^2, 2*x*y-w]; % Якобиан системы f(3:6)=Jac*Y;

END

FUNCTION [Texp,Lexp] = lyapunov(n, func, ode, TimeInt, stept, ystart, odeParam,

SParam);

% n - число ОДУ в системе

% func - название функции, описывающей расширенную систему % ode - название солвера ОДУ

% ystart - начальная точка траектории решения системы

% БЛОК 1

Lexp=[];

Texp=[]; tstart=0; tend=TimeInt; options =

odeset('RelTol',odeParam(1),'AbsTol',odeParam(2),'MaxStep',odeParam(3),'OutputFcn',@ode

outp,'Refine',0,'InitialStep',odeParam(4));

n2=n*(n+1) ; % общее количество ОДУ, n - число нелинейных ОДУ nit = round((tend-tstart)/stept)+1; % число шагов

% инициализация промежуточных переменных y=zeros(n2,1); cum=zeros(n2,1); y0=y; gsc=cum; znorm=cum;

% начальные значения y(1:n)=ystart(:);

for i=1:n

y((n+1)*i)=1.0;

end;

t=tstart;

% основной цикл

while t<tend

% БЛОК 2

tt = t + stept; if tt>tend, tt = tend; end;

% решение расширенной системы [T,Y] = ode(@(t,X) func (t,X,SParam),[t t+stept],y); t=tt;

y=Y(size(Y,1),:);

% БЛОК 3

% создание ортогонального базиса с помощью преобразования Г рама — Шмидта % нормализация первого вектора znorm(1)=0.0; for j=1:n

znorm(1)=znorm(1)+y(n+j )A2; end;

znorm(1)=sqrt(znorm(1)); for j=1:n

y(n+j )=y(n+j )/znorm(1); end;

for j=2:n

% создание нормировочных коэффициентов преобразования for k=1:(j-1) gsc(k)=0.0;

for l=1:n gsc(k)=gsc(k)+y(n*j+l)*y(n*k+l); end; end;

% построение нового вектора for k=1:n

for l=1:(j-1) y(n*j+k)=y(n*j+k)-gsc(l)*y(n*l+k); end; end;

% вычисление нормы векторов znorm(j)=0.0;

for k=1:n znorm(j)=znorm(j)+y(n*j+k)A2; end; znorm(j )=sqrt(znorm(j));

for k=1:n y(n*j+k)=y(n*j+k)/znorm(j); end; % нормализация вектора end;

for k=1:n cum(k)=cum(k)+log(znorm(k)); end; % обновление векторов % БЛОК 4

% нормализация экспонент for k=1:n

lp(k)=cum(k)/(t-tstart);

end;

% Выходные данные - показатели Ляпунова Lexp=[Lexp; lp];

Texp=[Texp; t]; end;

END