Компьютерное моделирование динамических систем

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

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

2.1.1 Математическое моделирование динамических систем

При математическом моделировании исследование динамики работы оборудования осуществляется на основе изучения поведения его математической модели. Построение такой модели базируется на расчетной схеме исследуемого объекта, а также принятых ограничениях и допущениях. Большинство математических моделей динамических систем представляют собой систему дифференциальных уравнений, отражающих физические процессы в объекте [8].

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

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

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

=

Как известно, функция аналитическая, и ее разложение в ряд Тейлора выглядит так:

= х - + ? … = ? ()

При малых x. С увеличением x требуется учет второго, третьего и т.д. членов ряда, чтобы с заданной точностью аппроксимировать . Поэтому в случае мы получаем самую простую модель математического маятника:

+x = 0.

Следующим приближением будет модель нелинейного маятника:

+x ? = 0.

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

2.1.2 Численные методы решения систем дифференциальных уравнений

Рассмотрим методы решения обыкновенных дифференциальных уравнений (далее ОДУ) первого порядка с одним начальным условием. Это так называемая задача Коши: найти частное решение y = y(x) ОДУ y' = f(x,y) (4.1), удовлетворяющее начальному условию [11]:

y(x0) = y0 (4.2)

Эти методы легко обобщаются на системы уравнений первого порядка. Кроме того, уравнения более высоких порядков можно свести к системе ОДУ первого порядка. Например, изменяя в уравнении y"= ц(y',y,x) y на y1, а y' на y2, получим систему:

y' = y2,

y'2 = ц(y2,y1,x) (4.3)

Наиболее распространенными являются методы конечных разностей (МКР). В разностных методах вводится последовательность точек x0, x1,…xi и шаг hi = xi+1 - xi В каждой точке xi, называемой узлом вместо значений функции y(xi) вводят числа yi, которые аппроксимируют точное решение y на данном множестве точек. Функцию y, заданную в виде таблицы i=0,1,2,…, называют сеточной функцией.

Изменяя значение производной в уравнении (4.2) отношением конечных разностей, осуществляют переход от дифференциальной задачи (4.2), (4.3) к разностной как то сеточной функции y:

yi+1 = F(xi, hi, yi-r, yi-r+1,…, yi, yi+5). (4.4)

Конкретное выражение правой части (4.4) зависит от способа аппроксимации производной, и функция F определяет вычислительную схему метода.

Можно привести такую классификацию методов. Если r=0 и 0?s?1, то численный метод называется одношаговым, если r?1или s>1, то многошаговым. Многошаговые и одношаговые методы называются явными, если s=0, и неявными, если s?0.

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

 = ѓ(t, X, Y,…)

= g(t, X,Y,…)

которые имеют решение:

X = X(t),

Y = Y(t),

где t - независимая переменная (например, время); X, Y и т.д. - искомые функции (зависимые от t переменные). Функции f, g - заданы. Также предполагаются заданными и начальные условия, т.е. значения искомых функций в начальный момент.

Одно дифференциальное уравнение - частный случай системы с одним элементом. Поэтому, как правило, речь идет для определенности о системе уравнений.

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

Метод Рунге-Кутта заключается в рекуррентном применении следующих формул:

Xk+1 = Xk + (k1 + 2k2 + 2k3 + k4),

Yk+1 = Yk + (m1 + 2m2 + 2m3 + m4),…,

где

k1 = ѓ(tk, Xk, Yk,…)?t,

m1 = g(tk, Xk, Yk,…)?t,…,

k2 = ѓ(tk + , Xk +, Yk + ,…)?t,

m2 = g(tk + , Xk +, Yk + ,…)?t,…,

k3 = ѓ(tk + , Xk +, Yk + ,…)?t,

m3 = g(tk + , Xk +, Yk + ,…)?t,…,

k4 = ѓ(tk + ?t, Xk +k3, Yk + m3,…)?t,

m4 = g(tk + ?t, Xk +k3, Yk + m3,…)?t,…

2.2 Динамические системы типа хищник-жертва, обзор литературных источников по данной проблеме

Модель совместного существования двух биологических видов (популяций) типа «хищник - жертва» называется также моделью Вольтерра - Лотки. Групповой образ жизни хищников и их жертв радикально меняет поведение модели, придает ей повышенную устойчивость.

