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

С.А. Катрич

О КОМПЬЮТЕРНОМ АНАЛИЗЕ ПОВЕДЕНИЯ СТАЦИОНАРНЫХ СОСТОЯНИЙ

СИСТЕМЫ ЛОРЕНЦА

В работе ставится следующая цель. С помощью разработанных в [1] условий устойчивости по Ляпунову решений систем обыкновенных дифференциальных уравнений (ОДУ), основанных на разностных схемах, исследовать устойчивость стационарных (иными словами, тривиальных) решений системы Лоренца, которая исторически известна как первый объект с крайне нетривиальным поведением своих решений.

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

Описание условий устойчивости решений систем ОДУ. Пусть требуется исследовать устойчивость в смысле Ляпунова решения задачи Коши для нелинейной системы (ОДУ) в нормальной форме

с1У

— =Р«,Г),¥«0)=¥0, (1)

а ^

где У=У(0, г о (МО, У2(0, - искомое решение, 70 = (у^-у2ф„Ф ^ -

вектор начальных данных, /■ ( /. 7) = (/] (/. ¥ ). / 2 ((■ ¥ )....(/. ¥) - заданная вектор-функция от п + 1 переменных: независимой переменной I и п зависимых переменных 1, у2>, . . ., у,

n •

Предполагается, что для системы (1) выполнены все условия существования и единственности решения 7 =7(7) на всей полупрямой [0. ос . Предполагается, что эти же условия выполнены для всех решений 7 = 7 (7) с начальными условиями Y (10) = 70, если только

о<| \y0- Y0||< b , (2)

где /о = | ",у2 Уп " ~ вектор возмущенных начальных данных и b - некоторое

постоянное число.

Здесь и в дальнейшем под нормой вектора || • || понимается каноническая норма вектора, определенная как сумма модулей его координат.

Пусть Е п+1 - пространство точек _>',. _>' 2 ■ • • • ■ „ или 7 . Пусть Е п+1 множество

всех точек 4[,Y(t) и ¥ (! ) в пространстве £и+1,где Y (7) - всевозможные решения, получаемые из (2), при t е [0. ^

На множестве Е п+\ предполагаются выполненными условия:

1) Функция F{t,Y) определена и непрерывно дифференцируема по t на [0, со при 7 = 7(7) и всех 7=7(0 из (2).

2) На всем множестве точек из Е j выполнено условие Липшица

||f (t, Y)-F (t, Y )||< l|| Y-y||, l= const. (3)

3) Существует константа C о , такая что

||F;(i,7)||< С0, *С0 V 4,Y(t)~jßEn+i а V (t)~ßE п+1. (4)

Определение устойчивости по Ляпунову заимствуется из [2] с некоторыми упрощениями, допустимыми в рассматриваемых условиях. Решение Y — Y(t) устойчиво (справа), если для любого сколь угодно малого числа s >0 существует 8 > 0, 5 < b, такое что || Yq — Yq || < 8 влечет

lfC-Гфв Vie I о, со . Решение Y =Y(t) асимптотически устойчиво (справа), если оно

устойчиво и найдется §0>0, 8 0 < 8, такое что II 70 - 70 II < 8 0 влечет lim tf С;г Y 0 •

II II КО

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

Предложенные условия устойчивости строятся на основе разностных методов приближенного решения ОДУ, в частности, ниже приводится пример использования для этой цели метода Эйлера-Коши. Существенной особенностью при этом является выбор шага численного интегрирования h : предполагается, что для каждого произвольно фиксированного / е [0, со каково бы ни было /=0,1,...,

t = tt+и h= . . , (5)

7 + 1

при этом изменение переменной t рассматривается как изменение правой границы промежутка |о , t _. Иными словами, h = h(i) на любом промежутке |0. / _. но шаг остается равномерным внутри промежутка. Если /+1- номер заключительного шага на |0. / _. то при любом к . 0<k<i, tk+1=tk+h.

В [1] на основе тейлоровских разложений невозмущенного решения Y — Y(t) задачи Коши (1) и возмущенного решения Y = Y (/), с начальными условиями Г (/ 0) = Y(], показано, что для

разности между соответственными компонентами возмущенного и невозмущенного решений для любого t из (5) имеет место равенство

]1 " У]1 +1 = П ^ + Л ^ 0 " У] 0J itl f1 + hdJ Л qjk+qjI +1 ' (6) t=0 ^ 2 ^ к=11=к ^ 2 J

_ 2

где величины qjk имеют порядок малости O(h ) и

, _fjitJt " f{ .У .,.У •/;/•( .У " f{ .,.)' •/;/•< .)' ^

"¡I--~-; (')

УJt-Уit

7=1,2, ...,п.

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

' ' ( h Л

В [1] доказывается, что ^ [ | 1 Н— d ,( \qjk+qJj+l—>0 при h —> 0, что равносильно

к=1е=к ^ 2 '

/ —>• оо; для любого t из (5). Предельный переход в равенстве (6) влечет

yJ(t')-yj(t) = НтР.г ф 0-yj0 ^ j = 1,2,..„и, (8)

г—>оо

где

pj =U l+hdJА, (9)

о ^ 2 J

при этом условия существования и единственности решения задачи Коши всегда обеспечивают

У j (О ^ У j (О Vie [0, со], j = 1, 2, ..., п.

Соотношения (6)-(9) получаются на основе метода Эйлера-Коши с учетом остаточных членов тейлоровских разложений решений. Аналогично, в случае метода Эйлера, равенство (7) принимает вид [3]

~ •

У Л-У Л

Пусть Рг = ^ ,.I'2 ......Р„ j , где I',, из (9), / = 1, 2, ...,п. Непосредственно на основании (8) формулируются условия устойчивости в виде следующей теоремы.

Теорема 1. В условиях 1) - 3) для устойчивости по Ляпунову решения задачи (1) необходимо и достаточно, чтобы существовало 8 > 0, 8 < b, такое, что одновременно для всех решений

Y = Y (t) при ограничении 0 < y0 — y0 < 5 выполняется условие

lim /',

<С , С = const We [о,«С (10)

Для асимптотической устойчивости необходимо и достаточно, чтобы выполнялось (10) и нашлось 80>0, S0< 8, такое, что 0 < 11Y — Y0 || ^ 8 0 влечет

lim lim Pi = б. (11)

i—>со i—>00

Условие (10) означает равномерную ограниченность бесконечных произведений lim Pi в

2-» 00

случае устойчивости, при этом t, i и h связаны соотношениями (5). Условие (11) означает стремление к нулю этих произведений при t —> <х> в случае асимптотической устойчивости.

Условия (10), (11) реализуются программно путем циклического накопления частичных

произведении

П 1+2 ">'

1=о^ 2

j = 1, 2, ...,п, поведение этих произведений будет определять

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

На практике бесконечные произведения не могут быть вычислены точно. Моделирующая их поведение программа с необходимостью остановится на их приближении конечным числом сомножителей. Возникает вопрос, как такое приближение отразится на достоверности оценки устойчивости. Поэтому необходимо исследовать поведение левых частей соотношений (10), (11) в зависимости от количества сомножителей. Решение этого вопроса дано в [1]. Утверждение представленной выше теоремы переносится на случай частичных произведений Р{ с некоторыми оговорками, в частности, в формулировку теоремы добавляется, помимо требования существования 8 > 0, требование существования номера / 0 (шага А ()). начиная с которого выполняется условие

||Рг||<С, С = сопэ!. VI > 10 а V / е [0,оо. Из предельного перехода в неравенстве ||Рг||<С с

учетом (10) следует, что это условие достаточно для устойчивости решения задачи (1). В [1] при естественных ограничениях показано, что замена предельных значений в выражении условий (10),

147

(11) на их конечные приближения сохраняет достоверность оценки устойчивости.

В [1] исследована зависимость сконструированных условий устойчивости от погрешности разностного метода, на основе которого эти условия построены, при этом рассматривались методы Эйлера-Коши, Эйлера, Рунге-Кутта и Адамса. Для всех этих методов обоснована возможность компьютерного моделирования устойчивости при помощи циклического накопления произведений 1'п = [ 1 + — с! ,, ], приближаемых по разностной схеме. е=о ^ ^ )

Описание системы Лоренца. Данная система Лоренца создавалась с целью упрощенного описания явления двумерной атмосферной конвекции для решения проблемы долгосрочного прогноза погоды и впоследствии стала первой моделью детерминированного хаоса [4]. Как отмечается в [5], система Лоренца является самым упоминаемым математическим объектом за последние двадцать лет.

Используя ряд упрощающих предположений, Лоренц из уравнения Навье-Стокса получил систему ОДУ

ёх ^ а г ё у

—— = —Х2 + ГХ—у. (12)

а 2 ,

-= ху~01,

аг

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

Кроме того, модель Лоренца полностью совпадает с некоторыми моделями лазера, моделью механической системы - гиростата и моделями многих других нелинейных систем (модели «типа Лоренца») [4].

Система Лоренца диссипативна и имеет три равновесных состояния (три стационарных решения) [4]:

1) х — у - г - 0 ;

2) х = у = у/ Ь ^ -1 ^ г = г -1;

3) х = у = Ь ^ -1 г = г -1.

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

1) При г < 1 имеется только одно стационарное решение х = у = г = 0 и оно устойчиво. Это состояние чистой теплопроводности без конвекции. Устойчивость этого состояния зависит только от значения г .

2) При г > 1 нулевое решение теряет устойчивость, но появляются два новых стационарных решения. Физически эти решения соответствуют возникновению конвекции, а каждое из них связано с одним из двух возможных направлений конвективных валов. Доказано, что новая пара

ст^ + 6 + 3

точек теряет устойчивость при г =-=". Для упрощения анализа значения ст и о как

и-Ь-1

g

правило фиксируются: с = 10, b = — - так называемые «классические» значения параметров.

Тогда при 1 < г < 24,74 точки (j b (i b Ц; — \ -1 , ( у] b 4г -\ yj b Ц' — \ -1

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

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

4) Исследования различных авторов свидетельствуют, что при дальнейшем увеличении управляющего параметра r появляются чередующиеся режимы турбулентности и периодического поведения, т.е. возникают т.н. «окна периодичности». Имеются основания считать, что таких окон бесконечно много [4].

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

х = у = -J b ^ — \ 2 z = r — 1 системы Лоренца на основе программных моделей, описание которых дано в [1]. Программные модели представляют собой стандартные подпрограммы на языке программирования Object Pascal системы Delphi. Само моделирование основано на вычислении нормы бесконечного произведения из (10), (11), получающегося путем преобразования разностного метода решения ОДУ в форму бесконечного произведения.

Возможность аппроксимации бесконечных произведений частичными доказана в [1], где также обоснована возможность использования разностных методов Эйлера, Эйлера-Коши, Рунге-Кутта и Адамса для вычислений значений невозмущенного и возмущенного решений, входящих в величину из (7).

Непосредственно ниже все результаты численных экспериментов даны при использовании

метода Рунге-Кутта 4-го порядка и выборе шага h = 10 5. Графическое моделирование осуществлено с помощью системы MathCAD также с использованием метода Рунге-Кутта 4-го порядка с

Таблица 1

фиксированным шагом h = 10 5.

Результаты экспериментальной оценки асимптотической устойчивости состояния = -\] Ъ ^ — 1 ^ г^ = Г — 1 системы Лоренца

при значениях параметров а = 10, Г = 24, Ъ — 2,666667

t Норма произведения Значения Значения 7» ~ Значения z, ~

1.00 3.0E+0000 7.8E+0000 7.8E+0000 7.8E+0000 7.8E+0000 2.3E+0001 2.3E+0001

10.00 2.2E+0000 7.8E+0000 7.8E+0000 7.8E+0000 7.8E+0000 2.3E+0001 2.3E+0001

100.00 3.2E-0001 7.8E+0000 7.8E+0000 7.8E+0000 7.8E+0000 2.3E+0001 2.3E+0001

600.00 3.3E-0006 7.8E+0000 7.8E+0000 7.8E+0000 7.8E+0000 2.3E+0001 2.3E+0001

При г = 24 наблюдается монотонное убывание нормы вектора Р^ к нулю, что в соответствие с теоремой 1 условием (11) означает асимптотическую устойчивость рассматриваемого состояния, это соответствует известным результатам моделирования поведения системы Лоренца при выбранных значениях параметрах.

Ниже представлены графические интерпретации поведения системы Лоренца в данном случае полученные в МаШСАБ.

Рис. 1

График поведения возмущенного решения у системы Лоренца при X о = у о = у] Ь ^ — 1 ^ г о = Г — 1 для значений параметров

с = 10, г = 24, Ъ = 2,666667

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

Рис. 2

Фазовый портрет системы Лоренца при X 0 = у 0 = у] Ь ^ — 1 ^ г 0 = Г — 1 для значений параметров (7 = 10, Г = 24, Ь — 2,666667

Из устойчивости по Ляпунову следует орбитальная устойчивость [2], рис. 2 иллюстрирует орбитально устойчивые фазовые траектории системы Лоренца, стягивающиеся к рассматриваемой точке.

Таблица 2

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

= У 0 = л/Й^Г; 2 0 = г ~ 1 системы Лоренца при значениях параметров С = 10 , Г — 24,74, Ь = 2,666667

г Норма произведения Значения Значения у» у Значения г, у

100.00 2.5Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

600.00 2.6Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

900.00 2.7Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

1000.00 2.3Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

При данном значении управляющего параметра г происходит качественное изменение поведения системы, а именно, свойство асимптотической устойчивости изменяется на свойство неустойчивости, при этом для значения 24,74 наблюдается переходное состояние, которое определяется как просто устойчивость. Данное известное обстоятельство полностью совпадает с полученными результатами моделирования, а именно, норма произведения || 1| не превосходит величины

2,7 на всем промежутке [). 1000 _, что в соответствие с теоремой 1 и условием (10) означает устойчивость.

Ниже приводимый рис. 3 иллюстрирует поведение возмущенного решения у в данном случае. Амплитуда колебаний этого решения находится в пределах от 7,95645 до 7,95665 и в отличие от предыдущего случая с течением времени не затухает, что соответствует поведению устойчивых решений. Графики поведения у и у аналогичны.

Рис. 3

График поведения возмущенного решения X системы Лоренца

при Х0 =Уо = — — Г — 1 для значений параметров <Т = 10 , Г — 24,74, Ъ — 2,666667

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

Рис. 4

Фазовый портрет системы Лоренца при X 0 = о = л/ ^ ~~ ^ -> 2 0 = ^ — 1 для значений параметров (7 = 10, Г = 24,74, Ь = 2,666667

Результаты экспериментальной оценки неустойчивости состояния Х0 = о = л/ ^ ~~ 1 ^ 20 = — 1 системы Лоренца при значениях параметров (7 = 10, Г = 25 , Ь — 2,666667

Таблица 3

г Норма произведения Значения х, у Значения у. у Значения г, у

100.00 5.6Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

200.00 1.0Е+0001 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

500.00 1.3Е+0002 8.0Е+0000 8.0Е+0000 8.0Е+0000 8.0Е+0000 2.4Е+0001 2.4Е+0001

1000.00 6.4Е+0003 8.0Е+0000 7.9Е+0000 8.0Е+0000 7.6Е+0000 2.4Е+0001 2.4Е+0001

Норма произведения монотонно возрастает, что означает нарушение условия (10) и,

следовательно, по теореме 1 имеет место неустойчивость.

Аналогично двум предыдущим случая ниже приводятся рис. 5, 6, дающие графическую интерпретация данного случая.

Рис. 5

График поведения возмущенного решения X системы Лоренца при Х0 = у о = у] Ь ^ — 1 ^ г0 = Г — 1 для значений параметров (Т = 10, Г = 25 , Ь = 2,666667

На рис. 5 амплитуда колебаний возмущенного решения X с течением времени возрастает, что является свидетельством неустойчивости. Графики поведения ~ , ~ аналогичны.

8.0002 8

7.9998 Ф.9996

7.9997 7.9998 7.9999 8 8.0001 8.0002

Рис. 6

Фазовый портрет системы Лоренца при Х0 = у о = у/ Ь — 1 3 = Г — 1 для значений параметров О = 10, Г = 25 , Ь — Т.,666661

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

Таким образом, практическое использование предложенных условий устойчивости показало их пригодность для компьютерного анализа устойчивости даже в случае такого нетривиального объекта как система Лоренца. Аналогичные результаты практического применения условий имели место для многочисленных систем нелинейных ОДУ [1], с которыми проводились численные эксперименты.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Катрич С.А. Разработка и исследование схем программного моделирования устойчивости решений нелинейных дифференциальных уравнений на основе разностных методов: дис. ... канд. техн. наук. Таганрог: Изд-во ТРТУ, 2006. 217 с.

2. Демидович Б.П. Лекции по математической теории устойчивости. М.: Изд-во Моск. ун-та, 1998. 480 с.

3. Ромм Я.Е. Мультипликативные критерии устойчивости на основе разностных решений обыкновенных дифференциальных уравнений // Кибернетика и системный анализ. Киев, 2006. № 1. С. 128-144.

4. Синергетика: процессы самоорганизации и управления: учеб. пособие / под общей редакцией А.А. Колесникова. В 2-х частях. Таганрог: Изд-во ТРТУ, 2004. 360 с. Ч. I.

5. Зубов И.В. Методы анализа динамики управляемых систем. М.: Физматлит, 2003. 224 с.

Я.Е. Ромм, О.И. Деева

ИДЕНТИФИКАЦИЯ КОНТУРНО ПРЕДСТАВЛЕННЫХ ИЗОБРАЖЕНИЙ НА

ОСНОВЕ СОРТИРОВКИ И МАТРИЧНО-АЛГЕБРАИЧЕСКИХ ПРЕОБРАЗОВАНИЙ

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

1. Выделение контура фигуры и ее переход к базовому положению в декартовой системе координат

Данная часть работы непосредственно заимствуется из [3, 4]. Задача выделения контуров состоит в построении границ объектов и очертаний связных областей. Контур отдельного объекта на бинарном изображении может быть получен при помощи следящей пары точек [1]: граничной точки объекта и смежной с ней точки фона. Суть данного способа заключается в следующем. Для следящей пары задается направление от точки объекта к точке фона. Процесс прослеживания состоит в последовательном перемещении одного конца следящей пары в новую точку, лежащую от пары слева. Это обеспечивает обход контуров таким образом, что объект оказывается слева от границы, а фон - справа. На каждом шаге прослеживания анализируется одна новая пробная точка, смежная с обеими точками пары. Пробная точка замещает соответствующую точку следящей пары. Пусть К/),С/ - следящая пара точек, р= 4рх,ру - граничная точка объекта, а

q = С/,,Чу ~ точка фона. Пробная точка г = 4х,гу ищется следующим образом: если точки следящей пары лежат в одном столбце или в одной строке матрицы /, то новая пробная точка выбирается с координатами гх = с/х + ру — с/у. ¡\. = с/у + с/х — рх: в противном случае пробная точка имеет координаты гх = 4рх+qx+ р - q 12, г = ip + q - рх+qx J2.