Якісний аналіз динамічних систем. Побудова фазових портретів ДС

КІНЕТИКА БІОЛОГІЧНИХ ПРОЦЕСІВ

Як можна описати динаміку біологічних систем? У кожний момент часу біологічна система має набір деяких характеристик. Наприклад, спостерігаючи за популяцією якогось виду, можна реєструвати її чисельність, площу займаної території, кількість доступного харчування, температуру довкілляі т. д. Перебіг хімічної реакції можна характеризувати концентраціями речовин, що беруть участь, тиском, температурою, рівнем кислотності середовища. Сукупність значень всіх параметрів, які дослідник вибрав для опису системи, є станом системи в кожний момент часу. При створенні моделі у зазначеній сукупності виділяють змінні та параметри. Змінні – це величини, зміни яких насамперед цікавить дослідника, параметри – умови «довкілля». Саме обраних змінних становлять рівняння, що відбивають закономірності зміни системи у часі. Наприклад, при створенні моделі зростання культури мікроорганізмів, як змінна зазвичай виступає її чисельність, а як параметр – швидкість розмноження. Можливо, суттєвою виявиться температура, при якій відбувається зростання, тоді цей показник також включається до моделі як параметр. А якщо, наприклад, рівень аерації завжди є достатнім і не впливає на ростові процеси, тоді його взагалі не включають у модель. Як правило, параметри залишаються незмінними під час експерименту, проте слід зазначити, що це не завжди так.

Описувати динаміку біологічної системи (тобто зміна її стану у часі) можна дискретними, і безперервними моделями. У дискретних моделях передбачається, що час є дискретну величину. Це відповідає реєстрації значень змінних через певні фіксовані інтервали часу (наприклад, раз на годину або на рік). У безперервних моделях біологічна змінна є безперервною функцією часу, що позначається, наприклад, x(t).

Часто велике значення мають початкові умовимоделі – стан досліджуваної показники початковий час, тобто. при t = 0.

При вивченні безперервної зміни певної характеристики x(t) нам може бути відома інформація про швидкість її зміни. Ця інформація може бути записана у вигляді диференціального рівняння:

Така формальна запис означає, що швидкість зміни деякої досліджуваної характеристики є функцією часу та величини цієї характеристики.

Якщо права частина диференціального рівняння виду не залежить від часу, тобто. справедливо:

то таке рівняння називається автономним(система, що описується таким рівнянням, називається автономною). Стан автономних систем у кожний момент часу характеризується однією єдиною величиною – значенням змінної xв даний момент часу t.

Задамося питанням: нехай дано диференціальне рівняння для x(t), чи можна знайти всі функції x(t), які задовольняють цьому рівнянню? Або: якщо відомо початкове значення деякої змінної (наприклад, початковий розмір популяції, концентрація речовини, електропровідність середовища тощо) і є інформація про характер зміни цієї змінної, то чи можна передбачити, яким буде її значення у наступні моменти часу? Відповідь на поставлене питання звучить наступним чином: якщо задані початкові умови при і для рівняння виконані умови теореми Коші (функція , задана в деякій області, і її похідна постійна безперервні в цій області), то є єдине рішення рівняння, що задовольняє заданим початковим умовам. (Нагадаємо, що будь-яка безперервна функція, яка задовольняє диференціальному рівнянню, називається рішенням цього рівняння.) Це означає, що ми можемо однозначно передбачати поведінку біологічної системи, якщо відомі характеристики її початкового стану, і рівняння моделі задовольняє умовам теореми Коші.

Стаціонарний стан. Стійкість

Розглянемо автономне диференціальне рівняння

У стаціонарному стані значення змінних у системі змінюються згодом, тобто швидкість зміни значень змінних дорівнює 0: . Якщо ліва частина рівняння (1.2) дорівнює нулю, то права також дорівнює нулю: . Коріння цього рівня алгебри є стаціонарними станамидиференціального рівняння (1.2).

Приклад1.1:Знайдіть стаціонарні стани рівняння.

Рішення: Перенесемо доданок, що не містить похідну, у праву частину рівності: . За визначенням у стаціонарному стані виконується рівність: . Значить, має виконуватись рівність . Вирішуємо рівняння:

,

Отже, рівняння має 3 стаціонарні стани: , .

Біологічні системи постійно відчувають різні зовнішні дії та численні флуктуації. У цьому вони (біологічні системи) мають гомеостаз, тобто. стійкі. Математичною мовою це означає, що змінні при малих відхиленнях повертаються до своїх стаціонарних значень. Чи відображатиме такий характер поведінки біологічної системи її математична модель? Чи стійкі стаціонарні стани моделі?

Стаціонарний стан є стійким, якщо при досить малому відхиленні від положення рівноваги система ніколи не піде далеко від особливої ​​точки. Стійкий стан відповідає стійкому режиму функціонування системи.

Стан рівноваги рівняння стійкий за Ляпуновим, якщо для будь-якого завжди можна знайти такий, що якщо, то для всіх.

Існує аналітичний метод дослідження сталості стаціонарного стану – метод Ляпунова. Для його обґрунтування нагадаємо формулу Тейлора.

Говорячи нестрого, формула Тейлора показує поведінку функції навколо певної точки. Нехай функція має у точці похідні всіх порядків до n-го включно. Тоді для справедлива формула Тейлора:

Відкинувши залишковий член, який представляє себе нескінченно малу вищого порядку, ніж отримаємо наближену формулу Тейлора:

Права частина наближеної формули називається багаточлен Тейлорафункції, його позначають як.

Приклад1.2:Розкладіть функцію в ряд Тейлора на околиці точки до 4 порядку.

Рішення:Запишемо ряд Тейлора до 4-го порядку загальному вигляді:

Знайдемо похідні заданої функції в точці:

,

Підставимо отримані значення у вихідну формулу:

Аналітичний метод дослідження стійкості стаціонарного стану ( метод Ляпунова) полягає у наступному. Нехай – стаціонарний стан рівняння. Задамо невелике відхилення змінної xвід її стаціонарного значення: , де . Підставимо вираз для точки xу вихідне рівняння: . Ліва частина рівняння набуде вигляду: , оскільки у стаціонарному стані швидкість зміни значення зміною дорівнює нулю: . Праву частину розкладемо в ряд Тейлора в околиці стаціонарного стану, враховуючи, що залишимо лише лінійний член у правій частині рівняння:

Отримали лінеаризоване рівнянняабо рівняння першого наближення. Розмір є певна стала величина, позначимо її a: . Загальне рішення лінеаризованого рівняння має вигляд: . Це вираз описує закон, яким змінюватиметься у часі задане нами відхилення від стаціонарного стану . Відхилення згодом згасатиме, тобто. при , якщо показник ступеня експоненті буде негативним, тобто. . За визначенням стаціонарний стан буде стійким. Якщо ж , то зі збільшенням часу відхилення тільки збільшуватиметься, стаціонарний стан – нестійке. У разі коли рівняння першого наближення відповіді на питання про стійкість стаціонарного стану дати не може. Необхідно розглядати члени вищого порядку у розкладанні до ряду Тейлора.

Крім аналітичного методу дослідження сталості стаціонарного стану існує і графічний.

Приклад1.3.Нехай. Знайти стаціонарні стани рівняння та визначити їх тип стійкості за допомогою графіка функції .

Рішення:Знайдемо особливі точки:

,

,

Будуємо графік функції (рис.1.1).

Мал. 1.1. Графік функції (Приклад 1.3).

Визначимо за графіком чи стійкий кожен зі знайдених стаціонарних станів. Задамо невелике відхилення зображувальної точки від особливої ​​точки вліво: . У точці з координатою функція набуває позитивного значення: або . Остання нерівність означає, що з часом координата повинна збільшуватися, тобто точка, що зображає, повинна повернутися до точки . Тепер задаємо невелике відхилення зображувальної точки від особливої ​​точки праворуч: . У цій галузі функція зберігає позитивне значення, отже, з часом координата xтакож збільшується, тобто точка, що зображає, буде віддалятися від точки . Таким чином, мале відхилення виводять систему зі стаціонарного стану, отже, за визначенням, особлива точка нестійка. Аналогічні міркування призводять до того, що будь-яке відхилення від особливої ​​точки з часом згасає, стаціонарний стан стійкий. Відхилення зображувальної точки в будь-якому напрямку від стаціонарного стану призводить до її віддалення від точки, це нестійкий стаціонарний стан.

Рішення системи лінійних диференціальних рівнянь

Перейдемо до вивчення систем рівнянь, спочатку лінійних. У загальному вигляді систему лінійних диференціальних рівнянь можна подати у вигляді:

Аналіз системи рівнянь починається із знаходження стаціонарних станів. У систем виду (1.3) особлива точка єдина, її координати – (0,0). Виняток становить вироджений випадок, коли рівняння можна подати у вигляді:

(1.3*)

І тут всі пари , які відповідають співвідношенню , є стаціонарними точками системи (1.3*). Зокрема, точка (0,0) також стаціонарної для системи (1.3*). На фазовій площині у разі маємо пряму з коефіцієнтом нахилу , що проходить через початок координат, кожна точка якої є особливою точкою системи (1.3*) (див. таблицю1.1, пункт 6).

Основне питання, на яке має відповідати результат дослідження системи рівнянь: чи стійкий стаціонарний стан системи, і який характер має це рішення (монотонний чи немонотонний).

Загальне рішеннясистеми двох лінійних рівнянь має вигляд:

Характеристичні числаможуть бути виражені через коефіцієнти лінійних рівнянь наступним чином:

Характеристичні числа можуть бути: 1) дійсними різних знаків; 2) дійсними одного знака; нулю. Ці випадки визначають тип поведінки рішення системи звичайних диференціальних рівнянь. Відповідні фазові портрети представлені таблиці1.1.

Таблиця 1.1. Типи стаціонарних станів системи двох лінійних диференціальних рівнянь та відповідні фазові портрети. Стрілками показано напрямок руху зображувальної точки

Побудова фазових та кінетичних портретів системи двох лінійних диференціальних рівнянь

Фазовою площиноюназивається площина з осями координат, на яких відкладені значення змінних xі yкожна точка площини відповідає певному стану системи. Сукупність точок на фазовій площині, положення яких відповідає станам системи в процесі зміни в часі змінних, згідно з заданими рівняннями системи, що досліджується, називається фазовою траєкторією. Сукупність фазових траєкторій за різних початкових значеннях змінних дає портрет системи. Побудова фазового портретадозволяє зробити висновки про характер змін змінних xі yбез знання аналітичних рішень вихідної системи рівнянь.

Розглянемо систему лінійних диференціальних рівнянь:

Побудова фазового портрета починаємо з побудови головних ізоклін(Ізокліна - лінія, на всьому протязі якої нахил фазової кривої (траєкторії), що визначається рівнянням, зберігає постійне значення). Для системи двох лінійних диференціальних рівнянь – це прямі, які проходять через початок координат. Рівняння ізокліни горизонтальних дотичних: . Рівняння ізокліни вертикальних дотичних: . Для подальшої побудови фазового портрета корисно побудувати ізокліну дотичних, що проходять під кутом. Для знаходження відповідного рівняння ізокліни необхідно вирішити рівняння . Можна знаходити і ізокліни щодо інших кутів, користуючись приблизними значеннями тангенсів кутів. У побудові фазового портрета також може допомогти у відповідь питання, під яким кутом фазові траєкторії повинні перетинати координатні осі. Для цього в рівняння ізокліни підставляємо відповідні рівності (для визначення кута перетину з віссю OY) та (для визначення кута перетину з віссю OХ).

Приклад1.4.Визначте тип особливої ​​точки системи лінійних рівнянь:

Побудуйте фазовий та кінетичний портрет системи.

Рішення:Координати особливої ​​точки – (0,0). Коефіцієнти лінійних рівнянь дорівнюють: , , , . Визначимо тип стаціонарного стану (див. розділ про характеристичні числа):

Отже, характеристичні коріння є уявними: , отже, особлива точка аналізованої лінійної системи має тип центр (рис. 1.2а).

Рівняння ізокліни горизонтальних дотичних: , рівняння ізокліни вертикальних дотичних: . Під кутом у 45° траєкторії системи перетинають пряму .

Після побудови фазового портрета необхідно визначити напрямок руху по знайденим траєкторіям. Це можна зробити в такий спосіб. Візьмемо довільну точку на будь-якій траєкторії. Наприклад, на ізоклині горизонтальних дотичних (1,1). Підставимо координати цієї точки до системи рівнянь. Отримаємо вирази для швидкостей зміни змінних x,yу цій точці:

Отримані значення показують, що швидкість зміни змінної x– негативна, тобто її значення має зменшуватись, а змінна yне змінюється. Зазначаємо отриманий напрямок стрілкою. Таким чином, у прикладі рух по фазових траєкторіях спрямовано проти годинникової стрілки. Підставляючи в систему координати різних точок можна отримати карту напрямків швидкостей, так зване векторне поле.

Рис. 1.2. Фазовий (а) та кінетичний (б) портрет системи, приклад 1.4

Зазначимо, що на ізоклині горизонтальних дотичних змінна yдосягає свого максимального чи мінімального значення даної траєкторії. Навпаки, на ізоклині вертикальних дотичних, свого максимального за модулем значення для обраної траєкторії досягає змінна x.

Побудувати кінетичний портрет системи означає побудувати графіки залежності величин змінних x,yвід часу. За фазовим портретом можна побудувати кінетичний і навпаки. Однією фазовою траєкторією відповідає одна пара кінетичних кривих. Виберемо на фазовому портреті довільну точку довільної фазової траєкторії. Це початкова точка, що відповідає моменту часу. Залежно від напрямку руху в системі значення змінних x,yабо зменшуються або збільшуються. Нехай координати початкової точки – (1,1). Згідно з побудованим фазовим портретом, стартуючи з цієї точки, ми повинні рухатися проти годинникової стрілки, координати xі yпри цьому зменшуватимуться. З часом координата xпроходить через 0, значення yу своїй залишається позитивним. Далі координати xі yпродовжують зменшуватися, координата yпроходить через 0 (значення xу своїй негативно). Величина xдосягає мінімального значення на ізоклині вертикальних дотичних, потім починає збільшуватися. Величина yсвого мінімального значення досягає на ізоклині горизонтальних дотичних (значення xу цей час негативно). Далі і величина x, і величина yзбільшуються, повертаючись до початкових значень (рис. 1.2б).

Дослідження стійкості стаціонарних станів нелінійних систем другого порядку

Нехай біологічна система описується системою двох автономних диференціальних рівнянь другого порядку загального виду:

Стаціонарні значення змінних системи визначаються алгебраїчних рівнянь:

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

Поза

Маємо, , , особлива точка груба. Характеристичні коріння системи першого наближення рівні , обидва дійсні та негативні, отже, в околиці нульової особливої ​​точки поведінка фазових траєкторій системи буде відповідати типу стійкий вузол.