Данная модель была впервые получена А. Лоткой в 1925 году (она использовалась для описания динамики взаимодействующих биологических популяций).



Рис. 4 Модель химических реакций Лотки. Фазовый портрет системы при значениях параметров, соответствующих затухающим колебаниям.

В 1926 году независимо от Лотки аналогичные и более сложные модели были разработаны итальянским математиком В. Вольтерра. Его книга "Математическая теория борьбы за существование" стала классической для рассмотрения математических моделей взаимодействия видов. Книга, построена как математический трактат, в ней постулированы в математической форме свойства биологических объектов и их взаимодействий, а затем эти взаимодействия исследуются как математические объекты, Именно с этой работы В. Вольтерра начались современная математическая биология и математическая экология.

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

Пища либо имеется в неограниченном количестве, либо ее поступление с течением времени жестко регламентировано.

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

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

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

Если вид питается пищей, имеющейся в неограниченном количестве, прирост численности вида за единицу времени пропорционален численности вида.

Если вид питается пищей, имеющейся в ограниченном количестве, то его размножение регулируется скоростью потребления пищи, т.е. за единицу времени прирост пропорционален количеству съеденной пищи [15].

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

Допустим, что на некоторой территории обитают два вида животных: кролики (питающиеся растениями) и лисы (питающиеся кроликами), те и другие не иммигрируют, и не эмигрируют. Также допустим, что еды для травоядных животных у нас имеется с избытком. Пусть число кроликов x, число лис y. Используя Модель Мальтуса с необходимыми поправками, учитывающими поедание кроликов лисами, приходим к следующей системе, носящей имя модели Вольтерра - Лотки:

Уравнение изменения количества жертв примет вид:

= бx

где: б -- это коэффициент рождаемости жертв, х-- это величина популяции жертв, -- это скорость прироста популяции жертв.

Так как хищники стабильным питанием не обеспечены, то они вымирают. Следовательно, уравнение для хищников примет вид:

= ? вy

где: в -- это коэффициент убыли хищников, y -- это величина популяции хищников, -- это скорость прироста популяции хищников.

Встречи хищников и жертв (которые ? xy), убивают жертв с коэффициентом с и рождают новых хищников с коэффициентом d. С учётом этого, получаем систему уравнений:

= бx - cxy = (б ?cy)x

= ? вy + dxy = (?в + dx)y

Найдем стационарную точку, вокруг которой происходят колебания. Для стационарной позиции изменение популяции равно нулю. Следовательно:

б ? c = 0 (5.1)

?в + d = 0

Из чего следует, что:



Теперь нам надо ввести в нашу систему колебания и . Из-за малой величины квадратами, кубами и т.д. можно пренебречь. Теперь популяция x и y будет равняться:

x=

y =

Далее расписываем уравнение (5.1):

= б () - c () () =

=

Похожий ответ получаем относительно хищников:

После чего дифференцируем одно уравнение и подставляем в него другое:

,

,

что является уравнением гармонического осциллятора с периодом T=.

Классическая запись системы Лотки - Вольтерра:

, (5.2)

или, что тоже самое:

= бx - cxy

= dxy - вy,

где: х-- количество жертв, у-- количество хищников, t-- время, б, в, c и d --параметры, отражающие взаимодействия между видами.

Эта система имеет равновесное состояние, когда число кроликов и лис постоянно. Отклонение от этого состояния приводит к колебаниям численности кроликов и лис, аналогичным колебаниям гармонического осциллятора. Как и в случае гармонического осциллятора, это поведение не является структурно устойчивым: малое изменение модели (например, учитывающее ограниченность ресурсов, необходимых кроликам) может привести к качественному изменению поведения. Например, равновесное состояние может стать устойчивым, и колебания численности будут затухать. Возможна и противоположная ситуация, когда любое малое отклонение от положения равновесия приведет к катастрофическим последствиям, вплоть до полного вымирания одного из видов. На вопрос о том, какой из этих сценариев реализуется, модель Вольтерра - Лотки ответа не дает: здесь требуются дополнительные исследования.

На рис. 5 представлены фазовый портрет системы, по осям которого отложены численности жертв и хищников - (а) и кинетика численности обоих видов - зависимость численности от времени-(б). Видно, что численности хищников и жертв колеблются в противофазе.

