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

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

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

Для эффективного использования моделей данного класса следует остановиться кратко на обзоре результатов исследований в этом направлении.

Мы остановимся на рассмотрении особенностей поведения нелинейной модели

, n=0,1,…. (7.1)

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

В приложениях к моделям биологических популяций переменная состояния x имеет тенденцию роста на последующем шаге для малых значений на предыдущем, и тенденцию уменьшения – для больших x. В начале координат f(0)=0; функция монотонно возрастает в диапазоне 0<x<A (f(x) достигает максимума в точке x=A); монотонное уменьшение происходит при значениях x>A. Функция может содержать один, либо несколько параметров, вариация которых оказывает существенное влияние на поведение (7.1). Иначе говоря, наблюдается высокая чувствительность модели к изменению параметров «настройки». Эти параметры, безусловно, имеют биологическую, экономическую, либо социологическую интерпретацию. Наиболее часто встречающиеся модели вида (7.1) содержат квадратичную нелинейность

, t=0,1,….. (7.2)

Это логистическое разностное уравнение. В пределе, когда b=0, с помощью (7.2) описывается чисто экспоненциальный рост популяции (для a>1). Для квадратичная нелинейность определяет функцию, обладающую свойством строгой выпуклости. Крутизна функции зависит от параметра . Вместо (7.2) часто используется ее каноническая форма, которая получается путем подстановки / a; t=n:

(7.3)

В дальнейшем мы будем также использовать запись (7.3), в которой a=4·r:

, (7.4)

поскольку в этом случае диапазон вариации r можно установить неравенство 0<r<1. Практические соображения, связанные с невозможностью получения отрицательных значений популяций на любом шаге, определяют ограничения диапазона изменения переменной состояния 0<x<1. В этой связи следует также отметить, что максимум f(x) соответствует x=0,5 и равен a/4. Уравнение (7.3) поэтому обладает нетривиальными свойствами, если а<4. Кроме того, все траектории сходятся в точке х=0, если а<1. Для уравнения (7.4) аналогичные свойства обеспечиваются при r<1 и r< .

В математической теории экологических процессов часто используется уравнение

. (7.5)

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

Конечно, моделями (7.2) и (7.5) список динамических систем, простых в аналитическом описании, но сложных в поведении, ввиду большого разнообразия динамических свойств при неизменных структурах, вовсе не ограничивается [35]. В этой же работе, а также в работах [46], [47] можно найти большое число различных приложений, свидетельствующих (по названию последней работы) о важном научном направлении, связанном с изучением хаотических процессов и фракталей.

Остановимся кратко на анализе динамических свойств процессов, описываемых уравнением (7.1). Прежде всего отметим, что возможно определить равновесные значения (или «фиксированные точки») для переменной состояния Х, если положить в установившемся режиме X_n+1=X_n=X*. В этом случае (7.1) становится алгебраическим нелинейным уравнением

, (7.6)

что соответствует нулевому значению роста популяции. Для решения (7.6) можно использовать графические методы, отображающие X_n в состояние X_n+1, которые представляются точками пересечения нелинейной функции с прямой, определяющей равенство X_n= X_n+1 (при одинаковых масштабах по осям абсцисс и ординат – это прямая в первом квадранте в декартовой системе координат под углом 45⁰ и проходящая через начало координат). Для простой строго выпуклой функции с одним максимумом, представляемой уравнениями (7.3) и (7.5), существуют две таких точки: тривиальное решение Х=0 и нетривиальное решение Х* (которое для модели (8.3) равно X*=1- ).

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

Наклон можно определить с помощью простого соотношения

(7.7)

Из (7.7) следует, что если наклон поддерживается в диапазоне от 45⁰ до -45⁰ (т.е. значения изменяются на интервале от 1 до -1), точка равновесия X* будет устойчивой. Она будет обладать свойством притяжения, т.е. являться аттрактором для любой из наблюдаемых траекторий. В частности, модель (7.3) с наклоном будет обладать свойством притяжения всех траекторий к точке в начале координат в интервале 0<x<1, если и только если 1<а<3. За пределами этого интервала параметр а определяет совершенно иные свойства модели.

Заметим, что координата состояния на (n+2)-ом интервале, согласно (7.1), равна

(7.8)

Используя вместо (7.8) выражение, соответствующее установившемуся режиму, мы получим алгебраическое уравнение

(7.9)

Пересечение кривой с прямой линией теперь будет происходить в двух точках, т.е. будет порождаться двухпериодический процесс. Если наклон характеристики (7.9), определяемый по уравнению (7.7),