Метою дослідження є розробка орієнтованого на застосування суперкомп'ютерів логічного методу (методу булевих обмежень) та сервіс-орієнтованої технології створення та застосування комп'ютерної системи для якісного дослідження динаміки поведінки траєкторій автономних двійкових динамічних систем на кінцевому інтервалі часу. Актуальність теми підтверджується безперервно зростаючим спектром додатків двійкових моделей у наукових та прикладних дослідженнях, а також необхідністю якісного аналізу таких моделей з великою розмірністю вектора станів. Наведено математичну модель автономної двійкової системи на кінцевому інтервалі часу та еквівалентне цій системі булеве рівняння. Специфікацію динамічної властивості пропонується записувати мовою логіки предикатів з використанням обмежених кванторів існування та загальності. Отримано булеви рівняння пошуку рівноважних станів та циклів двійкової системи та умови їх ізольованості. Специфіковані основні властивості типу досяжності (досяжність, безпека, одночасна досяжність, досяжність при фазових обмеженнях, тяжіння, зв'язковість, тотальна досяжність). Для кожної властивості побудована його модель у вигляді бульового обмеження (бульова рівняння або квантифікованої булевої формули), що задовольняє логічної специфікації властивості та рівнянь динаміки системи. Таким чином, перевірка здійсненності різноманітних властивостей поведінки траєкторій автономних двійкових динамічних систем на кінцевому інтервалі часу зведена до завдання здійсненності булевих обмежень з використанням сучасних SAT та TQBF вирішувачів. Наведено демонстраційний приклад використання цієї технології для перевірки здійсненності деяких із наведених раніше властивостей. У висновку перераховано основні переваги методу булевих обмежень, особливості його програмної реалізації в рамках сервіс-орієнтованого підходу та позначені напрями подальшого розвитку методу для інших класів двійкових динамічних систем.

двійкова динамічна система

динамічна властивість

якісний аналіз

булеви обмеження

завдання булевої здійсненності

1. Biere A., Ganesh V., Grohe M., Nordstrom J., Williams R. Theory and Practice of SAT Solving. Dagstuhl Reports. 2015. vol. 5. no. 4. Р. 98-122.

2. Marin P., Pulina L., Giunchiglia E., Narizzano M., Tacchella A. Twelve роки з QBF Evaluations: QSAT є PSPACE-Hard and It Shows. Fundam. Inform. 2016. vol. 149. Р. 133-58.

3. Бохман Д., Постхоф Х. Двійкові динамічні системи. М.: Вища школа, 1986. 400 с.

4. Маслов С.Ю. Теорія дедуктивних систем та її застосування. М.: Радіо та зв'язок, 1986. 133 с.

5. Jhala R., Majumdar R. Software model checking. ACM Computing Surveys. 2009. vol. 41. no. 4. Р. 21:1–21:54.

6. Васильєв С.М. Метод редукції та якісний аналіз динамічних систем. I-II // Вісті РАН. Теорія та системи управління. 2006. № 1. С. 21-29. № 2. С. 5-17.

7. DIMACS format [Електронний ресурс]. Режим доступу: http://www.cs.utexas.edu/users/moore/acl2/manuals/current/manual/index-seo.php/SATLINK____DIMACS (дата звернення: 24.07.2018).

8. QDIMACS standard [Електронний ресурс]. Режим доступу: http://qbflib.org/qdimacs.html (дата звернення: 24.07.2018).

9. Delgado-Eckert E., Reger J., Schmidt K. Discrete Time Systems з Event-Based Dynamics: Recent Developments in Analysis and Synthesis Methods. Mario Alberto Jordan (Ed). Discrete Time Systems. InTech. 2011. Р. 447-476.

10. Васильєв С.М. Досяжність і зв'язність у автоматній мережі із загальним правилом перемикання станів // Диференціальні рівняння. 2002. Т. 38. № 11. С. 1533-1539.

11. Бичков І.В., Опарін Г.А., Богданова В.Г., Горський С.А., Пашинін А.А. Мультіагентна технологія автоматизації паралельного розв'язання булевих рівнянь у розподіленому обчислювальному середовищі // Обчислювальні технології. 2016. Т. 21. № 3. С. 5-17.

12. Lonsing F., Biere A. DepQBF. A Dependency-Aware QBF solver. Journal on Satisfiability. Boolean Modeling and Computation. 2010. vol. 9. Р. 71-76.

13. Oparin G.A., Bogdanova V.G., Pashinin A.A., Gorsky S.A. Distributed Solvers of Applied Problems Базовані на Microservices and Agent Networks. Proc. З 41th Intern. Конвенція про інформаційну і комунікаційну технологію, електронні та мікроелектроніки (MIPRO-2018). Р. 1643-1648.

14. Богданова В.Г., Горський С.А. Scalable Parallel Solver of Boolean Satisfiability Problems. Proc. З 41th Intern. Convention on Information and Communication Technology. Electronics and Microelectronics (MIPRO-2018). Р. 244-249.

15. Бичков I.V., Oparin G.A., Bogdanova V.G., Pashinin A.A. The Applied Problems Solving Technology Based on Distributed Computational Subject Model: a Decentralized Approach // Паралельні обчислювальні технології XII міжнародна конференція, ПВТ'2018, м. Ростов-на-Дону, 2–6 квітня 2018 р. Короткі статті та описи Челябінськ: Видавничий центр ЮУрДУ, 2018. C. 34-48.

Спектр додатків двійкових динамічних моделей надзвичайно широкий і з кожним роком кількість об'єктів та завдань, де потрібне їх використання, лише зростає. Класичним прикладом є двійковий синхронний автомат, що є моделлю багатьох дискретних пристроїв у системах управління, обчислювальної техніки, телемеханіки. До сучасних додатків двійкових динамічних моделей ставляться завдання біоінформатики, економіки, соціології та інших, здавалося б далеких від застосування двозначних змінних, областей. У зв'язку з цим значною мірою підвищується актуальність розробки нових та вдосконалення існуючих методів якісного аналізу поведінки траєкторій двійкових динамічних систем (ДДС).

Як відомо, метою якісного аналізу динамічної системи (не тільки двійкової) є отримання позитивної чи негативної відповіді на запитання: Чи виконується в заданій системі потрібна динамічна властивість? Перефразуємо це питання так: Чи задовольняє поведінка траєкторій динамічної системи деякої сукупності обмежень, що характеризують властивість? Далі будемо використовувати саме це трактування мети якісного аналізу динамічних властивостей системи.

Для ДДС, функціонування яких розглядається на кінцевому інтервалі часу, такі обмеження є булевими та записуються мовою булевих рівнянь або булевих формул з кванторами. Перший тип обмежень призводить до необхідності вирішення SAT-завдання (завдання булевої здійсненності); другий тип обмежень пов'язаний з розв'язанням задачі TQBF (перевірки істинності квантифікованих булевих формул). Перше завдання є типовим представником класу складності NP, а друге завдання – класу складності PSPACE. Як відомо, PSPACE-повнота дискретної задачі дає сильніше свідчення про її труднорозв'язність, ніж NP-повнота. З огляду на це зведення завдання якісного аналізу ДДС до SAT-задачі більш переважно, ніж зведення до завдання TQBF. У загальному випадку дослідження не кожної властивості ДДС можна подати мовою булевих рівнянь.

Теоретична можливість використання булевих обмежень (а саме булевих рівнянь) у якісному аналізі ДДС вперше була продемонстрована в роботі. Слід, проте, відзначити, що застосування цього підходу практично тоді стримувалося відсутністю ефективних алгоритмів і програм розв'язання булевих рівнянь (особливо з великою кількістю невідомих змінних), дозволяють значною мірою скоротити простір пошуку. В останнє десятиліття в результаті інтенсивних досліджень у цій галузі з'явилася достатня кількість різноманітних ефективних вирішувачів булевих рівнянь (SAT-рішителів), що використовують сучасні досягнення (нові евристики, швидкі структури даних, паралельні обчислення та ін.) у вирішенні задачі булевої здійсненності. Аналогічні процеси (але з деяким запізненням) спостерігаються і в галузі створення дедалі ефективніших алгоритмів та програм вирішення задачі TQBF. Таким чином, до цього часу є всі необхідні передумови систематичного розвитку методу булевих обмежень у якісному аналізі ДДС, його програмної реалізації та застосування у вирішенні наукових та прикладних завдань.

Крім методу булевих обмежень, до ДДС застосовуються й інші методи якісного аналізу, до яких відносяться дедуктивний аналіз, model checking та метод редукції. Кожен із цих методів (включаючи і метод булевих обмежень) має свої обмеження, переваги та недоліки. Загальним недоліком є ​​те, що всі методи мають перебірний характер і проблема скорочення перебору є фундаментальною для цих методів.

Важливість дедуктивного аналізу , що передбачає застосування аксіом і правил виведення на підтвердження правильності функціонування системи, визнається широким колом фахівців, але це трудомісткий і тому метод, що рідко застосовується. У методі model checking як мови специфікації необхідної властивості використовується мова темпоральних логік, яка не звична для фахівців з автоматної динаміки. Метод редукції пов'язаний із побудовою спрощеної (у певному сенсі) моделі вихідної системи, дослідженні її властивостей та умов переносимості цих властивостей у вихідну складну систему. Умови переносимості властивостей мають лише достатній характер. Простота ідеї методу редукції у якісному аналізі ДДС стикається з проблемою вибору спрощеної системи, яка задовольняє всі умови методу.

Практичне використання методу булевих обмежень передбачає алгоритмізацію та автоматизацію наступних процесів:

1) розробка орієнтованого спеціаліста з динаміці систем логічного мови специфікації динамічних властивостей;

2) побудова моделі динамічної властивості у вигляді бульового обмеження того чи іншого типу, що задовольняє логічної специфікації властивості та рівнянь динаміки двійкової системи;

3) представлення отриманої моделі у міжнародному форматі DIMACS або QDIMACS;

4) вибір (розробка) ефективного паралельного (розподіленого) вирішувача завдання здійсненності булевих обмежень (SAT або TQBF вирішувача);

5) розробка інструментальних засобів створення програмних сервісів;

6) розробка сервісів для якісного вивчення різноманітних динамічних якостей ДДС.

Метоюданого дослідження є рішення лише перших двох завдань стосовно алгоритмізації якісних досліджень автономних (без керуючих входів) синхронних ДДС. Такі системи в англомовних публікаціях називають синхронними булевими мережами (Boolean network). Інші аспекти застосування методу булевих обмежень (у тому числі і для ДДС з входами, що управляють) є предметом наступних публікацій.

Математична модель автономної ДДС

Нехай X = Bn (B = (0, 1) - безліч двійкових векторів розмірності n (простір станів ДДС).Через t∈T = (1,...,k) позначимо дискретний час (номер такту).

Для кожного стану x0∈X, званого початковим станом, визначимо траєкторію x(t, x0) як кінцеву послідовність станів x0, x1,…, xk з множини X. Далі розглядатимемо ДДС, в яких кожна пара суміжних станів xt, x(t - 1) (t∈T) траєкторії пов'язана ставленням

xt = F(xt – 1). (1)

Тут F:X>X - векторна функція алгебри логіки, яка називається функцією переходів. Таким чином, для будь-яких x0∈X система булевих рівнянь (1) представляє модель динаміки поведінки траєкторій ДДС у просторі станів X на кінцевому інтервалі часу T = (1, 2, ..., k). Тут і далі величина k у визначенні множини T передбачається наперед заданою постійною. Таке обмеження є цілком природним. Справа в тому, що при якісному аналізі поведінки траєкторій ДДС практичний інтерес представляє питання про те, що можна сказати про здійсненність будь-якої динамічної властивості при фіксованому, не надто великому k. Вибір величини k у кожному конкретному випадку здійснюється виходячи з апріорних відомостей про тривалість протікання процесів в дискретній системі, що моделюється.

Відомо, що система булевих рівнянь (1) з початковим станом x0∈X для T = (1, 2,...,k) еквівалентна одному булевому рівнянню виду

При k = 1 (розглядаються тільки однокрокові переходи) рівняння (2) набуває вигляду

(3)

Розв'язання цього рівняння визначають спрямований граф, що складається з 2n вершин, відзначених одним із 2n станів множини X. Вершини x0 і x1 графа з'єднані дугою, спрямованою від стану x0 до стану x1. Такий граф теоретично двійкових автоматів називається діаграмою переходів. Подання поведінки ДДС як діаграми переходів дуже наочно як із побудові траєкторій, і дослідженні їх властивостей, але практично реалізовано лише невеликих розмірностей n вектора стану x∈X.

Мовні засоби специфікації динамічних властивостей

Найбільш зручно задавати специфікацію динамічної властивості мовою формальної логіки. Дотримуючись роботи , позначимо через X0∈X, X1∈X, X*∈X - множини початкових, допустимих і цільових станів.

p align="justify"> Основними синтаксичними елементами логічної формули динамічної властивості є: 1) предметні змінні (компоненти векторів x0, x1, ..., xk, час t); 2) обмежені квантори існування та загальності; 3) логічні зв'язки v, &; останні формули. Заключна формула представляє твердження про належність деяких станів множини траєкторій x(t, x0) (x0∈X0) оцінним множинам X* і X1.

Слід зазначити, що використання обмежених кванторів існування та загальності забезпечує звичний для фахівця з динаміки вид запису динамічної властивості. У процесі побудови булевої моделі властивості для системи (1) обмежені квантори замінюються на звичайні згідно з такими визначеннями:

де A(y) - предикат, що обмежує значення змінної y.

В силу кінцівки області зміни змінної t, обмежені квантори існування та загальності по цій змінній замінюються на рівносильні формули, що не містять кванторів.

Надалі припускатимемо, що елементи множин X0, X1, X* визначаються відповідно нулями наступних булевих рівнянь

або характеристичними функціями цих множин.

З урахуванням обмеження на початкові стани G0(x) = 0 поряд із рівняннями (2, 3) для скорочення запису будемо використовувати наступні булеві рівняння:

(4)

Попередній якісний аналіз автономних ДДС

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

Стан x1 в (3) називатимемо послідовником стану x0, а x0 - попередником стану x1. В автономній ДДС кожен стан має лише одного послідовника, а число попередників даного стану може змінюватися від нуля до 2n - 1. Усі безпосередні попередники x0 стану s∈X є нулями бульового рівняння

Якщо рівняння (6) немає рішень, то попередники стану s відсутні.

Рівноважні стани (якщо вони існують) є рішеннями бульова рівняння

Траєкторія x0, x1, ..., xk називається циклом довжини k, якщо стани x0, x1, ..., xk-1 попарно відмінні один від одного і xk = x0. Циклічна послідовність довжини k (якщо вона існує) є рішенням бульова рівняння

де = 0 ( ) - умови попарної відмінності множини станів C циклу довжини k. Якщо жоден із станів циклу не має попередників, що не належать множині C, то такий цикл називається ізольованим. Нехай елементи s множини C визначаються рішенням бульова рівняння Gc(s) = 0. Тоді нескладно показати, що умова ізольованості циклу еквівалентна відсутності нулів у наступного бульова рівняння:

Рішення рівняння (7) (якщо вони існують) визначають стани циклу, що мають попередників, що не належать множині C.

Оскільки рівноважний стан є циклом довжини k = 1, то умова його ізольованості аналогічна умові ізольованості з k ≥ 2 за тією відмінністю, що Gc(s) має вигляд повної диз'юнкції, що визначає цей рівноважний стан.

Неізольовані рівноважні стани та цикли надалі називатимемо атракторами.

Специфікація динамічних властивостей типу досяжності

Основною властивістю ДДС, необхідність у перевірці якого найчастіше виникає практично, є традиційно досліджуване теоретично графів (у разі таким графом є ​​діаграма переходів) властивість досяжності та її різні варіації. Досяжність визначається як класичне завдання аналізу поведінки траєкторій ДДС.

Визначення цієї властивості пов'язане із завданням введених раніше множин X0, X *, X1 (що відповідають цим множинам булевих рівнянь). Передбачається, що для множин X0, X*, X1 виконується обмеження

З огляду на кінцівки множини T властивість досяжності та її варіації далі будемо розуміти як властивість практичної досяжності (досяжності за кінцеве число тактів). Розглядаються такі властивості типу досяжності:

1. Основна властивість досяжності множини X* з множини X0 формулюється наступним чином: будь-яка траєкторія, випущена з множини початкових станів X0, досягає цільової множини X*. З використанням обмежених кванторів існування та загальності, формула цієї властивості має вигляд:

2. Властивість безпеки забезпечує для будь-якої траєкторії, випущеної з X0, недосяжність множини X*:

3. Властивість одночасної досяжності. У ряді випадків може виставлятися «жорсткіша вимога», яка полягає в тому, щоб кожна траєкторія досягала цільової множини рівно за k тактів (k∈T):

4. Властивість досяжності при фазових обмеженнях:

Ця властивість гарантує, що всі траєкторії, що випускаються з множини X0, до моменту попадання в цільову множину X* знаходяться в множині X1.

5. Властивість тяжіння. Нехай X є атрактором. Тоді логічна формула якості тяжіння збігається з формулою основної якості досяжності:

тобто. для кожної траєкторії, випущеної з множини X0, існує час t∈T, починаючи з якого траєкторія не виходить за межі множини X*. Безліч X0 у разі належить частини області тяжіння множини X*(X0∈Xa, де Xа - повна область тяжіння (басейн) атрактора).

Зазначимо, що це змінні в наведених формулах властивостей фактично пов'язані, оскільки траєкторія x0, x1,…, xk повністю визначається початковим станом. Так як квантори змінної t замінюються на операції багатомісної диз'юнкції або кон'юнкції відповідних предикатів, в кожній з формул залишається єдиний обмежений квантор загальності (), що дозволяє записати умови здійсненності цих властивостей мовою булевих рівнянь (у вигляді SAT-завдання).

Наведемо дві властивості, перевірка яких призводить до необхідності розв'язання задачі TQBF.

6. Властивість зв'язності цільової множини:

тобто. існує початковий стан x0∈X0 такий, що кожен цільовий стан x*⊆X* у певний момент t∈T можна досягти, що означає існування відповідної цьому стану траєкторії, такої, що всі цільові стани x*∈X* належать даній траєкторії.

7. Властивість тотальної досяжності множини X* з X0:

тобто. кожен цільовий стан є досяжним з X0.

Перевірка здійсненності динамічних властивостей

Для властивостей (1-5) перевірка їх здійсненності зводиться до пошуку нулів булева рівняння, технологія формування якого має стандартизований характері і докладно розглядається лише основного властивості досяжності. Властивості (6, 7) призводять до завдання перевірки істинності квантифікованої булевої формули.

1. Основне властивість досяжності. Його логічна формула має вигляд

З урахуванням (4) запишемо формулу (8) у вигляді

де - характеристична функція множини станів траєкторії, випущеної з початкового стану x0∈X0. Позбавимося квантора існування в (9). Тоді матимемо

де - характеристична функція множини X *. Замінимо обмежені квантори загальності на звичайні квантори. В результаті отримаємо

Формула (10) істинна і тоді, коли тотожно істинно подкванторное вираз, тобто.

Тотожна істинність імплікації означає, що булева функція є логічним наслідком функції , тобто. будь-яка траєкторія з початковим станом x0∈X0 досягає цільової множини X*.

Виконаність тотожності (11) еквівалентна відсутності нулів у бульова рівняння

При отриманні (12) ми позбулися імплікації та замінили ϕ*(x0, x1,..., xk) на . Якщо рівняння (12) має хоча одне рішення, то властивість досяжності немає місця. Таке рішення представляє (у певному сенсі) контрприклад для властивостей, що перевіряється, і може допомогти досліднику виявити причину виникнення помилки.

Далі для стислості викладу для кожної властивості (2-4) випишемо тільки рівняння типу (12), пропонуючи читачеві самостійно відтворити необхідні міркування, близькі до тих, які дано для основної властивості досяжності.

2. Властивість безпеки

3. Властивість одночасної досяжності

4. Властивість досяжності при фазових обмеженнях

5. Властивість тяжіння. Виконаність цієї якості перевіряється у два етапи. У першому етапі з'ясовується, чи є безліч X* атрактором. Якщо відповідь позитивна, то на другому етапі перевіряється основна властивість досяжності. Якщо X* можна досягти з X0, то всі умови якості тяжіння виконані.

6. Властивість зв'язності