Рис.5. Модель хищник-жертва Вольтерра, описывающая незатухающие колебания численности. А. Фазовый портрет. Б. Зависимость численности жертвы и хищника от времени.

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

Этот недостаток устранен в более реалистичных моделях. Модификация модели Вольтерра с учетом ограниченности субстрата в форме Моно:

R(S) = (5.3)

где КS - константа, равная концентрации субстрата, при которой скорость роста равна половине максимальной. 0 - максимальная скорость роста;

и самоограничения численности:

= rx (1 ? ), (5.4)

где K носит название "емкости популяции" и выражается в единицах численности (или концентрации), а r - собственная скорость роста, приводит к модели, подробно изученной А.Д.Базыкиным [12]:

= Ax ? -Ex2, (5.5)

= -Cy ? -My2,

Действительно, при малых численностях и в отсутствие хищника жертва (x) будет размножаться по экспоненциальному закону:

= Rx.

Хищник (y) в отсутствие жертв будут вымирать также по экспоненте. Если особей того или иного вида много, в соответствии с базовой моделью (5.4) срабатывает системный ферхюльстовский фактор (член -Ex2 в первом уравнении, и -My2 - во втором). Интенсивность взаимодействия видов считается пропорциональной произведению их численностей (как в модели (5.2)) и описывается в форме Моно (5.3). Здесь роль субстрата играет вид-жертва , а роль микроорганизмов - вид-хищник. Таким образом, модель (5.5) вобрала в себя свойства базовых моделей (5.2), (5.3), (5.4).

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

Представленные модели записаны в достаточно простой общепринятой форме для лучшего понимания. Рассмотрим системы уравнений представленные непосредственно в работе В. Вольтерра.

Запишем закон размножения:

(6.1)

где постоянное положительное число л - коэффициент рождаемости для рассматриваемой популяции.

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

Встречи хищников с жертвами пропорционально численности хищников и жертв бn1(t) n2(t). Эти встречи влияют отрицательно на увеличение популяции жертв (с коэффициентом -в1) и положительно на увеличение популяции (с коэффициентом в2).

Учитывая выше сказанное, запишем дифференциальные уравнения численности популяций n1(t) жертв и n2(t) хищников:

= л1n1?в1бn1n2, (6.2)

= ?л2n2?в2бn1n2.

Здесь величины n1 и n2 могут принимать только неотрицательные значения. При этом система (6.2) имеет следующие решения:

n1(t) ? 0, n2(t)=n2(0),

n2(t) ? 0, n1(t)=n1(0),

n1(t)? n2(t)?0.

Таким образом, фазовым пространством системы (6.2) является первый квадрант плоскости (см. рис. 5).

Покажем, что система (6.2) имеет первый интеграл:

V(n1,n2)=л2.

Действительно,

V(n1(t),n2(t))=л2.=

=

Отсюда следует, что и функция W(n1,n2)= =g(n1)h(n2), где g(n1)=, h(n2)=, так же является первым интегралом системы (6.2). Легко заметить, что функции g и h имеют максимумы соответственно в точках

, (6.3)

Таким образом, функция W(n1,n2) обладает следующим свойством:

W(n1,n2)? W(,)

Отсюда следует, что линии функции W(n1,n2) являются замкнутыми кривыми, окружающими точку (6.3). Таким образом, траектории системы (6.2) расположены на замкнутых линиях уровня и точка (6.3) является состоянием равновесия[15,16].

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

= (л1 ?)n1?в1бn1n2, (6.4)

= ?(л2?)n2?в2бn1n2,.

Здесь и - коэффициент интенсивности истребления.

Рассмотрим сначала случай л1 <(истребление жертв идет интенсивнее их размножения при отсутствии хищников). Здесь обе популяции обречены на вымирание.

В самом деле, рассмотрим функцию V(n1,n2)=. Для решений n1(t) и n2(t) уравнений (6.4) получим

(n1(t),n2(t))= (л1 ?)? (л2 ?)??еV(n1(t),n2(t)), (6.5)

где е=min().

Переписав неравенство (6.5) в виде

(V(n1(t),n2(t)))*? 0 (6.6)