(7.10)

находится в интервале 0<л⁽²⁾<1, что соответствует углу наклона от 0⁰ до 45⁰, будет наблюдаться устойчивый режим при нетривиальном решении (7.9). Если же углы наклона нелинейной характеристики в точке пересечения с прямой будут таковы, что по абсолютному значению |л⁽²⁾|>1, в системе (7.9) установятся двухпериодические колебания. Базовая точка при этом, соответствующая режиму л⁽¹⁾, является неустойчивой. В частности, двухпериодический цикл можно наблюдать в системе (7.3), если а=3.414. Тогда кривая пересекает прямую (т.е. получается решение (7.9)) в трех точках. Между точками, соответствующими устойчивому двухпериодическому процессу, лежит базисная фиксированная точка, в которой процесс неустойчив.

Дальнейшее увеличение параметра а порождает неустойчивость двухпериодического процесса и устойчивый четырехпереодический цикл, а затем порождает также устойчивые циклы с периодами 8,16,32,64,…,2ⁿ. Иначе говоря, каждый предшествующий процесс с периодом к с увеличением а становится неустойчивым, одновременно бифурцирует, в результате возникает устойчивый цикл периода 2к. Наконец, возникает режим детерминированного хаоса, который далее (при значениях а, близких к 4) будет иметь настолько «узкие» границы («окна») сохранения устойчивого t-периодического процесса, что они могут быть сравнимы с погрешностью «машинного нуля». В этом случае хаотический процесс практически становится стохастическим.

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

Классическим примером хаотических систем в гидромеханике являются турбулентные течения и влияние на них изменения параметра – числа Рейнольдса [4]. Система детерминированных уравнений Навье-Стокса, описывающих движения жидкости, при изменении параметра (числа Рейнольдса) позволяет получить резкие переходы из одного состояния в другое (например, переход от ламинарного течения к турбулентному и наоборот) [7]. Кроме традиционных уравнений Навье-Стокса, в работах Белоцерковского О.М. отмечается возможность использования моделей класса (7.1) для описания динамики турбулентности течений, что обеспечивается богатым спектром динамического поведения моделей и возможностью использования процедур пригонки их параметров к экспериментальным данным.

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

Представляет большой интерес логистическая модель следующего вида:

N_t+1=N_t(1 + r(1-N_t /K)), t=0, 1, 2, 3, …. (7.11)

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

В рассматриваемой же ситуации дискретное отображение в форме (7.11) непосредственно определяется природными свойствами экологических процессов динамики популяций.

Динамические свойства дискретной системы (7.11) необходимо исследовать при известных начальных условиях N₀, что не должно вызывать в общем серьезных вычислительных трудностей. Однако, несмотря на кажущуюся простоту, итерационные отображения обладают внутренней (структурной) сложностью, которая определяется параметрами (7.11), в частности, коэффициентами r и К. Этим нелинейным отображениям было посвящено большое число публикаций, в том числе-с совершенно неожиданными результатами. Один из наиболее интересных и показательных результатов состоит в обнаружении свойств универсальности отображений. Универсальность проявляется в переходе от регулярных к хаотическим процессам, осуществляемом по так называемому “сценарию Фейгенбаума”. Предположим, что К=1. Тогда мы можем записать два дискретных квадратных уравнения:

N_t+1 = N_t+ rN_t(1-N_t), t = 0, 1, 2, 3, …, (7.12)

y_t+1=y_t²+c, t=0, 1, 2, 3, …, (7.13)

где с - постоянный коэффициент.

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

с = , . (7.14)

Если начальное условие

- r · N₀,

то

- r · N_t, t=1,2,3, …. (7.15)

Произведем подстановку (7.15) в уравнение (7.13). В результате получим

N_t+1 = N_t + r N_t (1-N_t) + K,

где

К = - - .

Приравнивая коэффициент К к нулю, будем иметь

с = ,

при котором

N_t+1 = N_t + r N_t (1-N_t), t = 0, 1, 2, 3, …

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

, a = 1+r.

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

x_t+1 = ax_t(1-x_t), t = 0, 1, 2, 3, … (7.16)

Рассмотрим график, соответствующий уравнению (7.16) и представленный непрерывной функцией

y=ax(1-x)