7. Властивість тотальної досяжності`

Для властивостей (6, 7) скалярна форма рівності двох булевих векторів xt = x* має вигляд

Продемонструємо викладену вище технологію якісного аналізу автономних ДДС з використанням методу булевих обмежень при перевірці здійсненності деяких із перерахованих вище властивостей для моделі 3.2 з роботи:

Позначимо через x0∈X = B3 початковий стан моделі (13). Нехай T = (1, 2). Випишемо потрібні для специфікації властивостей функції однокрокового та двокрокового переходів моделі (13):

(14)

де знайомий «.» позначено операцію кон'юнкції.

Для перевірки здійсненності кожної властивості задаються початкова (X0) та цільова (X*) множини, що визначаються нулями рівнянь G0(x) = 0, G*(x) = 0 або характеристичними функціями цих множин (див. п. 2). Як SAT-рішувача використовується решатель інструментального комплексу (ІЧ) РЕБУС , а TQBF-рішателя - DepQBF. Кодування змінних у булевих моделях аналізованих нижче властивостей цих вирішувачів наведено в табл. 1, булеві моделі цих властивостей у форматах DIMACS та QDIMACS розташовані в табл. 2.

Таблиця 1

Кодування змінних

Таблиця 2

Булеві моделі властивостей

1. Основне властивість досяжності (k = 2). Нехай X0 = (x∈X: x1 = 0), X * = (x∈X: x1 = 1). Початкова і цільова множини визначаються відповідно до рівнянь G0(x) = x1 = 0 і . Булеве рівняння (12) у цьому випадку набуває вигляду

де функція ϕ(x0, x1, x2) визначена (14). Решатель ІК РЕБУС видає відповідь «unsat» (рівняння немає нулів), в такий спосіб властивість досяжності X* з X0 виконується, що видно з наступної діаграми переходів, наведеної малюнку.

2. Цикли довжини k = 2. Циклічна послідовність довжини 2 (якщо вона існує) є розв'язком бульова рівняння

Функція має вигляд

Вираз R(x0, x1) при знаходженні циклу не включалося до рівняння, оскільки цикли довжини k = 1 (рівноважні стани) моделі (13) відсутні. За допомогою вирішувача ІК РЕБУС отримано дві відповіді (у вихідному форматі DIMACS): 1 2 3 4 5 -6 0 та 1 2 -3 4 5 6 0, що відповідають циклічним послідовностям (рисунок): ((1 1 1), (1 1 0)) та ((1 1 0), (1 1 1)). Багато станів обох циклів збігаються, що означає наявність у моделі (13) одного циклу довжини k = 2.

Діаграма переходів системи (13)

3. Властивість ізольованості циклу. Якщо елементи s множини станів C циклу довжини k = 2 визначаються розв'язуванням бульова рівняння Gc(s) = 0, то умова ізольованості циклу еквівалентна відсутності нулів у наступного бульова рівняння:

Оскільки C = ((1 1 1), (1 1 0)) маємо

Для цього рівняння вирішувач ІК РЕБУС знаходить два рішення: -1 2 3 4 5 -6 0 і -1 2 -3 4 5 -6 0 (у двійковому поданні згідно з кодуванням змінних у табл. 1 це пари станів (0 1 1), (1 1 0) і ((0 1 0), (1 1 0)) Таким чином, стан циклу (1 1 0) має два попередники, (0 1 1) і (0 1 0), які не належать безлічі станів циклу Це означає, що властивість ізольованості циклу не виконується, тобто цей цикл є атрактором.

4. Властивість тяжіння. Нехай X * = C є атрактором. Логічна формула властивості тяжіння збігається з формулою основної властивості досяжності

а відповідне булеве рівняння для нашого випадку має вигляд

Випишемо функції G0(x0), ϕ(x0, x1, x2) та . Функція ϕ(x0, x1, x2) наведена у (14). Для X * = C вираз дорівнює . Розглянемо два варіанти завдання множини початкових станів X0, для випадків виконання (А) і невиконання (B) властивості тяжіння за k = 2 тактів.

A. Нехай. Тоді

У цьому випадку для булевого рівняння (15) видається відповідь unsat. Властивість тяжіння для заданої множини X0 виконується.

B. Нехай. Тоді

У цьому випадку ІЧ РЕБУС для рівняння (15) знаходить рішення: 1 -2 3 4 -5 -6 -7 8 9 0, яке відповідає траєкторії ((1 0 1),(1 0 0),(0 1 1)) . Ця траєкторія з початковим станом x0 = (1 0 1) за два такти не досягає множини X * = C, що означає нездійсненність властивості тяжіння для заданого X0.

5. Властивість зв'язності. Логічна формула якості зв'язності має вигляд такого висловлювання:

Для k = 2 ϕ*(x0, x1, x2) = G0(x0)∨ϕ(x0, x1, x2), де функція ϕ(x0, x1, x2) наведена у (14). Як початковий виберемо стан (1 0 1). Тоді. Нехай цільова множина X* = ((0 1 1), (1 0 0)). У цьому випадку функція G*(x*) має вигляд

Запишемо G*(x*) у форматі КНФ:

Використовуючи закон Де-Моргана, знайдемо заперечення функції x(x0, x1, x2). Підставивши в (16) усі отримані функції та з урахуванням кодування булевих змінних (табл. 1), отримаємо булеву модель у форматі QDIMACS (табл. 2). Решатель DepQBF видає відповідь «sat», що означає істинність висловлювання (16). Властивість зв'язності заданих X0, X*, T = (1, 2) виконано.

Висновок

До основних переваг методу булевих обмежень у якісному дослідженні ДДС відносяться:

1. Звична для фахівця з автоматної динаміки логічна мова специфікації динамічної властивості за рахунок використання обмежених кванторів існування та загальності.

2. За формулою властивості та рівняннями динаміки автоматично виконується побудова відповідного бульова рівняння або квантифікованої булевої формули.

3. Досить просто автоматизується процес перетворення одержуваних булевих виразів у кон'юнктивну нормальну форму з подальшою генерацією файлу у форматах DIMAX і QDIMAX, що є вхідними для SAT-рішувачів та QBF-рішувачів.

4. Проблема скорочення перебору тією чи іншою мірою вирішується розробниками цих вирішувачів та екранована від фахівців з якісного аналізу ДДС.

5. Забезпечується можливість вирішення завдання якісного аналізу ДДС для великих розмірностей вектора стану n на досить тривалому проміжку часу T. За кількістю станів метод булевих обмежень кількісно порівнянний з методом model checking. Через те, що в останні рокиспостерігається суттєве зростання продуктивності спеціалізованих алгоритмів рішення SAT і TQBF-завдань, загальна кількість змінних у булевій моделі властивості для сучасних вирішувачів може вимірюватися тисячами.

Програмне забезпечення процесу якісного аналізу ДДС на основі методу булевих обмежень реалізується в рамках сервіс-орієнтованого підходу з використанням спеціалізованих розв'язувачів булевих рівнянь. У роботі наведено приклад реалізації методу булевих обмежень на основі сервіс-орієнтованого підходу для пошуку циклів та рівноважних станів у генних регуляторних мережах.

Слід зазначити, що метод булевих обмежень є загальним методом якісного аналізу ДДС на кінцевому інтервалі часу. Він застосовується як до автономним системам, до систем з управляючими входами, до систем із глибиною пам'яті більше одиниці, до ДДС загального вигляду, коли функція переходів нерозв'язна щодо стану xt і має вигляд F(xt, xt-1) = 0. Для ДДС з входами особливого значення набуває властивість керованості та її різні варіації. Окрім завдань аналізу ДДС метод булевих обмежень застосуємо до задач синтезу зворотнього зв'язку(статичної або динамічної, за станом або за входом), що забезпечують синтезованій системі виконання необхідної динамічної властивості.

Дослідження виконано за підтримки РФФД, проект № 18-07-00596/18.

Бібліографічне посилання

Транскрипт

1 Якісний аналіз динамічних систем Побудова фазових портретів ДС

2 Динамічна система 2 Динамічна система математичний об'єкт, відповідний реальним фізичним, хімічним, біологічним та інших. системам, еволюція у часі, що у будь-якому інтервалі часу однозначно визначається початковим станом. Таким математичним об'єктом може бути система автономних диференціальних рівнянь. Еволюцію динамічної системи можна спостерігати у просторі станів системи. Диференціальні рівняння вирішуються аналітично у явному вигляді рідко. Використання ЕОМ дає наближене рішення диференціальних рівнянь на кінцевому часовому відрізку, що дозволяє зрозуміти поведінка фазових траєкторій загалом. Тому значної ролі набувають методи якісного дослідження диференціальних рівнянь.

3 3 Відповідь питанням, які режими поведінки можуть встановлюватися у цій системі, можна з так званого фазового портрета системи сукупності всіх її траєкторій, зображених у просторі фазових змінних (фазовому просторі). Серед цих траєкторій є кілька основних, які визначають якісні властивості системи. До них відносяться насамперед точки рівноваги, що відповідають стаціонарним режимам системи, і замкнуті траєкторії (граничні цикли), що відповідають режимам періодичних коливань. Чи буде режим стійкий чи ні, можна судити з поведінки сусідніх траєкторій: стійка рівновага або цикл притягує всі близькі траєкторії, нестійке відштовхує хоча б деякі з них. Таким чином, «фазова площина, розбита на траєкторії, дає легко доступний для огляду «портрет» динамічної системи, вона дає можливість відразу, одним поглядом охопити всю сукупність рухів, що можуть виникнути за всіляких початкових умов». (А.А. Андронов, А.А. Вітт, С.Е. Хайкін. Теорія коливань)

4 Частина 1 Якісний аналіз лінійних динамічних систем

5 5 Лінійна автономна динамічна система Розглянемо лінійну однорідну систему з постійними коефіцієнтами: (1) dx ax by dt dy cx dy. dt Координатну площину xoy називають її фазовою площиною. Через будь-яку точку площини проходить одна і лише одна фазова крива (траєкторія). У системі (1) можливі три типи фазових траєкторій: точка, замкнута крива, незамкнена крива. Крапка на фазовій площині відповідає стаціонарному рішенню (положенню рівноваги, точці спокою) системи (1), замкнута крива періодичному рішенню, а незамкнена неперіодичному.

6 Положення рівноваги ДС 6 Положення рівноваги системи (1) знайдемо, вирішуючи систему: (2) ax by 0, cx dy 0. Система (1) має єдине нульове положення рівноваги, якщо визначник матриці системи: det ab A ad cb 0. cd Якщо ж det A = 0, то, крім нульового положення рівноваги, є й інші, тому що в цьому випадку система (2) має безліч рішень. Якісна поведінка фазових траєкторій (тип положення рівноваги) визначається власними числами матриці системи.

7 Класифікація точок спокою 7 Власні числа матриці системи знайдемо, вирішуючи рівняння: (3) 2 λ (ad)λ ad bc 0. Зауважимо, що a + d = tr A (слід матриці) та ad bc = det A. Класифікація точок спокою у випадку, коли det A 0 наведена в таблиці: Коріння рівняння (3) 1, 2 - речові, одного знака (1 2 > 0) 1, 2 - речові, різного знака (1 2< 0) 1, 2 - комплексные, Re 1 = Re 2 0 1, 2 - комплексные, Re 1 = Re 2 = 0 Тип точки покоя Узел Седло Фокус Центр

8 Стійкість точок спокою 8 Власні значення матриці системи (1) однозначно визначають характер стійкості положень рівноваги: ​​Умова на речову частину коренів рівняння (3) 1. Якщо речові частини всіх коренів рівняння (3) негативні, то точка спокою системи (1) асимптотично стійка . 2. Якщо речова частина хоча б одного кореня рівняння (3) позитивна, то точка спокою системи (1) нестійка. Тип точки та характер стійкості Стійкий вузол, стійкий фокус Сідло, Нестійкий вузол, Нестійкий фокус 3. Якщо рівняння (3) має чисто уявне коріння, то точка спокою системи (1) стійка, але не асимптотично. Центр

9 Фазові портрети 9 Стійкий вузол 1 2, 1< 0, 2 < 0 Неустойчивый узел 1 2, 1 > 0, 2 >

10 Фазові портрети 10 Стійкий фокус 1,2 = i,< 0, 0 Неустойчивый фокус 1,2 = i, >0, 0 Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

11 Фазові портрети 11 Сідло 1 2, 1< 0, 2 >0 Центр 1,2 = i, 0 Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

12 Фазові портрети 12 Дикритичний вузол має місце для систем виду: dx ax, dt dy ay, dt коли a 0. При цьому 1 = 2 = a. Нестійкий дикритичний вузол Якщо a< 0, то узел асимптотически устойчив, если a >0, то нестійкий. Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

13 Фазові портрети 13 Вироджений вузол, якщо 1 = 2 0 і в системі (1) b 2 + c 2 0. Якщо 1< 0, то устойчивый Если 1 >0, то нестійкий Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

14 Безліч точок спокою 14 Якщо det A = 0, то система (1) має безліч положень рівноваги. При цьому можливі три випадки: Коріння рівняння (3) 1 1 = 0, = 2 = = 2 = 0 Визначення точок спокою Система (2) рівнозначна одному рівнянню виду x + y = 0 (2) рівносильна рівнянню x + y = 0 Геометричне місце точок спокою Пряма на фазовій площині: x + y = 0 Вся фазова площина Пряма x + y = 0 У другому випадку будь-яка точка спокою стійка за Ляпуновим. У першому випадку тільки, якщо 2< 0.

15 Фазові портрети 15 Пряма стійка точка спокою 1 = 0, 2< 0 Прямая неустойчивых точек покоя 1 = 0, 2 >0 Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

16 Фазові портрети 16 Пряма нестійких точок спокою 1 = 2 = 0 Фазові прямі будуть паралельні до прямої точок спокою (x + y = 0), якщо перший інтеграл рівняння dy cx dy dx ax by має вигляд x + y = C, де C довільна постійна . Напрямок на фазовій кривій вказує напрямок руху фазової точки по кривій при зростанні t.

17 Правила визначення типу точки спокою 17 Можна визначити тип точки спокою та характер її стійкості, не знаходячи власних значень матриці системи (1), а знаючи тільки її слід tr A та визначник det A. Визначник матриці det A< 0 tra 0 det A 2 tra det A 2 tra det A След матрицы tr A < 0 tr A >0 tr A< 0 tr A >0 tr A< 0 tr A = 0 tr A >0 Тип точки спокою Cедло Стійкий вузол (УУ) Нестійкий вузол (НУ) Дикритичний або вироджений УУ Дикритичний або вироджений НУ Стійкий фокус (УФ) Центр Нестійкий фокус (НФ)

18 Центр Біфуркаційна діаграма 18 det A det tra A 2 2 УУ УФ НФ НУ tr A З е д л о

19 19 Алгоритм побудови фазового портрета ЛДС (1) 1. Визначити положення рівноваги, розв'язавши систему рівнянь: ax by 0, cx dy Знайти власні значення матриці системи, розв'язавши характеристичне рівняння: 2 λ (ad)λ ad bc Визначити тип точки спокою та зробити висновок про стійкість. 4. Знайти рівняння головних ізоклін горизонтальної та вертикальної, та побудувати їх на фазовій площині. 5. Якщо положення рівноваги є сідлом або вузлом, знайти фазові траєкторії, які лежать на прямих, що проходять через початок координат. 6. Намалювати фазові траєкторії. 7. Визначити напрямок руху фазовими траєкторіями, вказавши його стрілками на фазовому портреті.

20 Головні ізокліни 20 Вертикальна ізокліна (ВІ) сукупність точок фазової площини, в яких дотична, проведена до фазової траєкторії, паралельна вертикальній осі. Так як у цих точках фазових траєкторій x (t) = 0, то для ЛДС (1) рівняння ВІ має вигляд: ax + by = 0. . Так як у цих точках фазових траєкторій y(t) = 0, то для ЛДС (1) рівняння ГІ має вигляд: cx + dy = 0. Зауважимо, що точка спокою на фазовій площині це перетин головних ізоклін. Вертикальну ізокліну на фазовій площині позначатимемо вертикальними штрихами, а горизонтальну горизонтальними.

21 Фазові траєкторії 21 Якщо положення рівноваги є сідлом або вузлом, то існують фазові траєкторії, які лежать на прямих, що проходять через початок координат. Рівняння таких прямих можна шукати у вигляді * y = k x. Підставляючи y = k x до рівняння: dy cx dy, dx ax by для визначення k отримаємо: (4) c kd () 0. Даємо опис фазових траєкторій в залежності від кількості та кратності коренів рівняння (4). * Рівняння прямих, що містять фазові траєкторії, можна шукати і у вигляді x = k y. ak b ck d Тоді знаходження коефіцієнтів слід розв'язати рівняння k.

22 Фазові траєкторії 22 Корені рівняння (4) k 1 k 2 Тип точки спокою Сідло Вузол Опис фазових траєкторій Прямі y = k 1 x та y = k 2 x називають сепаратрисами. Інші фазові траєкторії це гіперболи, для яких знайдені прямі є асимптотами. Прямі y = k 1 x і y = k 2 x. Інші фазові траєкторії утворюють параболи, які стосуються початку координат однієї зі знайдених прямих. Фазові траєкторії стосуються тієї прямої, яка спрямована вздовж власного вектора, що відповідає меншому по абсолютної величини(корінь рівняння (3))

23 Фазові траєкторії 23 Коріння рівняння (4) k 1 k 2! k 1 Тип точки спокою Вироджений вузол Седло Вузол Опис фазових траєкторій Пряма y = k 1 x. Інші фазові траєкторії - це гілки парабол, які стосуються на початку координат цієї прямої прямої * y = k 1 x і x = 0 це сепаратриси. Інші фазові траєкторії гіперболи, для яких знайдені прямі є асимптотами Прямі* y = k 1 x і x = 0. Інші фазові траєкторії утворюють параболи, які стосуються початку координат однієї зі знайдених прямих. * Якщо рівняння прямих шукаються як x = k y, тоді це будуть прямі x = k 1 y і y = 0.

24 Фазові траєкторії 24 Коріння рівняння (4) kr Тип точки спокою Дикритичний вузол Опис фазових траєкторій Усі фазові траєкторії лежать на прямих y = k x, kr. Якщо положення рівноваги є центром, фазові траєкторії є еліпсами. Якщо положення рівноваги є фокусом, фазові траєкторії є спіралями. У разі, коли ЛДС має пряму точку спокою, можна знайти рівняння всіх фазових траєкторій, вирішивши рівняння: dy cx dy dx ax by Його перший інтеграл x + y = C і визначає сімейство фазових прямих.

25 Напрям руху 25 Якщо положення рівноваги є вузлом або фокусом, то напрям руху по фазових траєкторіях визначається однозначно його стійкістю (на початок координат) або нестійкістю (від початку координат). Щоправда, у разі фокусу потрібно встановити ще й напрямок закручування (розкручування) спіралі за годинниковою або проти годинникової стрілки. Це можна зробити, наприклад, так. Визначити знак похідної y(t) у точках осі x. dy Коли cx 0, якщо x 0, то ордината точки, що рухається по фазовій траєкторії при перетині «позитивного променя осі x» зростає. Отже, «закручування (розкручування)» траєкторій відбувається проти годинникової стрілки. Коли dt dy dt y0 y0 cx 0 якщо x 0, то «закручування (розкручування)» траєкторій відбувається за годинниковою стрілкою.

26 Напрямок руху 26 Якщо положення рівноваги є центром, то напрям руху за фазовими траєкторіями (за годинниковою стрілкою або проти) можна визначити так само, як встановлюється напрямок «закручування (розкручування)» траєкторії у разі фокусу. У разі «сідла» рух однією з його сепаратрис відбувається у напрямку початку координат, по інший від початку координат. За рештою фазових траєкторій рух відбувається відповідно до руху по сепаратрисах. Отже, якщо положення рівноваги сідло, достатньо встановити напрямок руху по якій-небудь траєкторії. І далі можна однозначно встановити напрямок руху по всіх інших траєкторіях.

27 Напрямок руху (сідло) 27 Щоб встановити напрямок руху по фазових траєкторіях у разі сідла, можна скористатися одним із наступних способів: 1 спосіб Визначити, яка з двох сепаратрис відповідає негативному власному значенню. Рух нею відбувається до точки спокою. 2 спосіб Визначити, як змінюється абсциса точки, що рухається по будь-якій з сепаратрис. Наприклад, для y = k 1 x маємо: dx(abk1) t ax bk1x(a bk1) x, x(t) x(0) e. dt yk x 1 Якщо x(t) при t+, рух по сепаратрисі y = k 1 x відбувається до точки спокою. Якщо x(t) при t+, рух відбувається від точки спокою.

28 Напрямок руху (сідло) 28 3 спосіб Якщо вісь x не є сепаратрисою, визначити як змінюється ордината точки, що рухається по фазовій траєкторії при перетині осі x. Коли dy dt y0 cx 0, якщо x 0, то ордината точки зростає і, отже, рух фазовими траєкторіями, що перетинають позитивну частину осі x, відбувається знизу вгору. Якщо ж ордината зменшується, то рух відбуватиметься згори донизу. Якщо визначати напрямок руху по фазовій траєкторії, що перетинає вісь y, то краще аналізувати зміну абсциси точки, що рухається.

29 Напрямок руху 29 4 спосіб* Побудувати в довільній точці (x 0,y 0) фазову площину (відмінну від положення рівноваги) вектор швидкості: dx dy v, (ax0 by0, cx0 dy0). dt dt (x, y) 0 0 Його напрямок та вкаже напрямок руху по фазовій траєкторії, що проходить через точку (x 0,y 0) : (x 0, y 0) v * Цей спосіб може бути використаний при визначенні напрямку руху по фазовому траєкторій для будь-якого типу точки спокою.

30 Напрямок руху 30 5 спосіб* Визначити області «знакопостійності» похідних: dx dt dy ax by cx dy. dt Межами цих областей будуть головні ізокліни. Знак похідної вкаже на те, як змінюється ордината і абсцису точки, що рухається по фазовій траєкторії в різних областях. y y x (t)<0, y (t)>0 x (t)<0, y (t)<0 x x x (t)>0, y(t)>0 x (t)>0, y(t)<0 * Этот способ может быть использован при определении направления движения по фазовым траекториям для любого типа точки покоя.

31 Приклад dx dt dy dt 2x 2 y, x 2y 1. Система має єдине нульове положення рівноваги, оскільки det A = Побудувавши відповідне характеристичне рівняння 26 = 0, знайдемо його коріння 1,26. Отже, положення рівноваги сідло. 3. Сепаратриси сідла шукаємо як y = kx. 4. Вертикальна ізокліну: x + y = 0. Горизонтальна ізокліну: x 2y = 0. Коріння речове та різного знака. 1 2k 2 6 k k k k k k 2 2k ,2, 1 2, 22, 2 0, 22.

32 Приклад 1 (сідло) 32 Намалюємо на фазовій площині сепаратриси y = k 1 x та y = k 2 x і головні ізокліни. y x Решту площини заповнюють траєкторії – гіперболи, для яких сепаратриси є асимптотами.

33 Приклад 1 (сідло) 33 y x Знайдемо напрямок руху траєкторіями. Для цього можна визначити знак похідної y(t) у точках осі x. При y = 0 маємо: dy dt y0 x 0, якщо x 0. Таким чином, ордината точки, що рухається по фазовій траєкторії при перетині «позитивного променя осі x» зменшується. Отже, рух фазовими траєкторіями, що перетинають позитивну частину осі x, відбувається зверху вниз.

34 Приклад 1 (сідло) 34 Тепер легко встановити напрямок руху по інших траєкторіях. y x

35 Приклад dx 4x2 y, dt dy x3y dt 1. Система має єдине нульове положення рівноваги, тому що det A = Побудувавши відповідне характеристичне рівняння = 0, знайдемо його коріння 1 = 2, 2 = 5. Отже, положення рівноваги нестійкий вузол. 3. Прямі: y = kx. 1 3k 1 k k k k k 4 2k , Вертикальна ізокліну: 2x + y = 0. Горизонтальна ізокліну: x + 3y = 0.

36 Приклад 2 (нестійкий вузол) 36 yx Оскільки 1 = 2 є меншим за абсолютною величиною, то, знайшовши відповідний йому власний вектор = (a 1,a 2) т: 4 2 a1 a1 2 a1 a2 0, 1 3 aa 2 2 = (1,1) т, встановимо, що інші фазові траєкторії, що утворюють параболи, стосуються початку координат прямої y = x. Нестійкість положення рівноваги однозначно визначає напрямок руху від точки спокою.

37 Приклад 2 (нестійкий вузол) 37 Так як 1 = 2 є меншим за абсолютною величиною, то, знайшовши відповідний йому власний вектор = (a 1,a 2) т: 4 2 a1 a1 2 a1 a2 0, 1 3 aa 2 2 = (1,1) т, встановимо, що інші фазові траєкторії, що утворюють параболи, стосуються початку координат прямої y = x. Нестійкість положення рівноваги однозначно визначає напрямок руху від точки спокою. y x

38 Приклад dx x 4 y, dt dy 4x2y dt 1. Система має єдине нульове положення рівноваги, оскільки det A = Побудувавши відповідне характеристичне рівняння = 0, знайдемо його дискримінант D. Оскільки D< 0, то корни уравнения комплексные, причем Re 1,2 = 3/2. Следовательно, положение равновесия устойчивый фокус. 3. Вертикальная изоклина: x 4y = 0. Горизонтальная изоклина: 2x y 0. Фазовые траектории являются спиралями, движение по которым происходит к началу координат. Направления «закручивания траекторий» можно определить следующим образом.

39 Приклад 3 (стійкий фокус) 39 Визначимо знак похідної y(t) у точках осі x. При y = 0 маємо: dy 4x 0, якщо x 0. dt y0 y Таким чином, ордината точки, що рухається по фазовій траєкторії при перетині «позитивного променя осі x» зростає. Значить, закручування траєкторій відбувається проти годинникової стрілки. x

40 Приклад dx x4 y, dt dy x y dt 1. Система має єдине нульове положення рівноваги, оскільки det A = Побудувавши відповідне характеристичне рівняння 23 = 0, знайдемо його коріння 1,2 = i3. Отже, становище рівноваги центр. 3. Вертикальна ізокліна: x 4y = 0. Горизонтальна ізокліна: x y 0. Фазові траєкторії системи еліпси. Напрямок руху по них можна встановити, наприклад, так.

41 Приклад 4 (центр) 41 Визначимо знак похідної y(t) у точках осі x. При y = 0 маємо: dy dt y0 x 0, якщо x 0. y Таким чином, ордината точки, що рухається по фазовій траєкторії при перетині «позитивного променя осі x» зростає. Отже, рух еліпсами відбувається проти годинникової стрілки. x

42 Приклад 5 (вироджений вузол) 42 dx xy, dt dy x3y dt 1. Система має єдине нульове положення рівноваги, тому що det A = Побудувавши відповідне характеристичне рівняння = 0, знайдемо його коріння 1 = 2 = 2. Отже, положення рівноваги стійке вироджений вузол. 3. Пряма: y = kx. 13k k 2 k k k k1,2 4. Вертикальна ізокліну: x + y = 0. Горизонтальна ізокліну: x 3y = 0.

43 Приклад 5 (вироджений вузол) 43 y x Намалюємо на фазовій площині ізокліни і пряму фазову траєкторію. Решта площини заповнюється траєкторіями, що лежать на гілках парабол, що стосуються прямої y = x.

44 Приклад 5 (вироджений вузол) 44 Стійкість положення рівноваги однозначно визначає напрямок руху до початку координат. y x

45 Приклад dx 4x 2 y, dt dy 2x y dt Оскільки визначник матриці системи det A = 0, система має нескінченно багато положень рівноваги. Усі вони лежать на прямій y 2 x. Побудувавши відповідне характеристичне рівняння 25 = 0, знайдемо його коріння 1 = 0, 2 = 5. Отже, всі положення рівноваги стійкі за Ляпуновим. Побудуємо рівняння інших фазових траєкторій: dy 2x y dy 1 1 =, y x C. dx 4x 2y dx Таким чином, фазові траєкторії лежать на прямих y x C, C const. 2

46 Приклад Напрямок руху однозначно визначається стійкістю точок прямої y 2 x. y x

47 Приклад dx 2 x y, dt dy 4x2y dt Оскільки визначник матриці системи det A = 0, система має нескінченно багато положень рівноваги. Усі вони лежать на прямій y 2 x. Оскільки і слід матриці системи tr A, то коріння характеристичного рівняння 1 = 2 = 0. Отже, усі положення рівноваги нестійкі. Побудуємо рівняння інших фазових траєкторій: dy 4x 2 y dy, 2, y 2 x C. dx 2x y dx Таким чином, фазові траєкторії лежать на прямих y 2 x C, C const, і паралельні прямий точок спокою. Встановимо напрямок руху по траєкторіях наступним чином.

48 Приклад Визначимо знак похідної y(t) у точках осі x. При y = 0 маємо: dy 0, якщо x 0, 4 x dt y0 0, якщо x 0. Таким чином, ордината точки, що рухається по фазовій траєкторії при перетині «позитивного променя осі x» зростає, а «негативного» зменшується. Значить рух фазовими траєкторіями правіше прямої точок спокою буде знизу вгору, а лівіше зверху вниз. y x

49 Вправи 49 Вправа 1. Для заданих систем визначте тип та характер стійкості положення рівноваги. Побудуйте фазові портрети. 1. dx 3, 3. dx 2 5, 5. dx x y x y 2 x y dt dt dt dy dy dy 6x 5 y; 2x 2 y; 4x 2 y; dt dt dt 2. dx, 4. dx 3, 6. dx x x y 2x 2 y; dt dt dt dy dy dy 2 x y; x y; x y. dt dt dt Вправа 2. При яких значеннях параметра a система dx dy 2 ax y, ay 2ax dt dt має положення рівноваги і воно є сідлом? вузлом? фокусом? Який у своїй система має фазовий портрет?

50 Неоднорідні ЛДС 50 Розглянемо лінійну неоднорідну систему (НЛДС) з постійними коефіцієнтами: dx ax by; 5) положення рівноваги. Якщо det A 0, система має єдине положення рівноваги P(x 0,y 0). Якщо det A 0, то система або має нескінченно багато положень рівноваги точки прямої, що визначається рівнянням ax + by + = 0 (або cx + dy + = 0), або взагалі не має положень рівноваги.

51 Перетворення НЛДС 51 Якщо система (5) має положення рівноваги, то виконавши заміну змінних: xx0, y y0, де, у випадку, коли система (5) має безліч положень рівноваги, x 0, y 0 координати будь-якої точки, що належить прямій точок спокою, отримаємо однорідну систему: dab, (6) dt dc d. dt Ввівши на фазовій площині x0y нову систему координат із центром у точці спокою P, побудуємо у ній фазовий портрет системи (6). У результаті площині x0y отримаємо фазовий портрет системи (5).

52 Приклад dx 2x 2y12, dt dy x 2y 3 dt Оскільки 2x 2y 12 0, x 3, x 2y 3 0 y 3, то ДС має єдине положення рівноваги P(3;3). Виконавши заміну змінних x = + 3, y = + 3, отримаємо систему: d 2 2 dt d 2 dt нульове положення якої нестійке і є сідлом (див. приклад 1).

53 Приклад Побудувавши фазовий портрет на площині P, сумісний з фазовою площиною x0y, знаючи, які координати має у ній точка P. y P x

54 Фазові портрети НЛДС 54 При побудові фазових портретів у разі, коли система (5) не має положень рівноваги, можна використовувати такі рекомендації: 1. Знайти перший інтеграл рівняння dx dy, ax by cx dy і таким чином визначити сімейство всіх фазових траєкторій. 2. Знайти основні изоклини: ax by 0 (ВІ), cx dy 0 (ГІ). 3. Знайти прямі фазові траєкторії, що містять, у вигляді у = kx +. При цьому знаходження коефіцієнтів k і, враховуючи, що c: a d: b, побудувати рівняння: dy (ax by) k. dx y kx ax by (a kb) x b y kx

55 Фазові портрети НЛДС 55 Оскільки вираз (a kb) x b залежить від x, якщо a + kb = 0, то отримаємо такі умови перебування k і: a kb 0, k. b Рівняння прямої можна шукати у вигляді x = ky +. Умови визначення k і будуються аналогічно. Якщо існує лише одна пряма, то вона є асимптотою для інших траєкторій. 2. Для визначення напрямку руху фазовими траєкторіями визначити області «знакопостійності» правих частин системи (5). 3. Для визначення характеру опуклості (увігнутості) фазових траєкторій побудувати похідну y(x) та встановити області її «знакопостійності». Різні прийоми побудови фазових портретів розглянемо з прикладів.

56 Приклад dx dt dy dt 0, 1. y Розв'язавши рівняння: dx dy 0 0, 1 отримаємо, що всі фазові траєкторії лежать на прямих x C, C R. Оскільки y (t) = 1 > 0, то ордината точки, що рухається за будь-якою фазовою траєкторією зростає. Отже, рух фазовими траєкторіями відбувається знизу вгору. x

57 Приклад dx dt dy dt 2, 2. y Розв'язавши рівняння: dy dx 2 1, 2 отримаємо, що всі фазові траєкторії лежать на прямих y x + C, C R. Оскільки y (t)< 0, то ордината движущейся точки по любой фазовой траектории убывает. Следовательно, движение по фазовым траекториям происходит сверху вниз. x

58 Приклад dx 1, dt dy x 1. dt Розв'язавши рівняння: dy x 1, dx 2 (x 1) y C, CR, 2 отримаємо, що фазовими траєкторіями системи є параболи: осі яких лежать на горизонтальній ізоклині x 1 0, а гілки спрямовані нагору. Так як x (t) 1 > 0, то абсцису точки, що рухається по будь-якій фазової траєкторії зростає. Отже, рух лівою гілкою параболи відбувається зверху вниз до перетину з прямою горизонтальною ізоклиною, а далі знизу вгору.

59 Приклад y Визначити напрямок руху фазовими траєкторіями можна було б і встановивши області «знакопостійності» правих частин системи. y 1 x x"(t) > 0, y"(t)< 0 x"(t) >0, y"(t) > 0 x 1

60 Приклад dx y, dt dy y 1. dt Вертикальна ізокліну y = 0; горизонтальна ізокліна y 1 = 0. З'ясуємо, чи існують прямі, які містять фазові траєкторії. Рівняння таких прямих шукатимемо у вигляді y = kx + b. Оскільки k dy y , dx y y kx b ykxb ykxb ykxb, то останній вираз не залежить від x, якщо k = 0. Тоді для знаходження b отримаємо b 1. Таким чином, на прямій y = 1 лежать фазові траєкторії. Ця пряма є асимптотою на фазовій площині.

61 Приклад Встановимо, який характер опуклості (увігнутості) мають фазові траєкторії щодо осі x. Для цього знайдемо похідну y(x): y(x) > 0 y 1 1 "() 1 1, dx dx y dx y y 2 d y d y d y x y і визначимо області «знакопостійності» отриманого виразу. У тих областях, де y (x) >< 0, выпуклость «вверх». y (x) < 0 y (x) >0 x

62 Приклад З'ясуємо напрямки руху фазовими траєкторіями, визначивши області «знакопостійності» правих частин системи dx y, dt dy y 1. dt Кордонами цих областей будуть вертикальна і горизонтальна ізокліни. Отриманої інформації достатньо побудови фазового портрета. y x (t) > 0, y (t) > 0 y (x) > 0 x (t) > 0, y (t)< 0, y (x) < 0 x (t) >0, y(t)< 0 y (x) >0 x

63 Приклад x (t) > 0, y (t) > 0 y (x) > 0 y y x (t) > 0, y (t)< 0, y (x) < 0 x x x (t) >0, y(t)< 0 y (x) > 0

64 Приклад dx 2, dt dy 2 x y. dt Горизонтальна ізокліна: 2x y = 0. З'ясуємо, чи існують прямі, які містять фазові траєкторії. Рівняння таких прямих шукатимемо у вигляді y = kx + b. Оскільки dy 2 xy (2 k) xbk, 2 2 dx y kx by kx b останній вираз не залежить від x, якщо k = 2. Тоді для знаходження b отримаємо b 2 b 4. 2 Таким чином, на прямій y = 2x4 лежать фазові траєкторії. Ця пряма є асимптотою на фазовій площині.

65 Приклад Встановимо, який характер опуклості (увігнутості) мають фазові траєкторії щодо осі x. Для цього знайдемо похідну y(x): 2 d y d x y y x x y y x dx "() dx Визначимо області «знакопостійності» отриманого виразу. У тих областях, де y(x) > 0, фазові траєкторії мають опуклість «вниз», а де y (x)< 0, выпуклость «вверх». y (x) >0 y x y (x)< 0

66 Приклад З'ясуємо напрямок руху фазовими траєкторіями, визначивши області «знакопостійності» правих частин системи: dx 2, dt dy 2 x y. dt Кордоном цих областей буде горизонтальна ізокліну. x(t)>0, y(t)<0 y x (t)>0, y(t)>0 x Отриманої інформації достатньо для побудови фазового портрета.

67 Приклад y (x) > 0 y x y y (x)< 0 x x (t)>0, y(t)<0 y x x (t)>0, y(t)>0

68 Приклад dx x y, dt dy 2(x y) 2. dt Вертикальна ізокліну: x y = 0; горизонтальна ізокліну: x y + 1= 0. З'ясуємо, чи існують прямі, які містять фазові траєкторії. Рівняння таких прямих шукатимемо у вигляді y = kx + b. Так як dy 2(xy) k 2 2, dx xyxy (1 k) xb ykxb ykxb ykxb то останній вираз не залежить від x, якщо k = 1. Тоді для знаходження b отримаємо b 2. Таким чином, на прямий y = x +2 лежать фазові траєкторії. Ця пряма є асимптотою на фазовій площині.

69 Приклад Визначимо, як змінюються абсциса і ордината точки, що рухається по фазовій траєкторії. Для цього збудуємо галузі «знакопостійності» правих частин системи. y x (t)<0, y (t)<0 x (t)<0, y (t)>0 x x (t)>0, y (t)>0 Ця інформація буде потрібна для визначення напрямку руху по траєкторіях.

70 Приклад Встановимо, який характер опуклості (увігнутості) мають фазові траєкторії щодо осі x. Для цього знайдемо похідну y(x): 2(xy) () 2 2("() 1) xy 2(2) dx dx xy (xy) (xy) (xy) 2 dydxyyxxy Визначимо області «знакопостійності» отриманого виразу. У тих областях, де y(x) > 0, фазові траєкторії мають опуклість «вниз», а де y(x)< 0, выпуклость «вверх». y (x)>0 y y (x)< 0 x Полученной информации достаточно для построения фазового портрета. y (x)> 0

71 Приклад 14 (ФП) 71 y y x y x x

72 Вправи 72 Побудуйте фазові портрети для наступних систем: dx 3x 3, dt dy 2x y1; dt dx x; dt dy 2x 4; dt dx x y 2; dt dy 2x 2y1; dt dx 1, dt dy 2 x y; dt dx dt dy dt dx dt dy 2, 4; y 2, 2.

73 Література 73 Понтрягін Л.С. Прості диференціальні рівняння. М., Філіппов А.Ф. Збірник задач з диференціальних рівнянь. М., Пантелєєв А.В., Якімова А.С., Босов А.В. Звичайні диференціальні рівняння у прикладах та завданнях. М: Вища. шк., 2001.

4.03.07 Заняття 4. Існування та стійкість положень рівноваги лінійних динамічних (ЛДС) систем на площині. Побудувати параметричний портрет та відповідні фазові портрети ЛДС (x, yr, ar):

Семінар 4 Система двох звичайних диференціальних рівнянь (ОДП). Фазова площина. Фазовий портрет. Кінетичні криві. Особливі точки. Стійкість стаціонарного стану. Лінеаризація системи в

Математичні методи в екології: Збірник завдань та вправ / Упоряд. Є.Є. Семенова, Є.В. Кудрявцева. Петрозаводськ: Изд-во ПетрГУ, 005..04.09 Заняття 7 Модель «хижак-жертва» Лотки-Вольтерри 86 (побудова

РОСІЙСЬКИЙ ТЕХНОЛОГІЧНИЙ УНІВЕРСИТЕТ МИРЕА ДОДАТКОВІ ГЛАВИ ВИЩОЇ МАТЕМАТИКИ РОЗДІЛ 5. ТОЧКИ СПОКІВ Робота присвячена моделюванню динамічних систем з використанням елементів вищої математики

Система лінійних диференціальних рівнянь із постійними коефіцієнтами. Кільцов С.М. www.linis.ru Метод варіації довільних постійних. Розглянемо лінійне неоднорідне диференціальне рівняння:

Стор. Лекція 3 СТІЙКІСТЬ РІШЕННЯ СИСТЕМ ДК Якщо деяке явище описується системою ДК dx dt i = f (t, x, x...x), i =..n з початковими in умовами xi (t 0) = x i0, i =. n, які зазвичай є

4.04.7 Заняття 7. Стійкість положень рівноваги автономних систем (метод лінеаризації Ляпунова, теорема Ляпунова) x "(f(x, y), f, g C(). y"(g(x, y), D) Пошук положень рівноваги P (x*, : f

СЕМІНАР 5 І 6 Система двох автономних звичайних лінійних диференціальних рівнянь. Фазова площина. Ізокліни. Побудова фазових портретів. Кінетичні криві. Знайомство із програмою TRAX. Фазовий

Лекція 6. Класифікація точок спокою лінійної системи двох рівнянь із постійними дійсними коефіцієнтами. Розглянемо систему двох лінійних диференціальних рівнянь із постійними дійсними

СЕМІНАР 4 Система двох автономних звичайних лінійних диференціальних рівнянь (ОДП). Рішення системи двох лінійних автономних ОДУ. Типи спеціальних точок. РІШЕННЯ СИСТЕМИ ЛІНІЙНИХ ДИФЕРЕНЦІЙНИХ РІВНЯНЬ

Міністерство освіти та науки Російської ФедераціїФедеральне державне бюджетне освітня установа вищої освіти«Уфімський державний нафтовий технічний університет» Кафедра

Лекція 1 Елементи якісного аналізу динамічних систем з безперервним часом на прямій Розглянемо автономне диференціальне рівняння du = f(u), (1) dt яке може бути використане

СЕМІНАР 7 Дослідження стійкості стаціонарних станів нелінійних систем другого порядку. Класична система Ст. Вольтерра. Аналітичне дослідження (визначення стаціонарних станів та їх стійкості)

Особливі точки в системах другого та третього порядків. Критерії стійкості стаціонарних станів лінійних та нелінійних систем. План відповіді Визначення особливої ​​точки типу центр. Визначення особливої ​​точки

ПРАКТИЧНІ ЗАНЯТТЯ З ДИФЕРЕНЦІЙНИХ РІВНЯНЬ Методична розробкаУкладач: проф АН Саламатін На основі: АФ Філіппов Збірник завдань з диференціальних рівнянь Москва-Іжевськ НДЦ "Регулярна

1 ЛЕКЦІЯ 2 Системи нелінійних диференціальних рівнянь. Простір станів чи фазовий простір. Особливі точки та його класифікація. Умови сталості. Вузол, фокус, сідло, центр, граничний цикл.

7 ПОЛОЖЕННЯ РІВНОВАГИ ЛІНІЙНИХ АВТОНОМНИХ СИСТЕМ ДРУГОГО ПОРЯДКУ Автономною системою для функцій (t) (t) називається система диференціальних рівнянь d d P() Q() (7) dt dt де праві частини не залежать

Міністерство освіти і науки Російської Федерації Ярославський державний університетім. П. Г. Демидова Кафедра алгебри та математичної логіки С. І. Яблокова Криві другого порядку Частина Практикум

Розділ IV. Перші інтеграли систем ОДУ 1. Перші інтеграли автономних систем звичайних диференціальних рівнянь У цьому параграфі будемо розглядати автономні системивиду f x = f 1 x, f n x C 1

Лекція 9 Лінеаризація диффе6ренціальних рівнянь Лінійні диференціальні рівняння вищих порядків Однорідні рівняння властивості їх розв'язків Властивості розв'язків неоднорідних рівнянь Визначення 9 Лінійним

Побудова інтегральних кривих та фазового портрета автономного рівняння Маючи графік гладкої функції f(u), можна схематично побудувати інтегральні криві рівняння du dt = f(u). (1) Побудова спирається

7.0.07 Заняття. Динамічні системи з безперервним часом на прямий. Завдання 4. Побудувати біфуркаційну діаграму та типові фазові портрети для динамічної системи: d dt Розв'язування рівняння f (, 5 5,

Теорія стійкості Ляпунова. Багато завдань механіки і техніки буває важливо знати не конкретні значення рішення при даному конкретному значенні аргументу, а характер поведінки рішення при зміні

Стор. 1 із 17 26.10.2012 11:39 Атестаційне тестування у сфері професійної освітиСпеціальність: 010300.62 Математика. Комп'ютерні науки Дисципліна: Диференціальні рівняння Час виконання

Семінар 5 Моделі, що описуються системами двох автономних диференціальних рівнянь. Вивчення нелінійних систем другого порядку. Модель Лотки. Модель Вольтерри. У загальному вигляді моделі, що описуються системами

Семінар Диференціальне рівняння першого ладу. Фазовий простір. Фазові змінні. Стаціонарний стан. Стійкість стаціонарного стану за Ляпуновим. Лінеаризація системи на околиці

Математичний аналіз Розділ: диференціальні рівняння Тема: Поняття стійкості рішення ДК та рішення системи ДК Лектор Пахомова Є.Г. 2012 5. Поняття стійкості рішення 1. Попередні зауваження

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

РОСІЙСЬКИЙ ТЕХНОЛОГІЧНИЙ УНІВЕРСИТЕТ МИРЕА ДОДАТКОВІ ГЛАВИ ВИЩОЇ МАТЕМАТИКИ РОЗДІЛ 3. СИСТЕМИ ДИФЕРЕНЦІЙНИХ РІВНЕНЬ Робота присвячена моделюванню динамічних систем з використанням елементів

КВАДРАТНІ РІВНЯННЯ Зміст КВАДРАТНІ РІВНЯННЯ... 4. та дослідження квадратних рівнянь... 4.. Квадратне рівняння з числовими коефіцієнтами... 4.. Вирішити та дослідити квадратні рівняння щодо

7..5,..5 Заняття,. Дискретні динамічні системи на прямій Завдання Провести дослідження динаміки густини популяції (t), що описується рівнянням: t t, const. t Чи існують серед рішень рівняння

Дослідження функції та побудова її графіка Пункти Дослідження: 1) Область визначення, безперервність, парність/непарність, періодичність функції. 2) Асимптоти графіка функції. 3) Нулі функції, інтервали

ЛЕКЦІЯ 16 ЗАДАЧА ПРО СТІЙКІСТЬ ПОЛОЖЕННЯ РІВНОВАГИ У КОНСЕРВАТИВНІЙ СИСТЕМІ 1. Теорема Лагранжа про стійкість положення рівноваги консервативної системи Нехай є n ступенів свободи. q 1, q 2,

Криві другого порядку Окружність Еліпс Гіпербола Парабола Нехай на площині задана прямокутна декартова система координат. Кривий другого порядку називається безліч точок, координати яких задовольняють

Лекція 1 Диференціальні рівняння першого порядку 1 Поняття диференціального рівняння та його рішення Звичайним диференціальним рівнянням 1-го порядку називається вираз виду F(x, y, y) 0, де

Тема 41 «Завдання з параметром» Основні формулювання завдань з параметром: 1) Знайти всі значення параметра, при кожному з яких виконується певна умова.) Розв'язати рівняння чи нерівність з

Лекція 3. Фазові потоки на площині 1. Стаціонарні точки, лінеаризація та стійкість. 2. Граничні цикли. 3. Біфуркація фазових потоків на площині. 1. Стаціонарні точки, лінеаризація та стійкість.

Лекція 3 Стійкість рівноваги і руху системи При розгляді рухів рівняння збуреного руху, що встановилися, запишемо у вигляді d dt A Y де вектор-стовпець квадратна матриця постійних коефіцієнтів

5. Стійкість атракторів 1 5. Стійкість атракторів У минулому розділі ми навчилися знаходити нерухомі точки динамічних систем. Також ми з'ясували, що існує кілька різних типів нерухомих

4 лютого 9 г Практичне заняття Найпростіші завдання управління динамікою популяцій Завдання Нехай вільний розвиток популяції описується моделлю Мальтуса N N де N чисельність чи обсяг біомаси популяції

1) Привести рівняння кривої другого порядку x 4x y 0 до канонічного вигляду та знайти точки перетину її з прямою x y 0. Виконайте графічну ілюстрацію отриманого рішення. x 4x y 0 x x 1 y 0 x 1 y 4

РОЗДІЛ 4 Системи звичайних диференціальних рівнянь ЗАГАЛЬНІ ПОНЯТТЯ ТА ВИЗНАЧЕННЯ Основні визначення Для опису деяких процесів і явищ нерідко потрібно кілька функцій.

Семінар 9 Лінійний аналіз стійкості гомогенного стаціонарного стану системи двох рівнянь реакція дифузія Нестійкість Тьюринга Активатор та інгібітор Умови виникнення дисипативних структур

ЛЕКЦІЯ 17 КРИТЕРІЙ РАУСА-ГУРВІЦЯ. МАЛІ КОЛИВАННЯ 1. Стійкість лінійної системи Розглянемо систему двох рівнянь. Рівняння обуреного руху мають вигляд: dx 1 dt = x + ax 3 1 dx dt = x 1 + ax 3,

МІНІСТЕРСТВО ОСВІТИ І НАУКИ РФ НОВОСИБІРСЬКИЙ ДЕРЖАВНИЙ УНІВЕРСИТЕТ Фізичний факультет Кафедра вищої математики фізичного факультету Методи вирішення звичайних диференціальних рівнянь.

1. Що таке прості диференціальні рівняння та системи. Концепція рішення. Автономні та неавтономні рівняння. Рівняння та системи порядку вищі за перший та їх зведення до систем першого порядку.

Лекція 1 Дослідження руху на консервативної системі з одним ступенем свободи 1. Основні поняття. Консервативною системою з одним ступенем свободи ми називатимемо систему, що описується диференціальним

РОЗДІЛ. СТІЙКІСТЬ ЛІНІЙНИХ СИСТЕМ 8 ступінь зі знаком +, отриманого випливає, що () π зростає від до π. Отже, доданки ϕ i() і k () +, тобто вектор (i) ϕ монотонно ϕ монотонно зростають при

ФАЗОВА ПЛОЩІСТЬ ДЛЯ НЕЛІНІЙНОГО АВТОНОМНОГО РІВНЯННЯ -ГО ПОРЯДКУ.. Постановка задачі. Розглянемо автономне рівняння виду = f. () Як відомо, це рівняння еквівалентно наступній нормальній системі

ДИФЕРЕНЦІАЛЬНІ РІВНЯННЯ 1. Основні поняття Диференціальним рівнянням щодо деякої функції називається рівняння, що пов'язує цю функцію з її незалежними перемпними та з її похідними.

Математичні методи в екології: Збірник завдань та вправ / Упоряд. Є.Є. Семенова, Є.В. Кудрявцева. Петрозаводськ: Изд-во ПетрГУ, 2005. 2 семестр Заняття. Модель "Хижак-жертва" Лотки-Вольтерри Тема 5.2.

Геометричний змістпохідна, дотична 1. На малюнку зображено графік функції y=f(x) та дотична до нього в точці з абсцисою x 0. Знайдіть значення похідної функції f(x) у точці x 0. Значення

Графік функції y=f(x) називається опуклим на інтервалі (a; b), якщо він розташований нижче за будь-яку свою дотичну на цьому інтервалі Графік

Розділ 6 Основи теорії стійкості Лекція Постановка задачі Основні поняття Раніше було показано, що розв'язання задачі Коші для нормальної системи ОДУ = f, () безперервно залежить від початкових умов при

19.11.15 Заняття 16. Базова модель «Брюсселятор» До початку 70-х років. більшість хіміків вважало, що хімічні реакції не можуть йти в коливальному режимі. Експериментальні дослідження радянських вчених

Розділ 8 Функції та графіки Змінні та залежності між ними. Дві величини і називаються прямо пропорційними, якщо їх відношення постійно, тобто якщо =, де постійне число, що не змінюється зі зміною

Система підготовки учнів до ЄДІ з математики профільного рівня. (Завдання з параметром) Теоретичний матеріал Визначення. Параметром називається незалежна змінна, значення якої завдання вважається

Лекція Дослідження функції та побудова її графіка Анотація: Функція досліджується на монотонність, екстремум, опуклість-увігнутість, існування асимптот Наводиться приклад дослідження функції, будується

29. Асимптотична стійкість рішень систем звичайних диференціальних рівнянь, сфера тяжіння та методи її оцінки. Теорема В.І. Зубова про межі сфери тяжіння. В.Д.Ногін 1 о. Визначення

Лекція 13 Тема: Криві другого ладу Криві другого ладу на площині: еліпс, гіпербола, парабола. Виведення рівнянь кривих другого порядку виходячи з них геометричних властивостей. Дослідження форми еліпса,

ЗАТВЕРДЖЕНО Проректор з навчальної роботи та довузівської підготовкиО. О. Воронов 09 січня 2018 р. ПРОГРАМА з дисципліни: Динамічні системи за напрямом підготовки: 03.03.01 «Прикладна математика

Вступ

Оскільки концепція нелінійної динамічної системи досить багата, щоб охопити надзвичайно широке коло процесів, у яких майбутня поведінка системи визначається минулим, методи аналізу, розроблені в цій галузі, корисні у величезній різноманітності контекстів

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

По-друге, є випадки, коли нелінійна динамічна теорія може використовуватися для твердження, що деяка спрощена модель має демонструвати важливі особливостіданої системи, з чого випливає, що описову модель можна побудувати та вивчити у широкому діапазоні параметрів. Часто це призводить до моделей, які поводяться якісно по-різному при різних параметрах і демонструють, що одна область показує поведінку, дуже схожу на поведінку, що спостерігається в реальній системі. У багатьох випадках поведінка моделі досить чутлива до змін параметрів, тому якщо параметри моделі можуть бути виміряні в реальній системі, модель демонструє реалістичну поведінку при цих значеннях, і можна бути впевненим, що модель охопила істотні особливості система.

По-третє, бувають випадки, коли модельні рівняння будуються з урахуванням докладних описів відомої фізики. Потім чисельні експерименти можуть дати інформацію про змінні, недоступні фізичним експериментам.

Спираючись на другий шлях, ця робота є розширенням моєї попередньої роботи «Нелінійна динамічна модель взаємозалежних галузей виробництва», а також іншої роботи (Dmitriev, 2015)

Усі необхідні визначення та інші теоретичні відомості, необхідні у роботі, з'являтимуться у першому розділі, у міру їхньої необхідності. Тут же буде наведено два визначення, що необхідні для розкриття самої теми дослідження.

Спочатку дамо визначення системної динаміки. Відповідно до одного з визначень, системна динаміка - підхід імітаційного моделювання, який завдяки своїм методам та інструментам допомагає оцінити структуру складних систем та їх динаміку (Штерман). Варто додати, що системна динаміка також і метод моделювання, який використовують з метою відтворення вірних (з точки зору точності) комп'ютерних моделей для складних систем заради їх майбутнього використання для того, щоб створити ефективнішу компанію/організацію, а також покращити методи взаємодії з даною системою. Переважно, необхідність у системній динаміці виникає при зіткненні з довгостроковими, стратегічними моделями, а також варто відзначити, що вона досить абстрактна.

Говорячи про нелінійну диференціальну динаміку, ми розглядатимемо нелінійну систему, яка за визначенням, є системою, в якій зміна результату не пропорційна зміні вхідних параметрів, і в якій функція описує залежність зміни в часі та положенні точки у просторі (Boeing, 2016).

Виходячи з вищеназваних визначень, стає зрозуміло, що ця робота розглядатиме різні нелінійні диференціальні системи, що описують взаємодію компаній, а також побудовані на їх основі імітаційні моделі. Грунтуючись на цьому і буде визначено мету роботи.

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

Для досягнення поставленої мети було виділено такі завдання:

Визначення сталості системи.

Побудова фазових портретів.

Знаходження інтегральних траєкторій систем.

Побудова імітаційних моделей.

Кожним із цих завдань буде присвячений один із розділів кожного з розділів роботи.

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

Що стосується причини вибору якісного, а не кількісного аналізу, то варто зазначити, що в переважній кількості випадків результати та висновки з якісного аналізу динамічних систем виявляються значнішими за результати їх кількісного аналізу. У такій ситуації доречно зазначити висловлювання В.П. Мілованова, у якому він стверджує, що зазвичай вважають, що результати, очікувані при застосуванні математичних методів для аналізу реальних об'єктів, повинні зводитися до чисельного результату. У цьому сенсі якісні методи покладається дещо інше завдання. У ньому акцентується увагу на досягненні результату, що описує якості системи, пошуку характерних особливостей всіх явищ в цілому, на прогнозування. Зрозуміло, важливо розуміти, як зміниться попит при зміні цін на певний вид товарів, але не варто забувати, що набагато важливіше розуміння, чи в таких умовах настане дефіцит чи надлишок цих товарів (Дмитрієв, 2016).

Об'єктом цього дослідження є нелінійна диференціальна та системна динаміка.

У такому разі предмет дослідження - опис процесу взаємодії між компаніями через нелінійну диференціальну та системну динаміку.

Говорячи про практичне застосування дослідження, варто відразу розбити його на дві частини. А саме на теоретичну, тобто якісний аналіз систем, та практичну, в якій буде розглянуто побудову імітаційних моделей.

Теоретична частина цього дослідження дає базові поняття та явища. У ній розглянуті прості диференціальні системи, як і в роботах багатьох інших авторів (Teschl, 2012; Nolte, 2015), але дозволяють описати взаємодію між компаніями. Грунтуючись на цьому надалі можна буде проводити більш поглиблені дослідження, або ж починати своє знайомство з тим, що являє собою якісний аналіз систем.

Практичну частину роботи можна використовуватиме створення системи підтримки прийняття рішень. Система підтримки прийняття рішень - автоматизована інформаційна система, орієнтована підтримку бізнесу чи прийняття рішень у створенні, що дозволяє вибрати між безліччю різних альтернатив (Keen, 1980). Нехай в даний момент моделі і не мають високої точності, але змінивши їх під конкретну компанію, можна досягти більш вірних результатів. Таким чином, при зміні в них різних параметріві умов, здатних виникнути над ринком, можна отримати прогноз на майбутнє і прийняти заздалегідь вигідніше рішення.

1. Взаємодія підприємств за умов мутуалізму

У роботі будуть представлені двовимірні системи, які досить прості порівняно з системами вищого порядку, але при цьому дозволяють продемонструвати необхідні нам взаємини між організаціями.

Почати роботу варто з вибору виду взаємодії, яка надалі і буде описуватися, оскільки для кожного з видів системи, що їх описують, нехай небагато, але різні. На малюнку 1.1, зазначена модифікована під економічну взаємодію класифікація Юджима Одума для взаємодії популяцій (Одум, 1968), ґрунтуючись на якій надалі ми розглядатимемо взаємодію компаній.

Малюнок 1.1. Типи взаємодії між підприємствами

Ґрунтуючись на малюнку 1.1, виділимо 4 типи взаємодії і наведемо для кожного з них, що описує їх систему рівнянь, заснованої на моделі Мальтуса (Malthus, 1798). Відповідно до неї, швидкість зростання має пропорційну залежність від поточної чисельності виду, іншими словами, її можна описати наступним диференціальним рівнянням:

де a - певний параметр, що залежить від природного приросту популяції. Також варто додати, що в системах, що розглядаються далі, всі параметри, а також змінні приймають невід'ємні значення.

Виробництво сировини - виробництво продукції, що аналогічно моделі хижака-жертва. Модель хижак-жертва, також відома як модель Лотки-Вольтерри, - пара нелінійних диференціальних рівнянь першого порядку, що описують динаміку біологічної системи з двома видами, один із яких хижаки, а інший жертви (Llibre, 2007). Зміна чисельності цих видів описується наступною системою рівнянь:

(1.2)

де – характеризує зростання продукції першого підприємства без впливу другого (у разі моделі хижак-жертва, зростання популяції жертв без хижаків),

Характеризує зростання продукції другого підприємства без впливу першого (зростання популяції хижаків без жертв),

Характеризує зростання продукції першого підприємства з урахуванням впливом нього другого (зростання чисельності жертв при взаємодії з хижаками),

Характеризує зростання продукції другого підприємства, враховуючи вплив нього першого (зростання чисельності хижаків за її взаємодії з жертвами).

На одного, на хижака, як видно із системи, а також класифікації Одума, їхня взаємодія накладає сприятливий вплив. На іншого несприятливе. Якщо розглядати в економічних реаліях, то як видно на малюнку найпростішим аналогом є фірма-виробник та її постачальник ресурсів, які відповідають хижакові та жертві відповідно. Таким чином, без сировини випуск продукції експоненційно знижується.

Конкуренція - це суперництво між двома і більше (у разі ми розглядаємо двовимірні системи, тому беремо саме двовидову конкуренцію) видами, економічними групами за території, обмежені ресурси чи інші цінності (Elton, 1968). Зміни в кількості видів, або кількості продукції в нашому випадку, описується системою нижче:

(1.3)

В даному випадку види або компанії, що випускають один продукт, несприятливо впливають один на одного. Тобто, за відсутності конкурента, зростання продукції експоненційно зросте.

Тепер перейдемо до симбіотичного взаємодії, у якому обидва підприємства мають одне одного позитивний вплив. Спочатку розглянемо мутуалізм. Мутуалізм - тип відносин між різними видами, при якому кожен з них отримує вигоду від дії іншого, причому варто зазначити, що присутність партнера є обов'язковою умовою існування (Thompson, 2005). Такий тип відносин описується системою:

(1.4)

Оскільки взаємодія між компаніями необхідна їх існування, то відсутність товару однієї компанії, випуск товарів інший експоненційно знижується. Таке можливо, коли компанії просто не мають інших альтернатив для закупівель.

Розглянемо ще один вид симбіотичної взаємодії, протокооперацію. Протокооперація схожа на мутуалізм з єдиним винятком, немає потреби в обов'язковому існуванні партнера, оскільки, наприклад, існують інші альтернативи. Оскільки вони схожі, то їх системи виглядають практично аналогічно один одному:

(1.5)

Таким чином, відсутність продукту однієї компанії не заважає зростанню продукту іншої.

Звісно, ​​крім перелічених у пунктах 3 і 4, можна назвати й інші види симбіотичних відносин: комменсалізм і аменсалізм (Hanski, 1999). Але вони не згадуватимуться надалі, оскільки в коменсалізмі одному з партнерів байдуже його взаємодія з іншим, а ми таки розглядаємо випадки, коли вплив є. А аменсалізм не розглядається, оскільки з економічної точки зору таких відносин, коли одному їхня взаємодія шкодить, а іншому байдуже, просто не може бути.

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

Ця глава присвячена взаємодії підприємств за умов мутуалізму. У ній будуть розглянуті дві системи, що є подальшим розвитком систем, заснованих на моделі Мальтуса, а саме системами з обмеженнями на збільшення продукції, що накладаються.

Динаміку пари, пов'язаної мутуалістичними відносинами, як було зазначено вище, у першому наближенні можна описати системою:

(1.6)

Можна зауважити, що при великій початковій кількості продукції система необмежено зростає, а при малому виробництво продукції падає. У цьому полягає некоректність білінійного опису ефекту, що виникає при мутуалізмі. Щоб спробувати виправити картину, введемо фактор, що нагадує насичення хижака, тобто фактор, який дозволить зменшити швидкість зростання продукції, за його надлишку. У цьому випадку приходимо до наступної системи:

(1.7)

де - зростання виробництва продукту першої компанії при її взаємодії з другою з урахуванням насичення,

Зростання виробництва продукту другої компанії при її взаємодії з першою з урахуванням насичення,

Коефіцієнти насичення.

Таким чином ми і отримали дві системи: мальтузіанська модель зростання з насиченням і без нього.

1.1 Стійкість систем у першому наближенні

Стійкість систем у першому наближенні розглядається у багатьох, як іноземних (Hairer, 1993; Bhatia, 2002; Khalil, 2001; Strogatz, 2001 та інші), так і російськомовних роботах (Ахромєєва, 1992; Беллман, 1954;, Демид; 1959 та інші), та її визначення є базовим кроком для аналізу процесів, що відбуваються в системі. Для цього виконаємо такі кроки:

Знайдемо рівноважні точки.

Знайдемо матрицю Якобі системи.

Знайдемо власні значення матриці Якобі.

Класифікуємо рівноважні точки за теоремою Ляпунова.

Розглянувши кроки, варто докладніше зупинитися на їхнє роз'яснення, тому дам визначення та опишу методи, якими ми користуватимемося в кожному з цих кроків.

Перший крок – пошук рівноважних точок. Для їхнього знаходження прирівняємо кожну функцію до нуля. Тобто вирішимо систему:

де a і b маються на увазі всі параметри рівняння.

Наступний крок – пошук матриці Якобі. У нашому випадку це буде матриця 2 на 2 з першими похідними в деякій точці , як представлено нижче:

Після виконання перших двох кроків переходимо до знаходження коріння наступного характеристичного рівняння:

Де точка відповідає рівноважним точкам, знайденим у першому кроці.

Знайшовши і перейдемо до четвертого кроку і скористаємося наступними теоремами Ляпунова (Parks, 1992):

Теорема 1: Якщо всі корені характеристичного рівняння мають негативну дійсну частину, то рівноважна точка відповідна початковій та лінеаризованій системам – асимптотично стійка.

Теорема 2: Якщо хоча б один з коренів характеристичного рівняння має позитивну дійсну частину, то рівноважна точка, що відповідає початковій та лінеаризованій системам, - асимптотично нестійка.

Також, дивлячись на і можна і точніше визначити тип стійкості, ґрунтуючись на наведеному на рисунках 1.2 поділу (Lamar University).

Малюнок 1.2. Типи стійкості рівноважних точок

Розглянувши потрібні теоретичні відомості, перейдемо до аналізу систем.

Розглянемо систему без насичення:

Вона дуже проста та не підходить для практичного застосування, оскільки не має жодних обмежень. Але як перший приклад аналізу системи підходить до розгляду.

Спочатку знайдемо рівноважні точки, прирівнявши праві частини рівнянь до нуля. Таким чином, виявляємо дві рівноважні точки, назвемо їх A і B: .

Об'єднаємо крок із пошуком матриці Якобі, коренів характеристичного рівняння та визначенням типу стійкості. Оскільки вони елементарні, то одразу отримаємо відповідь:

1. У точці , , стійкий вузол.

У точці: , , сідло.

Як я вже писав, дана система занадто тривіальна, тому не потрібно ніяких пояснень.

Тепер проведемо аналіз системи з насичення:

(1.9)

Поява обмеження на взаємонасичення продукцією підприємствами наближає нас до реальної картини того, що відбувається, а також трохи ускладнює систему.

Як і раніше, прирівнюємо праві частини системи до нуля та вирішуємо отриману систему. Крапка залишилася без змін, а ось інша точка в даному випадку містить більше параметрів, ніж раніше: .

У цьому випадку матриця Якобі набуває такого вигляду:

Віднімемо від неї одиничну матрицю, помножену на і прирівняємо визначник отриманої матриці в точці А і B до нуля.

У точці аналогічну ранню картину:

Стійкий вузол.

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

2 Фазові портрети систем

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

У таких випадках, варто згадати про фазову площину і фазові портрети. У прикладній математиці, зокрема в контексті нелінійного системного аналізу, фазова площина є візуальним відображенням певних характеристик деяких видів диференціальних рівнянь (Nolte, 2015). Координатна площиназ осями значень будь-якої пари змінних, що характеризують стан системи – двовимірний випадок загального n-вимірного фазового простору.

Завдяки фазовій площині можна графічно визначити існування граничних циклів у рішеннях диференціального рівняння.

Рішення диференціального рівняння є сімейством функцій. Графічно це можна побудувати у фазовій площині як двовимірне векторне поле. На площині малюються вектори, що представляють похідні у характерних точках за яким-небудь параметром, у разі по часу, тобто (). При достатній кількості цих стрілок в одній області можна візуалізувати поведінку системи та легко ідентифікувати граничні цикли (Boeing, 2016).

Векторне поле є фазовим портретом, конкретний шлях вздовж лінії потоку (тобто шлях завжди дотичний до векторів) є фазовим шляхом. Потоки у векторному полі вказують на зміну системи у часі, що описується диференціальним рівнянням (Jordan, 2007).

Варто зазначити, що фазовий портрет можна побудувати, навіть без розв'язання диференціального рівняння, і в той же час хороша візуалізація може надати багато корисної інформації. До того ж нині існує безліч програм, здатних допомогти з побудовою фазових діаграм.

Отже, фазові площини корисні візуалізації поведінки фізичних систем. Зокрема, коливальних систем, таких як згадана вище модель хижак-жертва. У цих моделях фазові траєкторії можуть «закручуватися» у напрямі нуля, «виходити зі спіралі» в нескінченність або досягати стійкої нейтральної ситуації, званої центрами. Це корисно при визначенні того, чи стабільна динаміка чи ні (Jordan, 2007).

Представлені в цьому розділі фазові портрети будуть побудовані за допомогою інструментів WolframAlpha, або наведені з інших джерел. Мальтузіанська модель зростання без насичення.

Побудуємо фазовий портрет першої системи з трьома набором параметрів, щоб порівняти їхню поведінку. Набір А ((1,1), (1,1)), який надалі називатиметься одиничним набором, набір B ((10,0.1), (2,2)), при виборі якого в системі спостерігається різкий спад виробництва продукції , і набір C ((1,10), (1,10)), при якому навпаки виникає різке і необмежене зростання. Варто відзначити, що значення по осях у всіх випадках будуть перебувати в одних і тих же інтервалах від -10 до 10, для зручності порівняння фазових діаграм між собою. Звичайно, це не стосується якісного портрета системи, у якого осі безрозмірні.

Малюнок 1.3 Фазовий портрет із параметрами А

мутуалізм диференціальний граничний рівняння

На малюнку 1.3, наведених вище, продемонстровані фазові портрети системи при трьох зазначених наборах параметрів, а також фазовий портрет, що описує якісну поведінку системи. Не варто забувати, що найважливішою з практичної точки зору є перша чверть, оскільки кількість продукції, яка може бути лише невід'ємною, є наші осі.

На кожному з малюнків явно видно стійкість у рівноважній точці (0,0). І першому малюнку також помітно «сідло» у точці (1,1), інакше кажучи, якщо підставити значення набору параметрів у систему, то рівноважної точці У. При зміні меж побудови моделі, сідлова точка виявляється і інших фазових портретах.

Мальтузіанська модель зростання з насичення.

Побудуємо фазові діаграми для другої системи, в якій є насичення, з трьома новими наборами значень параметрів. Набір А, ((0.1,15,100), (0.1,15,100)), набір ((1,1,0.5), (1, 1,0.5)) і набір С ((20,1,100), (20,1,100) )).

Малюнок 1.4. Фазовий портрет із параметрами А

Як можна помітити, за будь-яких наборів параметрів, точка (0,0) - рівноважна, і до того ж стійка. Також на деяких малюнках можна побачити і сідлову точку.

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

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

Інтегральні криві відомі під іншими назвами, залежно від природи та інтерпретації диференціального рівняння або векторного поля. У фізиці інтегральні криві для електричного або магнітного поля відомі як лінії поля, а інтегральні криві для поля швидкостей рідини відомі як лінії струму. У динамічних системах інтегральні криві для диференціального рівняння називають траєкторіями.

Малюнок 1.5. Інтегральні криві

Рішення будь-якої системи можна розглядати і як рівняння інтегральних кривих. Очевидно, що кожна фазова траєкторія - проекція деякої інтегральної кривої в просторі x,y,tна фазову поверхню.

Для побудови інтегральних кривих є кілька способів.

Один з них – метод ізоклін. Ізокліна - це крива, що проходить через точки, в яких нахил функції буде завжди одним і тим же, незалежно від початкових умов (Hanski, 1999).

Він часто використовується як графічний метод вирішення звичайних диференціальних рівнянь. Наприклад, у рівнянні виду y"= f(x, y) ізокліни є лініями на площині (x, y), отриманими прирівнюванням f (x, y) до константи. Це дає ряд ліній (для різних констант), вздовж яких криві рішення мають один і той же градієнт, обчислюючи цей градієнт для кожної ізокліни, поле нахилу можна візуалізувати, що дозволяє порівняно легко намалювати наближені криві рішення.

Малюнок 1.6. Метод ізоклін

Цей метод не вимагає обчислень на комп'ютері, і був дуже популярним раніше. Зараз існують програмні рішення, які на комп'ютерах побудують інтегральні криві гранично точно та швидко. Однак, навіть такий метод ізоклін непогано зарекомендував себе як інструмент для дослідження поведінки рішень, оскільки дозволяє показати області типової поведінки інтегральних кривих.

Мальтузіанська модель зростання без насичення.

Почнемо з того, що, незважаючи на існування різних способів побудови, показати інтегральні криві системи рівнянь не так просто. Метод ізоклін, що згадується раніше, не підходить, оскільки він працює для диференціальних рівнянь першого порядку. А програмні засоби, що мають можливість побудови таких кривих, не знаходяться у відкритому доступі. Наприклад, Wolfram Mathematica, здатна цього, платна. Тому намагатимемося максимально використовувати можливості Wolfram Alpha, робота з якою описана в різних статтях та роботах (Orca, 2009). Навіть незважаючи на те, що картина буде явно не зовсім достовірною, але принаймні дозволить показати залежність у площинах (x, t), (y, t). Для початку розв'яжемо кожне з рівнянь щодо t. Тобто виведемо залежність кожної із змінних щодо часу. Для цієї системи отримуємо:

(1.10)

(1.11)

Рівняння симетричні, тому розглянемо лише одне їх, саме x(t). Нехай константа дорівнює 1. У разі скористаємося функцією побудови графіків.

Малюнок 1.7. Тривимірна модель для рівняння (1.10)

Мальтузіанська модель зростання з насичення.

Виконаємо аналогічні дії і для іншої моделі. Зрештою отримуємо два рівняння, що демонструють залежність змінних від часу.

(1.12)

(1.13)

Знову побудуємо тривимірну модель та лінії рівня.

Малюнок 1.8. Тривимірна модель для рівняння (1.12)

Оскільки значення змінних неотрицательны, то дроби при експоненті отримуємо негативне число. Таким чином згодом інтегральна крива зменшується.

Раніше було дано визначення системної динаміки для розуміння суті роботи, тепер зупинимося на цьому детальніше.

Системна динаміка - методологія та метод математичного моделювання для формування, розуміння та обговорення складних проблем, спочатку розроблена в 1950-х роках Джеєм Форрестером, та описана у його роботі (Forrester, 1961).

Системна динаміка одна із аспектів теорії систем як методу розуміння динамічного поведінки складних систем. Основою методу є визнання того, що структура будь-якої системи складається з численних відносин між її компонентами, які найчастіше настільки ж важливі визначення її поведінки, як і самі окремі компоненти. Прикладами є теорія хаосу та соціальна динаміка, описані у роботах різних авторів (Grebogi, 1987; Sontag, 1998; Кузнєцов, 2001; Табір, 2001). Також стверджується, що, оскільки у властивостях елементів часто не можуть бути знайдені властивості-цілого, у деяких випадках поведінка цілого не може бути пояснена з погляду поведінки частин.

Моделювання може по-справжньому показати всю практичну значущість динамічної системи. Хоча воно і можливе в електронних таблицях, існує безліч програмних пакетів, які були оптимізовані спеціально для цього.

Саме собою моделювання - це процес створення та аналізу прототипу фізичної моделі для прогнозування її продуктивності у світі. Імітаційне моделювання використовується, щоб допомогти проектувальникам та інженерам зрозуміти, за яких умов і в яких випадках процес може зазнати невдачі та які навантаження він може витримати (Хемді, 2007). Моделювання також допоможе передбачити поведінку потоків рідини та інші фізичні явища. У моделі аналізується приблизні умови роботи за рахунок застосовуваних імітаційних програмних засобів (Строгалев, 2008).

Обмеження можливості імітаційного моделювання мають загальну причину. Побудова та чисельний розрахунок точної моделі гарантує успіх лише в тих областях, де існує точна кількісна теорія, тобто коли відомі рівняння, що описують ті чи інші явища, і завдання полягає лише в тому, щоб вирішити ці рівняння з необхідною точністю. У тих областях, де кількісної теорії немає, побудова точної моделі має обмежену цінність (Базыкин, 2003).

Проте можливості моделювання не безмежні. Насамперед, це пов'язано з тим, що важко оцінити сферу застосування імітаційної моделі, зокрема, період часу, для якого прогноз може бути побудований з необхідною точністю (Law, 2006). Крім того, за своєю природою імітаційна модель прив'язана до конкретного об'єкта, а при спробі застосувати до іншого, навіть аналогічного йому об'єкта, вимагає радикального коригування або принаймні суттєвої модифікації.

Є загальна причина існування обмежень на імітаційне моделювання. Побудова і чисельний розрахунок «точної» моделі успішний лише за існуванні кількісної теорії, тобто у тому разі, якщо всі рівняння відомі, а завдання зводиться лише вирішення даних рівнянь з певною точністю (Базыкин, 2003).

Але навіть незважаючи на це, імітаційне моделювання – чудовий засіб візуалізації динамічних процесів, що дозволяє, за більш-менш правильної моделі, приймати рішення, засновані на її результатах.

У цій роботі моделі систем будуть побудовані за допомогою засобів системної динаміки, які пропонуються програмою AnyLogic.

Мальтузіанська модель зростання без насичення/

Перед побудовою моделі необхідно розглянути елементи системної динаміки, якими ми користуватимемося, та зв'язати їх із нашою системою. Наступні визначення були взяті із довідкової інформації програми AnyLogic.

Накопичувач – основний елемент діаграм системної динаміки. Вони застосовуються уявлення об'єктів реального світу, у яких накопичуються деякі ресурси: гроші, речовини, чисельності груп людей, деякі матеріальні об'єкти тощо. Накопичувачі відбивають статичний стан моделируемой системи, які значення змінюються згодом у відповідність до існуючими у системі потоками. Звідси випливає, що динаміку системи задають потоки. Вхідні та вихідні з накопичувача потоки збільшують або зменшують значення накопичувача.

Потік, як і вищезгаданий накопичувач, - основний елемент системно-динамічних діаграм.

Поки накопичувачі визначають статичну частину системи, потоки визначають швидкість зміни значень накопичувачів, тобто у часі відбуваються зміни запасів і, таким чином, визначають динаміку системи.

Агент може містити змінні. Змінні зазвичай використовуються для моделювання змінних характеристик агента або зберігання результатів роботи моделі. Зазвичай динамічні змінні складаються з функції накопичувачів.

Агент може мати параметри. Параметри часто використовуються для представлення деяких характеристик змодельованого об'єкта. Вони корисні, коли екземпляри об'єктів мають однакову поведінку, описану у класі, але відрізняються деякими значеннями параметрів. Між змінними та параметрами існує явна різниця. Змінна стан моделі і може змінюватися під час моделювання. Параметр зазвичай використовується для статичного опису об'єктів. Під час одного прогону моделі параметр зазвичай є константою і змінюється тільки тоді, коли потрібно переналаштувати поведінку моделі.

Зв'язок - елемент системної динаміки, що використовується для визначення залежності між елементами діаграми потоків і накопичувачів. Як приклад, якщо якийсь елемент A згадується в рівнянні або початковому значенні елемента B, спочатку необхідно з'єднати ці елементи зв'язком, що йде від A до B, і тільки потім ввести вираз у властивостях B.

Існують і деякі інші елементи системної динаміки, але вони не будуть задіяні в ході роботи, тому опустимо їх.

Для початку розглянемо з чого складатиметься модель системи (1.4).

По-перше, відразу відзначаємо два накопичувачі, які будуть містити в собі значення кількості продукції кожного з підприємств.

По-друге, оскільки у нас по два складові в кожному рівняння, то отримуємо по два потоки до кожного з накопичувачів, один вхідний, інший вихідний.

По-третє, переходимо до змінних та параметрів. Змінних лише дві. X та Y, відповідальні за зростання продукції. А також у нас є чотири параметри.

По-четверте, щодо зв'язків, кожен із потоків повинен бути пов'язаний зі змінними та параметрами, що входять до рівняння потоку, а також обидві змінні повинні мати зв'язок із накопичувачами для зміни значення з часом.

Детальний опис побудови моделі, як приклад роботи в середовищі моделювання AnyLogic, залишимо для наступної системи, оскільки вона дещо складніша і в ній використовується більше параметрів, і одразу перейдемо до розгляду готового варіанта системи.

Нижче на малюнку 1.9 представлена ​​побудована модель:

Малюнок 1.9. Модель системної динаміки для системи (1.4)

Усі елементи системної динаміки відповідають описаному вище, тобто. два накопичувачі, чотири потоки (два вхідних, два вихідних), чотири параметри, дві динамічні змінні, та необхідні зв'язки.

На малюнку видно, що більше продукції, тим сильніше її зростання, що призводить до різкого зростання кількості товарів, що відповідає нашій системі. Але, як вже було раніше сказано, відсутність обмеження на це зростання не дозволяють застосовувати дану модель на практиці.

Мальтузіанська модель зростання з насичення/

Розглядаючи цю систему, докладніше зупинимося на побудові моделі.

Першим кроком додаємо два накопичувачі, назвемо їх X_stock та Y_stock. Кожному з них поставимо початкове значення 1. Зазначимо, що відсутність потоків в класично заданому рівняння накопичувача нічого немає.

Малюнок 1.10. Побудова моделі системи (1.9)

Наступний крок – додавання потоків. Побудуємо для кожного накопичувача вхідний та вихідний потік за допомогою графічного редактора. Не можна забувати, що один із країв потоку повинен перебувати в накопичувачі, інакше вони не будуть пов'язані.

Можна помітити, що рівняння для накопичувача задалося автоматично, звичайно, користувач може і сам написати його, вибравши режим рівняння «довільний», але найпростіше залишити цю дію на програму.

Третім кроком у нас є додавання шести параметрів та двох динамічних змінних. Дамо кожному елементу ім'я у відповідність до його літерним виразом у системі, а також задамо початкові значення параметрів наступним чином: e1=e2=1, a12=a21=3, n1=n2=0,2.

Усі елементи рівнянь присутні, залишилося лише написати рівняння потоків, але цього спочатку необхідно додати зв'язки між елементами. Наприклад, вихідний потік, відповідальний за доданок , може бути пов'язані з e1 і x. А кожна динамічна змінна повинна бути пов'язана з відповідним їй накопичувачем (X_stock x, Y_stock y). Створення зв'язків відбувається аналогічно додавання потоків.

Після створення необхідних зв'язків можна переходити до написання рівнянь потоків, що демонструється на правому малюнку. Звичайно, можна піти і в зворотному порядку, але при існуванні зв'язків, під час написання рівнянь з'являються підказки для встановлення потрібних параметрів/змінних, що полегшує завдання в складних моделях.

Після виконання всіх кроків можна запускати імітаційну модель і подивитися на її результат.

Розглянувши системи нелінійних диференціальних рівнянь взаємодії підприємств за умов мутуалізму, можна зробити кілька висновків.

Існує два стани системи: різке необмежене зростання, або прагнення кількості продукції до нуля. Який із двох станів прийме система залежить від параметрів.

Жодна із запропонованих моделей, у тому числі модель з урахуванням насичення, не підходить для практичного застосування через відсутність ненульового сталого положення, а також причин, описаних у пункті 1.

У разі спроби подальшого дослідження даного типу симбіотичної взаємодії для створення моделі, що застосовується компаніями на практиці, необхідне подальше ускладнення системи та введення нових параметрів. Наприклад, Базикін у книзі наводить приклад динаміки двох мутуалістичних популяцій із запровадженням додаткового чинника внутрішньовидової конкуренції. За рахунок чого система набуває вигляду:

(1.15)

І в такому разі з'являється ненульове стійке положення системи, відокремлене від нульового «сідлом», що наближає її до реальної картини того, що відбувається.

2. Взаємодія підприємств в умовах протокооперації

Всі основні теоретичні відомості були представлені в попередньому розділі, тому при аналізі моделей, що розглядаються в цьому розділі, теорія здебільшого буде опущена, за винятком кількох моментів, з якими ми не стикалися у попередньому розділі, а також можливе скорочення обчислень. Розглянута в цьому розділі модель взаємодії організацій за умов протокооперації, що складається із систем двох рівнянь, заснованих на мальтузіанської моделі, виглядає як система (1.5). Проаналізовані в попередньому розділі системи показали, що для їх максимального наближення до діючих моделей необхідне ускладнення систем. Виходячи з даних висновків, відразу додамо на модель обмеження на зростання. На відміну від попереднього типу взаємодії, коли зростання, що не залежить від іншої компанії, негативний, у разі всі знаки позитивні, отже, маємо постійне зростання. Уникаючи недоліків, описаних раніше, постараємося обмежити його логістичним рівнянням, також відомим як рівняння Ферхюльста (Gershenfeld, 1999), що має такий вигляд:

, (2.1)

де P – чисельність популяції, r – параметр, що показує швидкість зростання, K – параметр, який відповідає за максимально можливу чисельність популяції. Тобто з часом чисельність популяції (у разі продукції), йтиме до певного параметра До.

Це рівняння допоможе стримати нестримне зростання продукції, яке ми спостерігали раніше. Таким чином система набуває наступного вигляду:

(2.2)

Не варто забувати, що обсяг товару, що зберігається на складі для кожної компанії різний, тому параметри, що обмежують зростання різні. Назвемо цю систему «», і надалі використовуватимемо цю назву, коли її розглядатимемо.

Другою системою, яку ми розглядатимемо, є подальший розвиток моделі з обмеженням Ферхюльста. Як і в попередньому розділі введемо обмеження на насичення, тоді система набуде вигляду:

(2.3)

Тепер кожен із доданків має власне обмеження, тому вже без подальшого аналізу можна помітити, що необмеженого зростання, як у моделях попереднього розділу, не буде. А оскільки кожен із доданків демонструє позитивне зростання, то й кількість продукції не впаде в нуль. Назвемо цю модель "модель протокооперації з двома обмеженнями".

Дані дві моделі розглядаються в різних джерелах про біологічні популяції. Тепер спробуємо дещо розширити системи. Для цього розглянемо наступний рисунок.

На малюнку продемонстровано приклад процесів двох компаній: сталеливарної та вугільної промисловості. В обох підприємствах є зростання продукції, що не залежить від іншої, а також є зростання продукції, яке виходить завдяки їхній взаємодії. Це ми вже враховували у ранніх моделях. Тепер варто звернути увагу, що фірми як виробляють продукцію, вони її ще й продають, наприклад, ринку чи взаємодіє з нею підприємства. Тобто. виходячи з логічних висновків, існує необхідність негативного зростання компаній за рахунок продажу продукції (на малюнку за це відповідають параметри 1 і 2), а також за рахунок передачі частини продукції іншому підприємству. Раніше ми враховували це лише з позитивним знаком у іншої компанії, але не розглядали те, що у першого підприємства під час передачі продукції її кількість зменшується. У такому разі отримуємо систему:

(2.4)

І якщо про доданок можна сказати, що якби в попередніх моделях було зазначено, що , характеризують природний приріст, а параметр може бути негативним, то різниці практично немає, то про доданок такого сказати не можна. До того ж надалі, розглядаючи подібну систему із запровадженим нею обмеженням, правильніше використовувати саме доданки позитивного і негативного зростання, оскільки у разі ними можуть накладатися різні обмеження, що неможливо для природного приросту. Назвемо її "розширена модель протокооперації".

І нарешті, четвертою моделлю, що розглядається, є розширена модель протокооперації з раніше згаданим логістичним обмеженням на зростання. І система для цієї моделі така:

, (2.5)

де - приріст продукції першого підприємства, який не залежить від другого, з урахуванням логістичного обмеження, - приріст продукції першої компанії, що залежить від другої, з урахуванням логістичного обмеження; - приріст продукції другого підприємства, що не залежить від першого, з урахуванням логістичного обмеження, - приріст продукції другої компанії, що залежить від першої, з урахуванням логістичного обмеження; - споживання товарів першого підприємства, не пов'язане з іншим; - споживання товарів другого підприємства, не пов'язане з іншим; - споживання товарів першої галузі другою галуззю; - споживання товарів другої галузі. першою галуззю.

Надалі ця модель позначатиметься як «розширена модель протооперації з логістичним обмеженням».

1 Стійкість систем у першому наближенні

Модель протокооперації з обмеженням Ферхюльста

Методи аналізу стійкості системи були зазначені в аналогічному розділі попереднього розділу. Насамперед знаходимо рівноважні точки. Одна з них, як завжди, нульова. Інша є точку з координатами .

Для нульової точки 1 = , 2 = , оскільки обидва параметри невід'ємні, то отримуємо нестійкий вузол.

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

Аналіз даної системи складніший за попередню в силу того, що додається фактор насичення, таким чином з'являються нові параметри, а при знаходженні рівноважних точок доведеться вирішувати не лінійне, а білінійне рівняння через змінну в знаменнику. Тому, як і попередньому випадку, залишимо визначення типу стійкості на фазові діаграми.

Незважаючи на появу нових параметрів, Якобіан у нульовій точці, так само, як і коріння характеристичного рівняння, виглядає аналогічно до попередньої моделі. Таким чином, у нульовій точці є нестійкий вузол.

Перейдемо до розширених моделей. Перша їх не містить жодних обмежень і набуває вигляду системи (2.4)

Зробимо заміну змінних, , і . Нова система:

(2.6)

У такому разі отримуємо дві рівноважні точки, точка А(0,0), В(). Точка лежить у першій чверті, оскільки змінні мають невід'ємне значення.

Для рівноважної точки А отримуємо:

. - Нестійкий вузол,

. - Сідло,

. - Сідло,

. - стійкий вузол,

У точці B коріння характеристичного рівняння є комплексними числами: λ1 = , λ2 = . Ми не можемо визначити тип стійкості покладаючись на теореми Ляпунова, тому проведемо чисельне моделювання, яке не покаже всі можливі стани, але дозволить дізнатися хоча б деякі з них.

Малюнок 2.2. Чисельне моделювання пошуку типу стійкості

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

Не вдаючись у подробиці обчислень, приходимо до наступних рівноважних точок. Точка А(0,0) та точка В з наступними координатами:

(), де a =

Для точки А визначення типу стійкості - тривіальна задача. Коріння характеристичного рівняння таке λ1 = , λ2 = . Таким чином отримуємо чотири варіанти:

1. λ1 > 0, λ2 > 0 - нестійкий вузол.

2. λ1< 0, λ2 >0 – сідло.

3. λ1 ​​> 0, λ2< 0 - седло.

4. λ1< 0, λ2 < 0 - устойчивый узел.

Говорячи про точку В, варто погодитися, що підстановка скорочень у вираз для неї ускладнить роботу з Якобіаном і знаходженням коренів характеристичного рівняння. Наприклад, після спроби їх пошуку за допомогою обчислювальних засобів WolframAlpha, виведення значення коренів зайняв близько п'яти рядків, що не дозволяє працювати з ними буквально. Звичайно, за наявності вже наявних параметрів, можна швидко знайти рівноважну точку, але це окремий випадокОскільки ми знайдемо стан рівноваги, якщо вона є, лише для даних параметрів, що не підходить для системи підтримки прийняття рішень, для якої модель і планується створюватися.

Через складність роботи з корінням характеристичного рівняння побудуємо взаємне розташування нуль-ізоклін за аналогією з розібраною в роботі Базикиною системою (Базикін, 2003). Це дозволить нам розглянути можливі стани системи, і надалі при побудові фазових портретів виявити рівноважні точки та типи їхньої стійкості.

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

(2.7)

Таким чином, ізокліни мають вигляд парабол.

Малюнок 2.3. Можливий варіант розташування нуль-ізоклін

Усього можливі чотири випадки їхнього взаємного розташування за кількістю загальних точок між параболами. Для кожного з них є свої набори параметрів, а значить і фазові портрети системи.

2 Фазові портрети систем

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

Як видно з наведених нижче малюнків, нульова точка - нестійкий вузол, а друга точка, якщо підставити числові значення параметрів, то отримаємо (-1.5, -1.5) - сідло.

Малюнок 2.4. Фазовий портрет для системи (2.2)

Таким чином, оскільки ніяких змін не повинно відбуватися, то для даної системи існують лише нестійкі стани, що пов'язано з можливістю необмеженого зростання.

Модель протокооперації із двома обмеженнями.

У цій системі є додатковий стримуючий фактор, тому фазові діаграми повинні відрізнятися від попереднього випадку, що і видно на малюнку. Нульова точка також - нестійкий вузол, але у цій системі з'являється стійке становище, саме стійкий вузол. При даних параметрах координати (5.5,5.5), він представлений малюнку.

Малюнок 2.5. Фазовий портрет для системи (2.3)

Таким чином, обмеження на кожен доданок дозволило отримати стійке положення системи.

Розширена модель протокооперації.

Побудуємо фазові портрети для розширеної моделі, але одразу використовуючи її модифікований вигляд:

Розглянемо чотири набори параметрів, причому таких, щоб розглянути всі випадки з нульовою рівноважною точкою, а також продемонструвати фазові діаграми чисельного моделювання, що використовується для ненульової рівноважної точки: набір А(1,0.5,0,5) відповідає стану , набір (1,0.5,-0.5) відповідає набір С(-1,0.5,0,5) а набір D(-1,0.5,-0,5) , тобто стійкого вузла в нульовій точці. Перші два набори продемонструють фазові портрети для параметрів, які ми розглядали при чисельному моделюванні.

Малюнок 2.6. Фазовий портрет системи (2.4) з параметрами А-D.

На малюнках необхідно звернути увагу на точки (-1,2) та (1,-2) відповідно, у них виникає «сідло». Для більш детального уявлення, малюнку представлений інший масштаб малюнка з сідловою точкою (1,-2). На малюнку в точках (1,2) та (-1,-2) видно стійкий центр. Що стосується нульової точки, то починаючи з малюнка по малюнок на фазових діаграмах чітко помітний нестійкий вузол, сідло, сідло та стійкий вузол.

Розширена модель протокооперації із логістичним обмеженням.

Як і в попередній моделі, продемонструємо фазові портрети для чотирьох випадків нульової точки, а також постараємося відзначити і ненульові рішення на цих діаграмах. Для цього візьмемо наступні набори параметрів із параметрами, вказаними в наступному порядку (): А(2,1,2,1), Б(2,1,1,2), С(1,2,2,1) та Д (1,2,1,2). Інші параметри для всіх наборів будуть наступними: , .

На поданих нижче малюнках можна спостерігати чотири рівноважні стани нульової точки, описаних у попередньому розділі для даної динамічної системи. А також на малюнках стійке положення крапки з однією ненульовою координатою.

Малюнок 2.7. Фазовий портрет системи (2.5) з параметрами А-B

3 Інтегральні траєкторії систем

Модель протокооперації з обмеженням Ферхюльста

Як і в попередньому розділі розв'яжемо кожне з диференціальних рівнянь окремо і явно виразимо залежність змінних від часового параметра.

(2.8)

(2.9)

З отриманих рівнянь видно, що значення кожної зі змінних зростає, що демонструється на тривимірної моделі нижче.

Малюнок 2.8. Тривимірна модель для рівняння (2.8)

Даний вид графіка на початку чимось нагадує тривимірне зображення мальтузіанської моделі без насичення, що розглядається в главі 1, оскільки має аналогічне їй швидке зростання, але надалі можна помітити зниження швидкості зростання через досягнення обмеження на обсяг продукції. Таким чином підсумковий зовнішній виглядінтегральних кривих схожий на графік логістичного рівняння, яке було використане, щоб обмежити одне з доданків.

Модель протокооперації із двома обмеженнями.

Вирішуємо кожне із рівнянь за допомогою засобів Wolfram Alpha. Таким чином, залежність функції x(t) зводиться до наступного виду:

(2.10)

Для другої функції ситуація аналогічна, тому опустимо її розв'язання. Чисельні значення виникли через заміну властивостей деякими відповідними їм значеннями, що не впливає на якісну поведінку інтегральних кривих. На наведених нижче малюнках помітно використання обмежень на зростання, оскільки згодом експоненційне зростання перетворюється на логарифмічний.

Малюнок 2.9. Тривимірна модель для рівняння (2.10)

Розширена модель протокооперації

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

(2.11)

(2.12)

Розширена модель протокооперації з логістичним обмеженням

Залежність x(t) виглядає так:

Без графіка складно оцінити, поведінка функції, тому скориставшись відомими нам засобами, побудуємо його.

Малюнок 2.10 Тривимірна модель рівняння

Значення функції зменшується при не малих значеннях іншої змінної, що пов'язано з відсутністю обмежень на негативний білінійний доданок, і є очевидним результатом

4 Системна динаміка взаємодіючих компаній

Модель протокооперації з обмеженням Ферхюльста.

Побудуємо систему (2.2). Використовуючи вже відомі нам інструменти, будуємо імітаційну модель. Цього разу на відміну від мутуалістичних моделей, у моделі буде логістичне обмеження.

Малюнок 2.11. Модель системної динаміки для системи (2.2)

Запустимо модель. У цій моделі варто відзначити той факт, що зростання взаємозв'язку нічим не обмежене, а зростання продукції без впливу іншого має специфічне обмеження. Якщо подивитися на вираз логістичної функції, можна помітити, що у разі, коли змінна (кількість товарів) перевищує максимально можливий обсяг зберігання, доданок стає негативним. У разі, коли існує лише логістична функція, таке неможливе, але при додатковому завжди позитивному факторі зростання таке можливо. І зараз важливо зрозуміти, що логістична функція впорається із ситуацією не надто швидкого зростання кількості продукції, наприклад, лінійного. Звернімо увагу на малюнки нижче.

Малюнок 2.12. Приклад роботи моделі системної динаміки системи (2.2)

На лівому малюнку продемонстровано 5 крок роботи програми відповідної запропонованої моделі. Але зараз варто звернути увагу на правий малюнок.

По-перше, для одного з вхідних потоків для Y_stock видалено зв'язок з х, виражена в доданку . Це зроблено для того, щоб показати різницю в роботі моделі при лінійному позитивному потоці, і білінійному зростанні, який представлений для X_stock. При лінійних необмежених потоках після перевищення параметра До система в якийсь момент приходить до рівноваги (у даній моделі рівноважний стан - 200 тисяч одиниць товару). Але набагато раніше білінійне зростання призводить до різкого зростання кількості товару, що переходить у нескінченність. Якщо ж залишити і правий і лівий постійно позитивні потоки білінійними, то вже приблизно на 20-30 кроці, значення накопичувача приходить до різниці двох нескінченностей.

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

Модель протокооперації із двома обмеженнями.

З'ясувавши недоліки минулої моделі і ввівши обмеження на другий доданок чинником насичення, побудуємо і запустимо нову модель.

Малюнок 2.13. Модель системної динаміки та приклад її роботи для системи (2.3)

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

Розширена модель протокооперації.

При розгляді системної динаміки цієї моделі буде продемонстровано можливості програмного середовища AnyLogic для яскравої візуалізації моделей. Усі попередні моделі були побудовані лише з використанням елементів системної динаміки. Тому самі моделі виглядали непомітно, вони не дозволяли відстежити динаміку зміни кількості продукції в часі та змінювати параметри під час роботи програми. При роботі з цією та наступною моделями постараємося скористатися ширшим спектром можливостей програми для зміни трьох зазначених вище недоліків.

По-перше, у програмі поряд із розділом «системна динаміка» у програмі також присутні розділи «картинки», «3D-об'єкти», що дозволяють урізноманітнити модель, що корисно за її подальшої презентації, оскільки робить вигляд моделі «приємнішою».

По-друге, для відстеження динаміки зміни значень моделі існує розділ статистика, який дозволяє додавати діаграми та різні інструменти збору даних, пов'язуючи їх з моделлю.

По-третє, для зміни параметрів та інших об'єктів під час виконання моделі є розділ «елементи керування». Об'єкти цього розділу дозволяють змінювати параметри під час роботи моделі (приклад, «бігунок»), вибирати різні стани об'єкта (приклад, «перемикач») і виконувати інші дії, що змінюють спочатку задані дані під час роботи.

Модель підходить для навчального знайомства з динамікою зміни продукції підприємств, але відсутність обмежень зростання не дозволяють використовувати її на практиці.

Розширена модель протокооперації із логістичним обмеженням.

Використовуючи вже готову попередню модель, додамо до неї параметри з логістичного рівняння для обмеження зростання.

Опустимо побудову моделі, оскільки на попередніх п'яти моделях, представлених у роботі, вже було продемонстровано всі необхідні інструментита принципи роботи з ними. Варто лише відзначити, що її поведінка схожа на модель протокооперації з обмеженням Ферхюльста. Тобто. відсутність насичення заважає її практичному застосуванню.

Після аналізу моделей за умов протокооперації визначимо кілька основних моментів:

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

Відповідні моделі повинні мати обмеження на кожному із доданків, оскільки в іншому випадку різке зростання білінійних множників «руйнує» всю імітаційну модель.

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

Висновок

В результаті проведеного дослідження було проведено аналіз шести систем, що описують динаміку виробництва продукції підприємствами, що взаємно впливають одна на одну. У результаті рівноважні точки і типи їх стійкості були визначені одним із наступних способів: аналітично, або завдяки побудованим фазовим портретам у випадках, коли аналітичне рішення з будь-яких причин неможливо. Для кожної із систем були побудовані фазові діаграми, а також побудовані тривимірні моделі, на яких при проектуванні можна отримати інтегральні криві в площинах (x, t), (y, t). Після використання середовища моделювання AnyLogic були побудовані всі моделі і розглянуті варіанти їх поведінки при певних параметрах.

Після проведення аналізу систем і побудови їх імітаційних моделей стає очевидним, що дані моделі можуть розглядатися лише як навчальні, або ж для опису макроскопічних систем, але ніяк не як система підтримки прийняття рішень для окремих компаній, через свою низьку точність і в деяких місцях не зовсім достовірне представлення процесів, що відбуваються. Але також не варто забувати, що якою б вірною не була описуюча модель динамічна система у кожної компанії/організації/галузі свої власні процеси та обмеження, таким чином створити та описати загальну модель неможливо. У кожному конкретному випадку вона видозмінюватиметься: ускладнюватиметься або навпаки спрощуватиметься для подальшої роботи.

Роблячи висновок з висновків до кожного розділу, варто загострити увагу на виявленому факті, що введення обмежень на кожне з доданків рівняння хоч і ускладнює систему, але дозволяє виявити стійкі положення системи, а також наблизити її до того, що відбувається насправді. І варто відзначити, що моделі протокооперації більше підходять для вивчення, оскільки мають ненульові стійкі положення на відміну від двох мутуалістичних моделей, які ми розглянули.

Таким чином, мету даного дослідження було досягнуто, а завдання виконано. У майбутньому як продовження даної роботи буде розглянуто розширену модель взаємодії виду протокооперації з трьома введеними на неї обмеженнями: логістичним, фактором насичення, нижньою критичною чисельністю, що дозволить створити більш точну модель для системи підтримки прийняття рішень, а також модель з трьома компаніями. Як розширення роботи можна розглянути і два інші типи взаємодії крім симбіозу, про які згадувалося в роботі.

Література

1. Bhatia Nam Parshad; Szegх Giorgio P. (2002). Стабільність теорії динамічних систем. Springer.

2. Blanchard P.; Devaney, R. L.; Hall, G. R. (2006). Differential Equations. London: Thompson. pp. 96-111.

Boeing, G. (2016). Visual Analysis of Nonlinear Dynamical Systems: Chaos, Fractals, Self-Similarity and Limits of Prediction. Systems. 4(4): 37.

4. Campbell, David K. (2004). Nonlinear physics: Fresh breather. Nature. 432 (7016): 455-456.

Elton C.S. (1968) reprint. Animal ecology. Велика Британія: William Clowes and Sons Ltd.

7. Forrester Jay W. (1961). Industrial Dynamics. MIT Press

8. Gandolfo, Giancarlo (1996). Economic Dynamics (Third ed.). Берлін: Springer. pp. 407-428.

9. Gershenfeld Neil A. (1999). The Nature of Mathematical Modeling. Cambridge, UK: Cambridge University Press.

10. Goodman M. (1989). Study Notes in System Dynamics. Pegasus.

Grebogi C, Ott E, та Yorke J. (1987). Chaos, Strange Attractors, і Фрактал Basin Boundaries в Nonlinear Dynamics. Science 238 (4827), pp 632-638.

12. Hairer Ernst; Nørsett Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York

Hanski I. (1999) Metapopulation Ecology. Oxford University Press, Oxford, pp. 43-46.

Hughes-Hallett Deborah; McCallum, William G.; Gleason, Andrew M. (2013). Calculus: Single and Multivariable (6 ed.). Джон wiley.

15. Libro J., Valls C. (2007). Global analytic перші integrals для реального планування Lotka-Volterra system, J. Math. Phys.

16. Jordan D.W.; Smith P. (2007). Non-Linear Ordinary Differential Equations: Introduction for Scientists and Engineers (4th ed.). Oxford University Press.

Khalil Hassan K. (2001). Nonlinear Systems. Prentice Hall.

Lamar University, Online Math Notes - Phase Plane, P. Dawkins.

Lamar University, Online Math Notes - Systems of Differential Equations, P. Dawkins.

Lang Serge (1972). Різні manifolds. Reading, Mass.-London-Don Mills, Ont.: Addison-Wesley Publishing Co., Inc.

Law Averill M. (2006). Simulation Modeling and Analysis with Expertfit Software. McGraw-Hill Science.

Lazard D. (2009). Тридцяти років з Polynomial System Solving, and now? Journal of Symbolic Computation. 44 (3): 222-231.

24. Lewis Mark D. (2000). Програми Dynamic Systems Approaches for Integrated Account of Human Development. Child Development. 71 (1): 36-43.

25. Malthus T.R. (1798). An Essay на Принципі популяції, в Оксфордському світі" Classics reprint. p 61, end of Chapter VII

26. Morecroft John (2007). Strategic Modelling and Business Dynamics: A Feedback Systems Approach. John Wiley & Sons.

27. Nolte D.D. (2015), Introduction to Modern Dynamics: Chaos, Networks, Space and Time, Oxford University Press.