и проинтегрировать обе части соотношения (6.6) от 0 до t, получим оценку

V(n1(t),n2(t) V(n1(0),n2(0)). (6.7)

Из оценки (6.7) следует , .

Последнее означает истребление жертв и вымирание хищников.

Далее рассмотрим наиболее интересный случай л1 >. Проинтегрировав систему (6.2) от 0 до Т получим:

 ,

.

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

Раздел 3. Компьютерное моделирование динамических систем

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

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

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

Несмотря на то, что современные универсальные пакеты визуального моделирования обладают рядом общих свойств (позволяют строить из блоков иерархические функциональные схемы, предоставляют пользователю схожие библиотеки численных методов, средства визуализации поведения и наборы анимационных возможностей, поддерживают технологию объектно-ориентированного моделирования), все же можно их разделить на три основные группы (схема 2):

пакеты, использующие язык блочного моделирования;

пакеты, использующие язык физического моделирования;

пакеты, ориентированные на использование схемы

гибридного автомата.

Схема 2.

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

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

Наиболее известными представителями первой группы являются:

подсистема Simulink пакета MATLAB;

пакет EASY5;

подсистема SystemBuild пакета MATRIXx ;

VisSim.

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

Среди пакетов, принадлежащих ко второй группе, можно отметить:

Dymola;

Omola и OmSim;

Smile;

Modelica.

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

К этой группе относятся:

пакет Shift;

пакет Model Vision Studium$

пакет WinSet.

3.2 Описание программного пакета WInSet

3.2.1 Окно программы WInSet

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

Рис. 6. Главное окно WinSet.

Для запуска программы WInSet достаточно вызвать главное меню Windows нажатием кнопки Пуск, перейти в меню Программы\WInSet и выбрать пункт WInSet. После этого на экране появится главное окно программы (рис. 6). Оно состоит из нескольких компонентов: главное меню, панель инструментов, рабочая область, панель координат, строка состояния.

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

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

Рис. 7. Строка главного меню и панель инструментов WinSet.

Панель координат точки находится справа от рабочей области. В нижней части панели отображаются значения координат точки на плоскости, соответствующей текущему положению указателя мыши. В верхней части панели отображаются координаты последней вычисленной точки (либо начальной точки, заданной по умолчанию). Эти значения используются в качестве начальных условий для построения, если оно запущено без указания начальной точки. Во время выполнения вычислений на панели координат отображается полоска индикатора, показывающая степень завершения процесса. Для того чтобы скрыть панель координат точки, нужно снять отметку с пункта меню Вид\Координаты точки. Над нижней кромкой окна расположена строка состояния, в которой отображается информация о работе программы. Когда программа выполняет какую-либо длительную операцию, в левой части строки состояния появляется текст пояснения [17].

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

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

а)

б)

Рис. 8. Работа с изображением: а) до выделения области, б) после выделения.

Все действия, выполняемые с помощью мыши, можно выполнить и с клавиатуры. Доступ к пунктам главного меню осуществляется с помощью комбинации клавиш Аlt+<подчеркнутая буква пункта меню>. Переход между элементами диалогового окна происходит с помощью клавиши табуляции Tab (Shift+Tab -- в обратном направлении). В диалоговых окнах истречаются поля ввода, содержащие выпадающий список для выбора значений. В таком поле для перемещения по списку используйте клавиши со стрелками вверх и вниз.

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

Вот некоторые из них:

Открыть файл слайда <Ctrl+O>

Сохранить построение в файле слайда WinSet <Ctrl+L>

Сохранить построение в графическом формате <Ctrl+S>

Загрузить графическую палитру из файла <Ctrl+T>

Распечатать построение <Ctrl+P>

Отменить последний запуск построения <Alt+BkSp>

Отменить предыдущее масштабирование <Ctrl+BkSp>

Задать начальные условия <F4>

Редактировать параметры системы <F7>

Изменить текущий цвет построения <F6>

Перерисовать изображение <F5>

Перерисовать оси координат <Shif t+F5>

Очистить построение <F8>

Выбрать систему <F3>

Добавить новую систему <Shif t+F3>

Запустить вычисление <Enter>

Прервать построение <Esc>

Редактировать установки построения <F9>

Редактировать установки вычислений <F11>