Это парабола, проходящая через точки (0.0) и (1.0), независимо от заданного параметра а. Максимальное значение y всегда расположено в точке с координатой по оси х=0.5. Ордината точки y=а/4. Парабола y=f(x) является верхней границей для любых процессов, описываемых логистическим уравнением (7.16), где х_t принимает значения на рабочем интеграле 0<x_t<1, при условии, что х₀ также принадлежит этому интервалу. Теперь для построения итеративного процесса в виде последовательности x₀, x₁, x₂, …на плоскости можно предложить простой графический метод, состоящий в следующем. Построим график параболы y=ax(1-x) и изобразим на плоскости y=f(x) прямую y=x, разделяющую плоскость пополам и проходящую через точки (0,0), (1,1). Предположим, что начальное условие x₀ 0. Проведем прямую, параллельную оси y и проходящую через x₀, до пересечения с параболой. Через полученную точку проведем прямую, параллельную оси x, до пересечения с прямой y=x. В результате мы получим абсциссу точки процесса на первой итерации x_t=x₁. Вновь проведем через эту точку прямую, параллельную y, до пересечения с параболой и т.д. В процессе построения траектории все точки на прямой должны располагаться на одинаковом расстоянии от осей y и х.

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

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

%File "SAH46.M".

%Period-Doubling Bifurcations.

%Iteration of logistic map.

r=0.95;

%---------------------------

x=0:0.01:1;

y1=4*r.*x.*(1-x);

y2=x;

plot(x,y1,x,y2)

hold on,

%---------------------------

x=0.01;

y=0.0;

s=[];

s1=[];

for n=1:500;

s=[s;x;x];

y=4*r.*x.*(1-x);

s1=[s1;x;y];

x=y;

end

plot(s,s1),grid

pause,

hold off

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

Вторая часть начинается с выбора начального условия х=0.01 (на нулевом шаге). Здесь вводятся «пустые» вектора S и S₁, которые заполняются на каждом шаге цикла, организованного с помощью операторов for, end. По завершении расчетов представляются графики динамических процессов в координатах x_t+1=f(x_t) для различных значений параметра а, который в программе записан в виде равенства а=4r.

Для моделирования процесса с помощью приведенного выше файла требуется задать значения r, х во второй части программы, а затем назначить число шагов t. В представленном варианте файла r=0.95, x=0.01 и число шагов t_max=500.

Рис. 7.1. Устойчивый процесс динамики логистической системы

r=0.7, x₀=0.01, f_max=500

Программа завершается оператором hold off, который «отменяет» условие нанесения различных графиков на один рисунок.

Приведем несколько характерных решений, полученных в процессе моделирования с помощью файла sah46.m.

На рис. 7.1 изображен переходный процесс в динамической системе (7.16) при следующих значениях параметров: r=0.7, x₀=0.01. t_max=500.

Стрелками показано направление движения изображающей точки по траектории. Процесс начинается с пошагового движения в направлении точки устойчивого равновесия. На 14-ой итерации переменная состояния достигает значения x₁₄=0.6364, а на 450-ой итерации с точностью до четвертого знака после запятой

x_t=x_t+1=0.6429.

Для представленного на графике процесса характерно

Рис. 7.2. Процесс с устойчивым предельным циклом

r=0.8, x₀=0.01, f_max=500

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

Увеличение коэффициента до значения r=0.8 приводит к заметным изменениям динамических свойств системы. Представленный на рис.7.2 переходный процесс (при неизменных начальных условиях) содержит устойчивый предельный цикл, на который изображающая точка переходит “изнутри”, т.е. из области пересечения граничных характеристик.

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

Рис. 7.3. Увеличение амплитуды предельного цикла

r=0.85, x₀=0.01, f_max=500

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

Дальнейшее увеличение параметра системы до величины r=0.90 приводит к возникновению хаоса (рис. 7.4). Если вывести на дисплей вектор S1 (см. файл sah46.m), то можно убедиться в сложности определения периодичности решений. Это свидетельствует о неравномерных интервалах времени достижения граничных значений переменными состояния. Хаотическое поведение характеризуется крайней степенью неустойчивости движения.

Рис. 7.4. Хаотический процесс

r=0.90, x₀=0.01, f_max=500

Дискретные модели логистических систем не являются тривиальным разностным аналогом соответствующих дифференциальных уравнений, а возникают при изучении заметно синхронизированных процессов [35]. На рис. 7.5 представлен хаотический процесс, полученный в системе, когда r=95. Как и на предшествующем рисунке, мы наблюдаем процесс с возрастающим числом решений (t-периодический процесс с числом периодов, составляющим несколько тысяч), которые чрезвычайно чувстви-тельны к изменению начальных условий и флуктуации параметров модели. Хаотический режим напоминает шумоподобный процесс в системе.

Рис. 7.5. Расширение области аттрактора: r=0.95, x₀=0.01, f_max=500

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

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