3.2.2 Приемы работы с графическими построениями для отображений и систем ОДУ

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

Рис. 9. Окно установок построения.

С помощью выпадающих списков для каждой из осей координат назначается переменная системы. Граничные значения по осям задаются с помощью полей ввода максимального и минимального значений. Расстановка отметок по каждой из осей координат управляется комбинированным выпадающим списком Отметки. Возможны следующие варианты: Auto -- автоматическая расстановка, None -- не делать отметок по осям. Если в этом поле ввести числовое значение, то интервал между отметками будет задан явно.

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

Рис. 10. Окно установок построения для трехмерных графиков.

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

Для всех систем ОДУ определен тип построения Traject - фазовые траектории и графики решений. Для некоторых типов систем ОДУ возможно построение отображения Пуанкаре. Для этих систем определен тип построения Map. Для итеративных отображений -- только Map.

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

Рис. 11. Окно редактирования начальной точки.

Для запуска построения из текущей точки можно использовать клавишу Enter или пункт меню «Запуск». Это приведет к выполнению одного этапа построения. После окончания вычислений в качестве начальной точки фиксируется последняя вычисленная точка.

Если после завершения одного этапа построения снова нажать «Запуск» или клавишу Enter, не меняя начальную точку, то второй этап построения будет продолжением первого. Это важно, например, при построении сепаратрис двумерных отображений а также траекторий неоднозначных отображений.

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

Рис. 12. Окно установок счета.

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

Для того чтобы отменить последний выполненный этап построения, нужно использовать пункт меню Правка\Отмена построения или нажать комбинацию клавиш <Alt> + <BackSpace>. Построение будет перерисовано без последнего этапа.

Кроме случая отмены предыдущего запуска, активное построение перерисовывается в следующих ситуациях:

При изменении размеров окна программы, например при развертывании на весь экран.

При переназначении переменных системы по осям координат с помощью окна Установки построения (рис. 10).

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

Для того чтобы сохранить построение в файле слайда WInSet, нужно выбрать пункт меню Файл\Сохранить слайд... и в появившемся окне сохранения файла ввести имя файла слайда. Файлы слайдов WInSet имеют тип .ini. Чтобы сохранить образ рабочей области в файле графического формата, нужно выбрать пункт меню Файл\Сохранить в графическом формате.... Появится окно сохранения графического файла (рис. 13), в котором требуется выбрать тип файла: растровый формат BMP или векторный формат EPS, и ввести имя файла.

Рис. 12. Окно сохранения файла.

Иногда необходимо выполнить несколько последовательных увеличений фрагментов построения и сохранить получившийся ряд изображений в графических файлах. При этом на каждом изображении обычно требуется выделить в виде прямоугольника область, которая будет увеличена на следующем изображении. Если при увеличении фрагмента с помощью мыши удерживать нажатой клавишу <Shift>, то непосредственно перед перерисовкой изображения в новом масштабе будет выведено окно сохранения построения в графическом файле. Это даст возможность указать имя и тип файла, в котором будет сохранен образ рабочей области с отмеченным для увеличения фрагментом.

3.2.3 Сценарий работы

Непосредственно после запуска программы рабочая область пуста, и некоторые пункты меню неактивны. Необходимо выбрать объект исследования -- динамическую систему. Выберем пункт главного меню Система\Выбор систем.

Открылось диалоговое окно выбора системы, показанное на на рис. 13. В списке категорий выберем Diff Eq, а из списка систем этой категории ? строчку Volterra. Нажмем кнопку «ОК». Система выбрана. В заголовке окна программы появилось название системы, а в рабочей области ? оси координат.

Рис. 13. Окно выбора системы.

После выбора системы автоматически вызывается диалоговое окно Параметры системы (рис. 14). В нем отображается запись уравнений системы и предоставляется возможность ввести значения изменяемых параметров. Оставим их без изменения и нажмем «ОК».

Для того чтобы запустить построение системы Вольтерра, выберем пункт главного меню «Запуск» или нажмем на клавиатуре клавишу Enter.

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

Рис. 14. Окно параметров системы.

Рис. 15. Модель хищник-жертва с стандартными параметрами WInSet.

Если выбрать пункт меню Правка\Параметры системы и ввести новое значение параметра, например р2=0,1, то внешний вид модели изменится, хотя задача остается той же (рис.16).

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

Выбрав в меню Система\Добавить новую..., откроем диалоговое окно Ввод новой системы (рис. 17а). Оно содержит две вкладки: «Общее» и «Текст».

Рис. 16. Модель хищник-жертва после изменения параметра р2=0,001 на р2=0,1.

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



Рис. 17а Диалоговое окно ввода новой системы.

Назовем систему «система 1» и выберем для нее тип Diff Eq (дифференциальные уравнения). Число 2 в поле Размерность означает, что построение будет задано системой из двух дифференциальных уравнений.

Система Вольтерра - Лотки для WinSet будет иметь вид:

x'1 = p1*x1 - p2*x1*x2

x'2 = -p3*x2+p4*x1*x2.

После окончания ввода уравнений нажмем кнопку «Компилировать текст» и кнопку «ОК» (рис.17б).



Рис. 17б. Окно ввода новой системы, вкладка Текст.

На рис. 18 видно, что теперь в окне выбора системы появилась заданная нами система.

Рис. 18. Выбор системы, заданной пользователем.

Нажав «ОК» мы попадем в окно параметров новой системы. Примем, что p1 = p2 = p3 = p4 = 1 (рис.19а).



Рис. 19а. Ввод данных в окно Параметры системы.

Исходя из предмета модели (численность популяций) необходимо, чтобы построение располагалось в первом квадранте фазовой плоскости, то есть имела неотрицательные значения. Для задания этого условия выберем пункт меню Установки - Построение и в полях минимального значения поставим 0. Затем снова нажмем «ОК». Теперь в рабочей области отображается только первый квадрант (рис.20).

Рис. 19б. Установка минимальных значений.



Рис. 20. Рабочая область программы после изменения установок.

Записав систему как

и, применив метод линеаризации, получим:

тогда

a11 = 0 a12 = ?1 и - матрица линеаризации системы.

a21 = 1 a22 = 0

Найдем детерминант:

л2+1 = 0, л1,2 =

,

Проинтегрируем полученную систему уравнений от 0 до T.

0 =л1T ? ,

0 =?л2T ?

C( ) C(1,1) -точка равновесия.

Колебания в малой окрестности C(1,1) близки к гармоническим с периодом 2р. По мере удаления от положения равновесия период будет возрастать.

Продемонстрируем это явление наглядно. Для этого зафиксируем численность популяции хищников x2 = 1. Затем, используя диалоговое окно Правка - Начальная точка и нажимая кнопку <Запуск>, будем двигаться по оси x1 от 1 к 0 с шагом 0,1 и временем интегрирования t = 6,28 (рис. 21а, 21б).



Рис. 21а. Ввод начальной точки построения.

Рис. 21б. Результаты построения при x2 = 1, x1 = 0,9;0,8;…;0,2.

Достроим траектории (рис.22) и представим изменение периода Т в виде диаграммы.



Рис. 22. Изменение периода восстановления исходной численности популяции жертв.

Диаграмма 1.

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

Рассмотрим построение в WinSet еще нескольких моделей типа хищник-жертва.

Модель А.Д.Базыкина (формула 5.5) для WinSet будет иметь вид (рис.23):

x'1 = p1*x1 - p2*x1*x2/(1+p3*x1)?p4*x1^2

x'2 = ?p5*x2?p6*x1*x2/(1+ p3*x1)?p7*x2^2.

Рис.23. Модель хищник-жертва в интерпретации А.Д.Базыкина.

Модель хищник-жертва с насыщением (рис.24), рассмотренная В.А. Треногиным в работе «Обыкновенные дифференциальные уравнения» [18]:

x'1 = p1*p2*x1 - p3*x1*x2/(1+p4*x1)

x'2 = p5* x1*x2/(1+ p6*x1)?p5*p6*x2.

Рис.24. Модель хищник-жертва с насыщением.

Модель Колмогорова-Вольтерра (рис. 25):

x'1 = x1*(1 - x1^2)/(1 + x1^2)- p1*x2*x1/(1 + x1^2)

x'2 = x2*(p2* x2^2 - 1)/(x1^2 + 1)

Рис.25 Модель Колмогорова-Вольтерра.

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

3.3 Описание программного пакета Mathlab

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

После запуска Mathlab на экране появляется основное окно системы, показанное на рис. 26.

Рис.26 Окно программы Mathlab.

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

Остальное пространство окна занято четырьмя областями, которые можно переместить или закрыть по усмотрению пользователя. Окно Current Directory (текущей папки) отображает папку, в которой расположены системные файлы Mathlab. В ту же папку будут по умолчанию сохранены созданные файлы и файлы ошибок.

Слева от Current Directory расположено окно Command Window (окно команд) со списком выбора функций (рис.27).

Рис.26 Окно Command Window.

В командном режиме действует простейший строчный редактор. Для редактирования используются следующие клавиши:

Перемещение курсора вправо на один символ <Ctrl+>

Перемещение курсора влево на один символ <Ctrl+f>

Перемещение курсора вправо на одно слово <Ctrl +r>

Перемещение курсора влево на одно слово <Ctrl +1>

Перемещение курсора в начало строки <Home или Ctrl+a>

Перемещение курсора в конец строки <End или Ctrl +e>

Перелистывание предыдущих команд вверх или вниз для подстановки в строку ввода <Ctrl+рu Ctrl+n>

Стирание символа справа от курсора <Del или Ctrl+d>

Стирание символа слева от курсора <Ctrl+h>

Стирание до конца строки <Ctrl+k>

Очистка строки ввода <Esc>

Включение/выключение режима вставки <Ins>

Перелистывание страниц сессии вверх <PgUp>

Перелистывание страниц сессии вниз <PgDn>

Эти возможности кажутся примитивными, но позволяют пользователю работать быстрее. Они обеспечивают важное свойство новых версий систем -- их совместимость со старыми версиями в части преемственности навыков работы.

Справа расположена рабочая область для просмотра и импорта данных и окно Command History, в котором отображаются последние проведенные операции с указанием даты и времени их выполнения (рис. 27).



Рис.27 Окно Command History и рабочая область.

Как упоминалось выше, любое из этих окон может быть закрыто, или перемещено для удобства работы. Закрытые панели можно восстановить, выбрав соответствующий пункт в меню Window.

3.3.1 Типы данных Mathlab

Число -- простейший объект языка Mathlab, представляющий количественные данные. Числа можно считать константами, имена которых совпадают с их значениями. Числа используются в общепринятом представлении о них. Они могут быть целыми, дробными, с фиксированной и плавающей точкой. Возможно, представление чисел в хорошо известном научном формате с указанием мантиссы и порядка числа. Ниже приводятся примеры представления чисел:

0

2

-3

2.301 0.00001 123.45бе-24

-234.456е10

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

Константа -- это предварительно определенное числовое или символьное значение, представленное уникальным именем. Числа (например 1, -2 и 1.23) являются безымянными числовыми константами.

Другие виды констант в Mathlab принято назвать системными переменными, поскольку, с одной стороны, они задаются системой при ее загрузке, а с другой -- могут переопределяться. Основные системные переменные, применяемые в системе Mathlab, указаны ниже:

i или j -- мнимая единица (корень квадратный из -1);

pi - число п - 3.1415926...;

eps -- погрешность операций над числами с плавающей точкой (2-52);

realmin -- наименьшее число с плавающей точкой (2-1022);

realmax -- наибольшее число с плавающей точкой (21023);

inf -- значение машинной бесконечности;

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

NaN -- указание на нечисловой характер данных (Not-a-Number).

Поскольку Mathlab используется для достаточно сложных вычислений, важное значение имеет наглядность их описания. Она достигается, в частности, с помощью текстовых комментариев. Текстовые комментарии вводятся с помощью символа %.

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

=: Имя_переменной - Выражение

Типы переменных заранее не декларируются. Они определяются выражением, значение которого присваивается переменной. Так, если это выражение -- вектор или матрица, то переменная будет векторной или матричной. Имя переменной (ее идентификатор) может содержать сколько угодно символов, но запоминается и идентифицируется только 31 начальный символ. Имя любой переменной не должно совпадать с именами других переменных, функций и процедур системы, т. е. оно должно быть уникальным. Имя должно начинаться с буквы, может содержать буквы, цифры и символ подчеркивания _. Недопустимо включать в имена переменных пробелы и специальные знаки, например +,.-, *, / и т. д., поскольку в этом случае правильная интерпретация выражений становится невозможной.

Оператор -- это специальное обозначение для определенной операции над данными -- операндами. Например, простейшими арифметическими операторами являются знаки суммы +, вычитания -, умножения * и деления /. Операторы используются совместно с операндами. Например, в выражении 2+3 знак + является оператором сложения, а числа 2 и 3 -- операндами.

3.3.2 Построение системы хищник-жертва

При построении моделей в Mathlab R2009 нам приходится столкнуться с несколькими проблемами. Начнем с того, что программа довольно громоздкая (2,17 Gb) и поставляется «в сборе» с полным комплектом прикладных пакетов, что существенно замедляет загрузку Mathlab и работу с ним. Неприятным моментом также является проблема совместимости с операционной системой Windows 7 - через случайные промежутки времени происходит крах системы, кроме того полученные «скриншоты» достаточно плохого качества из-за неудачного дизайна рабочих окон. Поэтому в работе использованы изображения, взятые из руководств по использованию программы и результат построения, обработанный средствами Adobe Photoshop.

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

Со временем число карасей и щук меняется, но так как рыбы в пруду много, то будем считать х и у непрерывными функциями времени t. Будем называть пару чисел (х, у) состоянием модели. Попробуем найти, как меняется состояние (х, у). Рассмотрим dx/dt - скорость изменения численности карасей. Если щук нет, то число карасей увеличивается и тем быстрее, чем больше карасей. Будем считать, что эта зависимость линейная : dx/dt ~ a1 x, причем коэффициент a1 зависит только от условий жизни карасей, их естественной смертности и рождаемости. Скорость изменения dy/dt числа щук (если нет карасей), зависит от числа щук y. Будем считать, что dy/dt ~ -a2y . Если карасей нет, то число щук уменьшается (у них нет пищи) и они вымирают. В экосистеме скорость изменения численности каждого вида также будем считать пропорциональной его численности, но только с коэффициентом, который зависит от численности особей другого вида. Так, для карасей этот коэффициент уменьшается с увеличением числа щук, а для щук увеличивается с увеличением числа карасей. Будем считать эту зависимость также линейной. Тогда получим систему из двух дифференциальных уравнений:

= a1x ? b1yx

= ? a2y + b2yx

Числовые коэффициенты a1, a2, b1, b2 - называются параметрами модели. Очевидно, что характер изменения состояния (x, y) определяется значениями параметров. Изменяя эти параметры и решая систему уравнений модели, можно исследовать закономерности изменения состояния экологической системы.

Реализация данной задачи средствами Mathlab производится с использованием прикладного пакета System Identification Toolbox. Это пакет расширения Mathlab, содержащий инструменты создания математических моделей динамических систем на основе наблюдаемых входных/выходных данных. Пакет снабжен гибким графическим интерфейсом, помогающим организовывать данные и создавать модели. Методы идентификации, входящие в пакет, применимы для широкого класс задач от проектирования систем управления и обработки сигналов до анализа временных рядов и вибраций.

На рис. 28 и 29 модель хищник-жертва, построенная средствами программного пакета Mathlab (на рис.29 добавлена ось времени t).



Рис.28 Реализация модели хищник-жертва в программе Mathlab

Рис.29 Реализация модели хищник-жертва c учетом времени t

Заключение

Таким образом, мы рассмотрели особенности построения динамических систем с использованием программных пакетов WInSeT и Mathlab, и проанализировали особенности математического моделирования средствами этих пакетов. Дополнительно были изучены литературные источники по данной проблеме и приведены классификация моделей и способов моделирования

Так же было произведено моделирование динамических систем типа хищник-жертва и приведен анализ полученных моделей.

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

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

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

Библиография

1. А.В. Морозов, И.А. Бригаднов - “Математические основы теории систем. Динамические системы” - Спб, СЗТУ,

2. Бенькович Е.С., Колесов Ю.Б., Сениченков Ю.Б. - “Практическое моделирование динамических систем” - Спб, БХВ-Петербург, 2002.

3. Ю.Б. Колесов, Ю.Б. Сениченков. - “Имитационное моделирование сложных динамических систем” - СПб.: БХВ-Петербург, 2002.