---
title: "Модели популяций"
author: "Varykhin A."
date: "22 11 2020"
output:
  html_document: default
  word_document: default
  pdf_document: default
---

```{r setup, include=FALSE,error=FALSE,message =FALSE,dpi=150}
knitr::opts_chunk$set(echo = TRUE)
```


















# **КОМПЬЮТЕРНЫЙ ПРАКТИКУМ ПО КУРСУ**
## "Системный анализ и основы компьютерного моделирования экосистем"











#### Составил: к.э.н. Варюхин А.М. 

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













#### Саратов 2020

* * *

## МЕТОДИЧЕСКИЕ УКАЗАНИЯ И ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТ

### *"Исследование моделей роста популяций"*

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

**Введение**

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

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

После ознакомления с теорией предлагается изучить несколько типов математических моделей. Логика работы с каждой из моделей следующая, пользователь строит модель, задает соответствующие параметры модели и получает в графическом представлении различные типы поведения системы.
Наборы параметров, обеспечивающие характерные режимы динамики численности приведены в заданиях, приложенных к каждой из моделей. После выполнения очередного задания пользователь отвечает на контрольные вопросы.
Модельные параметры, которые должны задаваться пользователем, следующие: начальная численность популяции, скорость роста популяции, емкость среды, коэффициент смертности, коэффициент внутривидовой конкуренции и.т.д..Работы основывается на последовательном выполнении заданий, в каждом из которых изучается влияние одного или группы параметров на режим динамики численности популяции. 
В качестве программного инструментария предлогается использовать язык
**R**, представляющий собой  свободную программную среду вычислений с открытым исходным кодом в рамках проекта $GNU$. Язык программирования     **R **довольно прост. Его наиболее сильная сторона- возможность неограниченного расширения с помощью пакетов. В базовую поставку **R** включен основной набор пакетов, а всего по состоянию на 2020 год доступно более 15 тысяч пакетов. В **R** реализованы практически все актуальные средства универсальных статистических вычислений, такие как регрессионный анализ и анализ временных рядов, а также множество специфических алгоритмов для решения узкоспециализированных задач и исследований в отдельных областях. Ещё одна особенность языка — возможность создания качественной графики типографского уровня, которая может быть экспортирована в распространённые графические форматы и использована для презентаций или публикаций. Все это обусловило выбор языка **R** как базового программного инструментария для изучения математических моделей развития и взаимодействия популяций.
Для выполнения работ небходимо установить на свой компьютер актуальную версию языка **R ** с $guiR$ (он устанавливается автоматически), $RStudio$ и некоторые допонительные пакеты. Изучить основы программной среды **R**, уделив особое внимание пакету $RMarkdown$, который должен быть использован при создании отчетов о проделанной работе.




* * *


### *Лабораторная работа 1 "Модель экспоненциального роста"*.

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

$$\Delta X = r X \Delta t $$

Здесь $r$- константа собственной скорости роста популяции. Совершив предельный переход, получим линейное дифференциальное уравнение:

$$X(t) = X_0 exp(r t)$$

где $X_0$ - начальная численность популяции.
Примером применения модели Мальтуса может служить описание развития однородной популяции в условиях неограниченных ресурсов питания (рост клеточной культуры до начала истощения питательной среды)

#### Порядок выполнения работы 
 
 В данном задании исследуют влияние начальной численности и параметра r на динамику роста популяции.
**Цель задания** - изучить влияние начальной численности и скорости роста на характер кривой динамики численности и освоить общие принципы работы с программными средствами R и RMarkdown

* * *

#### ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ

**скрипт 1**

```{r,message =FALSE,fig.height=4,fig.width=5}
#========ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ=================

library(deSolve) # решение дифф. уравнений с начальными условиями
library(ggplot2)
library(ggthemes)

p_r <<-0.2
p_y0 <<-4

p_t0 <<-0
p_dt <<-0.1
p_tn <<-20

f1<- function (t,y, parms ){
dY.dt <-y [1]*r
return (list (dY.dt ))}

t0<- seq (p_t0 ,p_tn ,p_dt)
r <-p_r
y0 <-p_y0

out <- ode (y=y0,t=t0,func =f1,parms = NULL )
out <- data.frame(out)
#head(out)

g40 <-ggplot(data=out, aes(x=out[,1],y=out[,2]))
g40 <-g40 + geom_line(data=out,aes(x=out[,1],y=out[,2]),
colour="red",linetype=1,size=0.8)
g40 <-g40 +labs(x = "t",y="Y",title = "ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ")
#g40 <-g40 + theme_stata() + scale_colour_stata()

g40

```


* * *

1.Скопируйте  к себе на компьютер, выше привенный скрипт.

Запуститу его в $guiR$ или $Rstudio$. Задайте $r$ – скорости роста популяции $(r < 0,5)$ и начальное значение численности $X_0$ $(X_0 < 5)$. Задайте интервал моделирования 20. Запустите скрипт на исполнение. Посмотрите на результаты. Поэкспериментируйте с моделью, меняйте параметры, и посмотрите, как меняются графики. Предложенный выше скрипт, позволяют анализировать отклик модели на задаваемые вручную параметры $r$ и $X_0$.

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

#### CКРИПТ С АВТОМАТИЧЕСКИМ ЗАДАНИЕМ ЗНАЧЕНИЙ ПАРАМЕТРА r

**скрипт 2**

```{r,fig.height= 4,fig.width=5}

#========ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ 2=================

library(deSolve) # решение дифф. уравнений с начальными условиями
library(ggplot2)
library(ggthemes)

n_r <-3
p_r0 <-0.5
d_p_r0 <-0.2
rtext <-c('r =','r =','r =')
p_y0 <<-50
y0text <-paste('y0 =',as.character(p_y0))
p_t0 <<-0
p_dt <<-0.1
p_tn <<-7

t0<- seq (p_t0 ,p_tn ,p_dt)
r <-p_r0
y0 <-p_y0
n_j <-as.integer((p_tn/p_dt)+1)
dat_r <-matrix(1:n_j*n_r, nrow = n_j, ncol =n_r)
 

for(i in 1:n_r){
r <-r+d_p_r0
rt <-as.character(r)
rtext[i] <-paste(rtext[i],rt)

f1<- function (t,y, parms ){
dY.dt <-y [1]*r
return (list (dY.dt ))}

out <- ode (y=y0,t=t0,func =f1,parms = NULL )
dat_r[,i] <-out[,2]
}
out1 <- data.frame(dat_r)
out1 <-cbind(out[,1],out1)

#head(out1)

g42 <-ggplot(data=out1, aes(x=out1[,1],y=out1[,2]))
g42 <-g42 + geom_line(data=out1,aes(x=out1[,1],y=out1[,2]),
colour="red",linetype=1,size=0.8)
g42 <-g42 + geom_line(data=out1,aes(x=out1[,1],y=out1[,3]),
colour="blue",linetype=1,size=0.8)
g42 <-g42 + geom_line(data=out1,aes(x=out1[,1],y=out1[,4]),
colour="black",linetype=1,size=0.8)
g42 <-g42 + geom_text(aes(label=rtext[1], 0.5*max(out1[,1]),max(out1[,4])*0.8),
size=7,colour="red")
g42 <-g42 + geom_text(aes(label=rtext[2], 0.5*max(out1[,1]),max(out1[,4])*0.7),
size=7,colour="blue")
g42 <-g42 + geom_text(aes(label=rtext[3], 0.5*max(out1[,1]),max(out1[,4])*0.6),
size=7,colour="black")
g42 <-g42 + geom_text(aes(label=y0text, 0.5*max(out1[,1]),max(out1[,4])*0.9),
size=7,colour="green")
g42 <-g42 +labs(x = "t",y="Y",title = "ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ")
#g42 <-g42 + theme_stata() + scale_colour_stata()

g42

```

* * *

2.Скопируйте  к себе на компьютер, данный скрипт.

Запуститу его в $guiR$ или $Rstudio$. Используйте этот скрипт для исследования поведения модели Мальтуса. Поэкспериментируйте с программой, наблюдая отклики модели при различных планах изменения параметра r.На основе этого скрипта, самостоятельно разработайте свой скрипт, который позволяет автоматически получить произвольное число значений $r$ и строить соотвествующее семейство графиков $X(t)$ при разных значениях $X_0$.

**скрипт 3**


```{r,fig.height=5,fig.width=6}
#==ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ-произвольное число значений r ====

library(deSolve) 


yy0=55 #начальная численность популяции
ri <-seq (0.1 ,2.5 ,0.2)#диапозон изменения параметра r
n_r=length(ri)# чило значений параметра r

rtext =vector('character',length(ri))
for(l in 1:n_r){
rtext[l]= 'r ='
}
p_y0<- seq (5 ,60 ,5)#диапозон изменения параметра y0
j = which(p_y0 == yy0)
y0text <-paste('y0 =',as.character(p_y0[j]))

p_t0 <<-0
p_dt <<-0.1
p_tn <<-1.0#интервал исследования
t0<- seq (p_t0 ,p_tn ,p_dt)
y0 <-p_y0[j]
n_j <- vector('numeric',length(t0))
dat_r <-matrix(n_j, nrow = length(t0), ncol =length(ri))
#dat_r 

for(i in 1:n_r){
r <-ri[i]
rt <-as.character(r)
rtext[i] <-paste(rtext[i],rt)

f1<- function (t,y, parms ){
dY.dt <-y [1]*r
return (list (dY.dt ))}

out <- ode (y=y0,t=t0,func =f1,parms = NULL )
dat_r[,i] <-out[,2]
}
#dat_r
out1 <- data.frame(dat_r)
out1 <-cbind(out[,1],out1)
#head(out1)

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r)



#x11()
plot(out1[,1],out1[,n_r+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=3,font = 1,font.axis=2,col = rcol[n_r],
main = "ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ")


for(i in 2:n_r){
#rcol=50+10*i
lines(out1[,1],out1[,i],col=rcol[i-1],lwd=3)
}

legend ("topleft", lwd = 2, bty='o',title='параметр r',
c(y0text,rtext), col = c("white",rcol),cex = 0.9)

```



3.Создайте $Rmd$ документ, который содержит:

   - Наименование доумента и фамилию студента;
   - Постановку задачи (историческая справка, математическая
     модель   Мальтуса и т. д.);
   - Всроенные скрипты и результаты их работы.
   - Интерпретацию полученных результатов, ответы на 
     контрольные вопросы и выводы.

Контрольные вопросы:

- чем отличаются полученные кинетические кривые?
- как получить эти кривые в логарифмическом масштабе?

4.Сформируйте на основе созданного Вами  $Rmd$  файла  $HTML$  файл


* * *

### Лабораторная работа.2 " Непрерывная модель логистического роста".

Модель логистического роста была предложена Ферхюльстом в 1838 г. описания развития популяции в условиях ограниченных ресурсов питания. В основу модели положено уравнение:

$dX/dt = r\cdot X -b\cdot X^{2}$ , которое приводится к виду: $dX/dt = r\cdot X\cdot(1-X/K)$

Член - $b\cdot X^{2}$ , пропорциональный количеству встреч между особями учитывает  "самоотравление" популяции, объяснимое многими причинами (конкуренция за ресурсы питания, выделение в среду вредного метаболита и др.).Коэффициент $b$ называется коэффициентом внутривидовой конкуренции. Величина $К = г/b$ соответствует устойчивому стационарному состоянию с максимально возможной в данных условиях численностью популяции и называется "емкостью среды". Параметр $r$ называется скоростью роста и характеризует способность популяции к увеличению численности. Решением дифференциального уравнения является функция:

$$X(t)=X_0/(X_0+(K-X_0)\cdot exp(-r\cdot t))$$

$X_0$ - начальная численность популяции.
Характер логистической кривой зависит oт величин параметров $r$ и $К$ и начальной численности $X_0$

* * *

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

**скрипт 4**

```{r,fig.height=4,fig.width=5}

#========логистическая модель =================

library(deSolve) # решение дифф. уравнений с начальными условиями
library(ggplot2)
#library(ggthemes)

#-----------ПАРАМЕТРЫ МОДЕЛИ--------------------

p_r <<-0.2  # скорость роста популяции
K <<-1000   # емкость среды
p_y0 <<-100 # начальная численность популяции

p_t0 <<-0   # начальная точка интервала исследования
p_dt <<-0.1 # шаг 
p_tn <<-200 # конечная точка интервала иссдедования

#------------------------------------------------

f1<- function (t,y, parms ){
dY.dt <-y [1]*r*(1-y[1]/K)
return (list (dY.dt ))}

t0<- seq (p_t0 ,p_tn ,p_dt)
r <-p_r
y0 <-p_y0

out <- ode (y=y0,t=t0,func =f1,parms = NULL )
out <- data.frame(out)
head(out)

g41 <-ggplot(data=out, aes(x=out[,1],y=out[,2]))
g41 <-g41 + geom_line(data=out,aes(x=out[,1],y=out[,2]),
colour="red",linetype=1,size=0.8)
g41 <-g41 +labs(x = "t",y="Y",title = "ЛОГИСТИЧЕСКАЯ МОДЕЛЬ")
#g41 <-g41 + theme_stata() + scale_colour_stata()


print(g41)

```
* * *

#### Порядок выполнения работы

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

1.Скопируйте  к себе на компьютер, выше привенный скрипт.Задайте  значения $r$ – скорости роста популяции ( $0 < r < 0,5$), начальное значение численности $X_0$ ($10 < X_ < 100$) и емкость среды $K$ ($100 < K < 10000$). Задайте интервал моделирования в диапазоне (**20-200**). Запустите модель на исполнение в  $guiR$ или $Rstudio$. Поэкспериментируйте с моделью, оцените полученные результаты.

2.Самостоятельно разработайте свой скрипт, который позволяет получать два семейства графиков.Семейство графиков зависимости $Xср(r)$ при разных значениях $K$ и семейство графиков $X(t)$ при разных значениях $r$ и заданном значении $K$.Поэкспериментируйте со скриптом, наблюдая отклики модели при различных планах изменения параметров.

**скрипт 5**

```{r,fig.height=5,fig.width=6}

#======== дискретная  логистическая модель =================

library(deSolve) # решение дифф. уравнений с начальными условиями

p_t0 <<-0.0   #начальная точка интервала исследования
p_dt <<-0.2   #шаг приращения t
p_tn <<-20 #конечная точка интервала исследования 

K2 <- seq (50 ,300 ,50) #диапозон емкостей среды


K <<-K2[2]   #емкость среды
Ktext <-paste('K =',as.character(K))
p_y0 <<-50 #начальная численность популяции
y0text <-paste('y0 =',as.character(p_y0))


t0<- seq (p_t0 ,p_tn ,p_dt)
#print(t0)

ri1 <-seq (0.1 ,1.2 ,0.1)#  изменения параметра r
n_r1=length(ri1)# чило значений параметра r
n_j <- vector('numeric',length(t0))
dat_r1 <-matrix(n_j, nrow = length(t0), ncol =length(ri1))
dat_K <-array(0,c(length(t0),length(ri1),length(K2)))


rtext =vector('character',length(ri1))
for(l in 1:n_r1){
rtext[l]= 'r ='
}

for(i in 1:n_r1){
r <-ri1[i]
rt <-as.character(r)
rtext[i] <-paste(rtext[i],rt)
}

nk=length(K2)-1
nkk=nk+1

K2text=vector('character',nkk)
for(l in 1:nkk){
K2text[l]= 'K ='
}

for(i in 1:nkk){
KK2 <-K2[i]
K2t <-as.character(KK2)
K2text[i] <-paste(K2text[i],K2t)
}

y0 <-p_y0

for(j in 1:n_r1){
r <-ri1[j]
for(i in 1:length(t0)){
f1<- function (t,y, parms ){
dY.dt <-y [1]*r*(1-y[1]/K)
return (list (dY.dt ))}

out <- ode (y=y0,t=t0,func =f1,parms = NULL )
}
dat_r1[,j]=out[,2]
}
dat_r1 <-cbind(out[,1],dat_r1)
#head(dat_r1) 

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r1)

#x11()
plot(dat_r1[,1],dat_r1[,n_r1+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=3,font = 1,font.axis=2,col = rcol[n_r1],cex.axis=1.2,
cex.lab=1.5,cex.main=1.3,main = "Непрерывная логистическая модель")

for(i in 2:n_r1){
lines(dat_r1[,1],dat_r1[,i],col=rcol[i-1],lwd=3)
}
legend ("bottomright", lwd = 2, bty='n',title='Параметры: у0,К,r',
c(y0text,Ktext,rtext), col = c("white","white",rcol),cex = 1.0)

#===========================================================================
for(l in 1:length(K2)){
K=K2[l]
y0 <-p_y0
for(j in 1:n_r1){
r <-ri1[j]
for(i in 1:length(t0)){
f1<- function (t,y, parms ){
dY.dt <-y [1]*r*(1-y[1]/K)
return (list (dY.dt ))}
out <- ode (y=y0,t=t0,func =f1,parms = NULL )
}
dat_K[,j,l]=out[,2]
}
}

#dat_K

ty=vector('numeric',length(ri1))

for(l in 2:length(K2)){
mean_datr=matrix(dat_K[,,l],nrow =length(t0) ,ncol =length(ri1))
mean_datr=colMeans(mean_datr)
mean_datr <-rbind(ty,mean_datr)
ty=mean_datr
}

mean_datr <-mean_datr[-1,]

nk=length(K2)-1
pal <- colorRampPalette(c('red','green'))
rcol <-pal(nk)

mas1 =max(mean_datr)
mas2=min(mean_datr)  
rmas=c(min(ri1),max(ri1))
mas=c(mas1,mas2)

#x11()
plot(rmas,mas,xlab = "r", ylab = "Yср(t)",
type = "l", pch = 19,lwd=3,font = 1,font.axis=2,col ="white" ,cex.axis=1.2,
cex.lab=1.5,cex.main=1.3,main = "Непрерывная логистическая модель")

for(m in 1:nk){
lines(ri1,mean_datr[m,],col=rcol[m],lwd=3)
}

legend ("topleft", lwd = 2, bty='n',title='Параметр: К',
c(K2text[-1]), col = c(rcol),cex = 1.0)



```


3.Создайте $Rmd$ документ, который содержит:

   - Наименование доумента и фамилию студента;
   
   - Постановку задачи (историческая справка, математическая
     модель Ферхюльста и т. д.);
     
   - Всроенные скрипты и результаты их работы.
   
   - Интерпретацию полученных результатов, ответы 
     на контрольные вопросы и выводы.

Контрольные вопросы:

- Какие характеристики логистической кривой определяются 
  константой скорости, емкостью среды и начальной численностью?
- Какие начальные условия приводят к решению с точкой перегиба?
- Как меняется форма кривой, если $X_0< 0.5\cdot K$; 
  $0.5\cdot K< X_0 < K$; $X_0>K$;  $X_0 = K$

4.Сформируйте на основе созданного Вами  $Rmd$  файла  $HTML$  файл


* * *


### Лабораторная работа 3 "Дискретная модель логистического роста".

Дискретные модели применяются для описания развития популяций, численность которых в момент времени t зависит численности в k предшествующих моментов времени: 

$$X(t)= F(X(t-1),X(t-2),...,X(t-k))$$

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

$$dX/dt = r\cdot X\cdot(1-X/K)$$

Заменив   $dX/dt$  на   $\Delta X/\Delta t$  получим  $\Delta X =X(t+1)-X(t)$ и  $\Delta t =1$ 

$$X(t+1=X(t)\cdot (1+r\cdot (1-X(t)/K))$$

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

$$X(t+1=X(t)\cdot exp(r\cdot (1-X(t)/K))$$

Это уравнение можно считать разностным аналогом логистического уравнения.
При различных соотношениях параметров $r$ и $К$, пользуясь этой моделью, можно получать различные режимы динамики численности популяции:

-	$0 < r < 1$		- монотонное приближение численности к стационарной;
-	$1 < r < 2$	- затухающие колебания;
-	$2 < г < 3$		- стационарные колебания;
-	$r > 3,1$		- нерегулярное поведение ( хаос ).

###### *Порядок выполнения*

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

1.Разработайте скрипт исследования дискретной модели логистического роста с построением соответствующих графиков. Скрипт для этой модели может иметь следующий вид:

**скрипт 6**

```{r,fig.height=4,fig.width=5}

#======== дискретная  логистическая модель =================

library(deSolve) # решение дифф. уравнений с начальными условиями

p_t0 <<-1   #начальная точка интервала исследования
p_dt <<-1   #шаг приращения t
p_tn <<-50 #конечная точка интервала исследования 

p_r <<-2.8

K <<-200   #емкость среды
Ktext <-paste('K =',as.character(K))
p_y0 <<-50 #начальная численность популяции
y0text <-paste('y0 =',as.character(p_y0))


t<- seq (p_t0 ,p_tn ,p_dt)
#print(t)
n=as.integer(p_tn/p_dt)
x=numeric(n)

ri1 <-seq (0.2 ,1.0 ,0.2)# первый диапозон изменения параметра r
n_r1=length(ri1)# чило значений параметра r
n_j <- vector('numeric',length(t))
dat_r1 <-matrix(n_j, nrow = length(t), ncol =length(ri1))

for(j in 1:n_r1){
x[1]=p_y0
r <-ri1[j]
for(i in 1:n){
x[i+1]=x[i]*exp(r*(1-x[i]/K))
x
}
dat_r1[,j]=x[1:n]
}
dat_r1 <-cbind(t,dat_r1)
#head(dat_r1) 

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r1)

#x11()
plot(dat_r1[,1],dat_r1[,n_r1+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=3,font = 1,font.axis=2,col = rcol[n_r1],
main = "Дискретная модель	  0 < r < 1")

for(i in 2:n_r1){
lines(dat_r1[,1],dat_r1[,i],col=rcol[i-1],lwd=3)
}
legend ("bottomright", lwd = 2, bty='n',title='Параметры: у0,К',
c(y0text,Ktext), col = c("white","white"),cex = 1.2)

#=================================================================
ri2 <-seq (1.2 ,2.0 ,0.2)# второй диапозон изменения параметра r
n_r2=length(ri2)# чило значений параметра r
n_j <- vector('numeric',length(t))
dat_r2 <-matrix(n_j, nrow = length(t), ncol =length(ri2))

for(j in 1:n_r2){
x[1]=p_y0
r <-ri2[j]
for(i in 1:n){
x[i+1]=x[i]*exp(r*(1-x[i]/K))
x
}
dat_r2[,j]=x[1:n]
}
dat_r2 <-cbind(t,dat_r2)
#head(dat_r2) 

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r2)

#x11()
plot(dat_r2[,1],dat_r2[,n_r2+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=2,font = 1,font.axis=2,col = rcol[n_r1],
main = "Дискретная модель	  1 < r < 2")

for(i in 2:n_r2){
lines(dat_r2[,1],dat_r2[,i],col=rcol[i-1],lwd=2)
}
legend ("bottomright", lwd = 2, bty='n',title='Параметры: у0,К',
c(y0text,Ktext), col = c("white","white"),cex = 1.2)

#========================================================================

ri3 <-seq (2.2 ,3.0 ,0.2)# третий диапозон изменения параметра r
n_r3=length(ri3)# чило значений параметра r
n_j <- vector('numeric',length(t))
dat_r3 <-matrix(n_j, nrow = length(t), ncol =length(ri3))

for(j in 1:n_r3){
x[1]=p_y0
r <-ri3[j]
for(i in 1:n){
x[i+1]=x[i]*exp(r*(1-x[i]/K))
x
}
dat_r3[,j]=x[1:n]
}
dat_r3 <-cbind(t,dat_r3)
#head(dat_r3) 

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r3)

#x11()
plot(dat_r3[,1],dat_r3[,n_r3+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=2,font = 1,font.axis=2,col = rcol[n_r1],
main = "Дискретная модель	  2 < r < 3")

for(i in 2:n_r3){
lines(dat_r3[,1],dat_r3[,i],col=rcol[i-1],lwd=2)
}
legend ("top", lwd = 2, bty='n',title='Параметры: у0,К',
c(y0text,Ktext), col = c("white","white"),cex = 1.2)

#========================================================================

ri4 <-seq (3.5 ,8.0 ,1.5)# четвертый диапозон изменения параметра r
n_r4=length(ri4)# чило значений параметра r
n_j <- vector('numeric',length(t))
dat_r4 <-matrix(n_j, nrow = length(t), ncol =length(ri4))

for(j in 1:n_r4){
x[1]=p_y0
r <-ri4[j]
for(i in 1:n){
x[i+1]=x[i]*exp(r*(1-x[i]/K))
x
}
dat_r4[,j]=x[1:n]
}
dat_r4 <-cbind(t,dat_r4)
#head(dat_r4) 

pal <- colorRampPalette(c('red','green'))
rcol <-pal(n_r4)

#x11()
plot(dat_r4[,1],dat_r4[,n_r4+1],xlab = "t", ylab = "Y(t)",
type = "l", pch = 19,lwd=2,font = 1,font.axis=2,col = rcol[n_r1],
main = "Дискретная модель	  3 < r < 8")

for(i in 2:n_r4){
lines(dat_r4[,1],dat_r4[,i],col=rcol[i-1],lwd=2)
}
legend ("top", lwd = 2, bty='n',title='Параметры: у0,К',
c(y0text,Ktext), col = c("white","white"),cex = 1.2)




```



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

2.Создайте $Rmd$ документ, который содержит:

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

Контрольные вопросы:

- При каких значениях $r$ дискретная и непрерывная 
  логистическая модели дают одинаковые решения? 
- При каких значениях $r$ получаются устойчивые решения
  и в чем их отличие друг от друга?
- При каких значениях $r$ модель может описывать хаотические
  вспышки численности насекомых?

4.Сформируйте на основе созданного Вами  $Rmd$  файла  $HTML$  файл


* * *
### Лабораторная работа 4  Исследование динамики двух видовых популяций

**Основная цель** работы: исследование взаимоотношений «хищник-жертва» на фазовой плоскости и во времени методом математического моделирования.

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

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

$$dx_1/dt = a_1\cdot x_1 + b_{12}\cdot x_1\cdot x_2 - c_1\cdot x_1^{2}$$
$$dx_2/dt = a_2\cdot x_2 + b_{21}\cdot x_1\cdot x_2 - c_2\cdot x_2^{2}$$

Здесь параметры $a_i$  - константы собственной скорости роста видов, $c_i$  - константы самоограничения численности (внутривидовой конкуренции),$b_{ij}$   - константы взаимодействия видов, $(i,j=1,2)$. Соответствие знаков этих последних коэффициентов различным типам взаимодействий приведено в таблице:

![Таблица типов взаимодействий популяций](tab01.png)

Исследование свойств моделей подобного типа приводит к некоторым важным выводам относительно исхода взаимодействия видов.
Уравнения конкуренции $( b_{12}>0 , b_{21} <0)$ предсказывают выживание одного из двух видов, в случае если собственная скорость роста другого вида меньше не-которой критической величины. Оба вида могут сосуществовать, если произведение коэффициентов межпопуляционного взаимодействия меньше произведе-ния коэффициентов внутри популяционного взаимодействия: $b_{12}b_{21} < c_2c_1$ .

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

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

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

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

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

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

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

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

---

##### Классическая модель Лотки и Вольтера.

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

Рассмотрим модель взаимодействия хищников и их добычи, когда между особями одного вида нет соперничества. Пусть $x_1$  и  $x_2$  число жертв и хищников соответственно. Предположим, что относительный прирост жертв $x_1^{'}/x_1$ = $\alpha-\beta x_2$, $\alpha>0$, $\beta>0$, где $\alpha$ — скорость размножения жертв в отсутствие хищников, $\beta x_2$ - потери от хищников.
Развитие популяции хищников зависит от количества пищи (жертв), при отсутствии пищи ( $x_1=0$ ) относительная скорость изменения популяции хищников $x_2^{'}/x_2 =-\gamma$, $\gamma >0$ , наличие пищи компенсирует убывание, и при $x_1>0$ имеем:
$$x_2^{'}/x_2 = (-\gamma + \delta x_1), \delta >0$$

Систему Вольтера - Лотка можно представить в следующем виде:

$$dx_1/dt =\alpha x_1-\beta x_2x_1$$
$$dx_2/dt =-\gamma x_2 + \delta x_1x_1$$

$\alpha,\beta , \gamma , \delta > 0$

Скрипт реализующий решение этой системы  дифференциальных уравнений  может иметь следующий вид:

**скрипт 7**

```{r,fig.height=5,fig.width=6}
#----классическая модель Лотки и Вольтера----
library (deSolve)

#--Параметры модели--------------------------

alfa =4.0
beta=2.5
gamma=2.0
delta=1.0

p_tn =20
p_dt =0.1
p_t0=0

x10=3
x20=1


#Решение системы диффуров классической модели---

f4 <- function (t, y, parms ){
with(as.list(y),{
dX.dt <- alfa*X - beta*Y*X
dY.dt <- -gamma*Y + delta*X*Y
list(c(dX.dt,dY.dt))})}

t0<- seq (p_t0 ,p_tn ,p_dt)
y0 <- c(X=x10,Y=x20)

out <- ode (y=y0,t=t0,func =f4,parms = NULL )

#x11()
plot (out [,"X"], out [,"Y"],col ="blue",xlab ="x1",ylab ="x2",
type = "l", pch = 19,lwd=3,cex.axis=1.2,cex.lab=1.5,cex.main=1.8,
main ="ФАЗОВЫЙ ПОРТРЕТ x1,x2")

#x11()
plot (out [,1], out [,2],col ="cyan4",xlab ="t",ylab ="x1,x2",
type = "l", pch = 19,lwd=3,cex.axis=1.2,cex.lab=1.5,cex.main=1.0,
main ="Численность популяции жертв x1 и хищников x2")
lines(out [,1],out [,3],col="brown4",lwd=3)
legend ("top", lwd = 2, bty='n',title='',
c("жертвы","хищники"), col = c("cyan4","brown4"),cex = 1.8)

```

Решение получено при  следующих значениях параметров: $\alpha=4 ,\beta = 2.5 ,\gamma =2 ,\delta =1$ и с начальным условием $x_1(0)=3$, $x_2(0)=1$.

Видно, что процесс имеет колебательный характер. При заданном начальном соотношении числа особей обоих видов 3 : 1, обе популяции сначала растут. Когда число хищников достигает величины $\beta$, популяция жертв не успевает восстанавливаться и число жертв начинает убывать. Уменьшение количества пищи через некоторое время начинает сказываться на популяции хищников и когда число жертв достигает величины $x_1=\gamma/\delta$ (в этой точке $x_2^{'}=0$), число хищников тоже начинает сокращаться вместе с сокращением числа жертв. Сокращение популяций происходит до тех пор, пока число хищников не достигнет величины $x_2= \alpha/\beta$ (в этой точке $x_1^{'}=0$).С этого момента начинает расти популяция жертв, через некоторое время пищи становится достаточно, чтобы обеспечить прирост хищников, обе популяции растут, и процесс повторяется снова и снова. На графике четко виден периодический характер процесса. Периодичность процесса явственно видна на фазовой плоскости — фазовая кривая ($x_1(t), x_2(t)$) — замкнутая линия. Самая левая точка, этой кривой, - это точка, в которой число жертв достигает наименьшего значения. Самая правая точка, - точка пика популяции жертв. Между этими точками количество хищников сначала убывает, до нижней точки фазовой кривой, где достигает наименьшего значения, а затем растет до верхней точки фазовой кривой. Если в начальный момент система находилась в стационарной точке, то решения $x_1(t), x_2(t)$ не будут изменяться во времени, останутся постоянными. Всякое же другое начальное состояние приводит к периодическому колебанию решений. Неэллиптичность формы траектории, охватывающей центр, отражает негармонический характер колебаний.


---

#### Уравнения Вольтерра-Лотка с логистической поправкой.

Рассмотрим модель конкурирующих видов с "логистической поправкой":

$$dx_1/dt =\alpha x_1-\beta x_2x_1 - \mu x_1^{2}$$
$$dx_2/dt =-\gamma x_2 + \delta x_1x_1 -\mu x_2^{2}$$

Внесем в скрипт необходимые изменения и получим следующие результаты ($\mu =0.1$)

**скрипт 8**

```{r,fig.height=5,fig.width=6}
#----модель Лотки и Вольтера с поправкой------
library (deSolve)

#--Параметры модели--------------------------

alfa =4.0
beta=2.5
gamma=2.0
delta=1.0
mu = 0.1

p_tn =20
p_dt =0.1
p_t0=0

x10=3
x20=1

#Решение системы диффуров модели с логистической поправкой

f5 <- function (t, y, parms ){
with(as.list(y),{
dX.dt <- alfa*X - beta*Y*X - mu*X^2
dY.dt <- -gamma*Y + delta*X*Y - mu*Y^2
 
list(c(dX.dt,dY.dt))})}

t0<- seq (p_t0 ,p_tn ,p_dt)
y0 <- c(X=x10,Y=x20)

out <- ode (y=y0,t=t0,func =f5,parms = NULL )

#x11()
plot (out [,"X"], out [,"Y"],col ="blue",xlab ="x1",ylab ="x2",
type = "l", pch = 19,lwd=3,cex.axis=1.2,cex.lab=1.5,cex.main=1.0,
main ="Фазовый портрет с логистической поправкой")

#x11()
plot (out [,1], out [,2],col ="cyan4",xlab ="t",ylab ="x1,x2",
type = "l", pch = 19,lwd=3,cex.axis=1.2,cex.lab=1.5,cex.main=1.2,
main ="Численность популяции жертв x1 и хищников x2")
lines(out [,1],out [,3],col="brown4",lwd=3)
legend ("top", lwd = 2, bty='n',title='',
c("жертвы","хищники"), col = c("cyan4","brown4"),cex = 1.8)

```

В этом случае поведение решений в окрестности стационарной точки меняется в зависимости от величины и знака параметра $\mu$ . Рассмотрим фазовый портрет системы Вольтерра—Лотка для     $\mu  =0.1, \alpha =4, \beta =2.5, \gamma =2, \delta =1$ и графики ее решения с начальным условием $x_1(0)=3, x_2(0)=1$.

Видно, что в этом случае стационарная точка превращается в устойчивый фокус, а решения — в затухающие колебания. При любом начальном условии состояние системы через некоторое время становится близким к стационарному и стремится к нему при $t\longrightarrow \infty$ . 
В случае при отрицательном значении параметра $\mu$, стационарная точка является неустойчивым фокусом и амплитуда колебаний численности видов растет. В этом случае как бы близко ни было начальное состояние к стационарному, с течением времени состояние системы будет сильно отличаться от стационарного.
Рассмотренные модели могут описывать поведение конкурирующих фирм, рост народонаселения, численность воюющих армий, изменение эколо-гической обстановки, развитие науки и пр.

---



#####  Порядок выполнения задания.

**Цель задания**: Разработка (использование) скриптов для решения соотвествующих систем дифференциальных уравнений и построение на этой основе графиков динамики переменных $x_1$  (жертв) и $x_2$ (хищников) во времени и фазовом пространстве.

1. Задайте в первом скрипте интервал моделирования в диапазоне $(0-20)$, значения вектора коэффициентов правых частей классической системы  $(0< \alpha, \beta, \gamma, \delta, <10)$,начальные значения численности $x_1 ,x_2$  $(0 < x_1 < 5, 0 < x_2< 5)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.

2. Задайте во втором скрипте дополнительно параметр $\mu$ в диапазоне $(0.1-1.0)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.

3. Используя пакет $shiny$, самостоятельно создайте интерактивную модель Вольтера Лотки с логистической поправкой с web инфейсом.

4. Создайте Rmd документ, который содержит:

  - Наименование доумента и фамилию студента;
   
  - Постановку задачи (историческая справка, математическая 
    модель Вольтера, основные принципы взаимодействия 
    популяций и т. д.);
  - Всроенные скрипты и shiny скрипт и результаты их работы.
   
  - Интерпретацию полученных результатов, ответы на контрольные 
    вопросы и выводы.

Контрольные вопросы:

-  При каких значениях параметра   кривые системы (2) и (3) подобны?

-  Как значения параметра $\mu$  для системы с поправкой влияют 
   на форму кривых?
-  Как начальные значения $x_1$ и $x_2$ влияют на форму кривых?

-  Дайте определение константам $\alpha,\beta , \gamma , \delta$, 
   как их значения влияют на формы кривых классической системы и
   системы с поправкой.

5.Сформируйте на основе созданного Вами  $Rmd$  файла  $HTML$  файл

---

## Литература
1. Гиляров А.М. Популяционная экология, М.: Из-во МГУ, 1990.
2. Демидович Б.П. Лекции по математической теории устойчивости. М.: Наука, 1967
3. Зарядов И.С. Введение в статистический пакет R. Учебно-методическое пособие, РУДН, 2010.
4. Кабаков Р. R в действии. Анализ и визуализация данных в программе R/пер.с англ. М.: ДМК Пресс, 2016.-588 с.: ил. 
5. Мастицкий С. Э.Визуализация данных с помощью ggplot2.М.: ДМК Пресс, 2017. 
6. Пискунов Н.С. Дифференциальное и интегральное исчисления. Для втузов. Учеб-ник. М.: Наука, 1985.
7. Полуэктов PA., Пых ЮА., Швытов ИА. Динамические модели экологи-ческих систем. Л.: Гидрометеоиздат. 1980.
8. Ризниченко Г.Ю., Рубин А.Ь. Математические модели биологических продукционных процессов: Учебное пособие. М.: Из-во МГУ, 1993.
9. Смит Дж. М, Модели в экологии, М.: Мир, 1976.
10. Федоров В.Д., Гильманов Г. Г. Экология. М.: Из-во МГУ. 1980.
11. Yihui Xie, J.J. Allaire, Garrett Grolemund (2019).R Markdown
The Definitive Guide. Published July 31, 2018 by Chapman and Hall/CRC
304 Pages.[ссылка на книгу](https://bookdown.org/yihui/rmarkdown/)
12. Learn Shiny [ссылка на сайт](https://shiny.rstudio.com/tutorial/)

---

###  Приложение

#### Элементы математические теории устойчивости систем 

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


$$\bar X^{'} = f(t,\bar X)\ldots(1)$$

с начальными условиями $\bar X(t_0)=\bar X_0$ где $\bar X=(X_1,X_2,\ldots,X_n)$

$n$ - мерный вектор; $t \in I = [t_0, + \infty)$ - независимая переменная, по которой производится дифференцирование;

$\bar X^{'}=(dX_1/dt,dX_2/dt,\ldots,dX_n/dt)$ - $n$ -мерный вектор производных по $t$

$\bar f(t,\bar X) = (f_1(t,\bar X),f_2(t,\bar X),\ldots,f_n(t,\bar X))$   $n$-мерная вектор - функция.

Если начальные данные $(t_0 ,\bar X_0)$ изменяются, то изменяется и решение. Тот факт, что решение зависит от начальных данных, обозначается следующим образом:$\bar X(t) =\bar X(t;t_0,\bar X_0)$   . Естественно, что в качестве математической модели пригодна лишь та задача Коши, которая устойчива к малым изменениям начальных данных.
Определим понятие устойчивости, асимптотической устойчивости и неустойчивости в смысле Ляпунова.
Для этого отклонение решения $\bar X(t) =\bar X(t;t_0,\bar X_0)$, вызванное отклонением $\bar \Delta X_0$ начального значения $\bar X_0$, будем записывать следующим образом:
$$\mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid = \mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t;t_0,\bar X_0)\mid $$

**Определение 1.** Решение $\bar X(t) =\bar X(t;t_0,\bar X_0)$ системы (1) называется устойчивым по Ляпунову в положительном направлении (или устойчивым), если оно непрерывно по $\bar X_0$ на интервале $I = [t_0, + \infty)$, т.е. $\forall \epsilon >0$  $\exists \delta >0$ такое, что если для $\lor \bar \Delta X_0$    $\mid \bar \Delta X_0 \mid \leq \delta$ $\Rightarrow$ $\mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid \leq \delta$  $\      \forall t \geq t_0$. Если, кроме того, если отклонение решения $\bar X(t)$  стремится к нулю при $t → + ∞$ для достаточно малых $\bar \Delta X_0$ , т.е. $\exists \Delta > 0$, что $\forall \bar \Delta X_0$   $\ \mid \bar \Delta X_0\mid \leq \Delta \Rightarrow \mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid \longrightarrow 0 ,t → + ∞ $.  То решение $\bar X(t)$ системы (1) называется асимптотически устойчивым в положительном направлении (или асимптотически устойчивым). Аналогично определяются различные типы устойчивости решения в отрицательном направлении.

**Определение 2.**Решение $\bar X(t)=\bar X(t;t_0,\bar X_0)$ системы (1) называется неустойчивым по Ляпунову в положительном направлении (или неустойчивым), если оно не является устойчивым в положительном направлении. Аналогично определяется неустойчивость в отрицательном направлении.



---


**Геометрическая интерпретация**

1.  Геометрически устойчивость по Ляпунову решения $\bar X(t)$ можно интерпретировать следующим образом (рисунок 1): все решения ) $\bar X(t;t_0,\bar X_0 +\bar \Delta X_0)$, близкие в начальный момент $t0$ к решению $\bar X(t)$ (т.е. начинающиеся в пределах $δ$ - трубки ), не выходят за пределы $ε$ - трубки при всех значениях $t ≥ t_0$ .
2. Асимптотическая устойчивость есть устойчивость с дополнительным условием. Геометрически это означает, что любое решение $\bar X_1(t)$ , начинающееся в момент $t_0$ в $Δ$ - трубке, с течением времени неограниченно приближается к решению $\bar X(t)$. Трубка радиуса $Δ$ называется областью притяжения решения $\bar X(t)$.
3. Неустойчивость по Ляпунову геометрически означает, что среди решений, близких в начальный момент $t_0$ к решению $\bar X(t)$ найдется хотя бы одно, которое в некоторый момент $t_1$ (свой для каждого такого решения) выйдет за пределы $ε$ - трубки (рисунок 1)/
 
Исследование устойчивости произвольного решения $\bar X(t)$  системы (1) всегда можно свести к исследованию устойчивости нулевого решения некоторой преобразованной системы.
В дальнейшем будем предполагать, что система (1) имеет нулевое решение, т.е. $\bar f(t,0)= 0 \ ∀ t ≥ t_0$ , и ограничимся исследованием устойчивости нулевого решения. Переформулируем определения различных типов устойчивости для нулевого решения $\bar X(t) ≡ 0$ системы (1).

**Определение 3.** Нулевое решение $\bar X(t) ≡ 0$ системы (1) называется устойчи-вым по Ляпунову в положительном направлении (или устойчивым),  если $∀ ε > 0 \ ∃ δ = δ ( ε ) > 0$ такое, что $∀ \bar X_0$

$\mid \bar \Delta X_0 \mid  ≤ δ \Rightarrow \mid \bar X(t;t_0,\bar X_0)\mid ≤ ε$ для $∀ t ≥ t_0$

Если кроме того, $∃ Δ > 0$ , что для $∀ X_0$
$\mid \bar \Delta X_0 \mid  ≤ δ \Rightarrow \mid \bar X(t;t_0,\bar X_0)\mid  → 0,t → + ∞$, то решение $\bar X(t) ≡ 0$ системы (1) называется асимптотически устойчивым в положительном направлении ( или асимптотически устойчивым ).

**Определение 4.** Нулевое решение $\bar X(t) ≡ 0$ системы (1) называется неустойчивым по Ляпунову в положительном направлении (или неустойчиво), если оно не явля-ется устойчивым в положительном направлении.


---

![Рисунок 1 Виды решений по Ляпунову](ris_01.png)


---

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

$$dX/dt =P(X,Y)\ldots(2)$$
$$dY/dt =Q(X,Y)$$

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

Такой подход допускает наглядное представление поведения переменных на фа-зовой плоскости $(X,Y)$. Совокупность траекторий точек $М(X,Y)$ плоскости $(X,Y)$, изображающих (представляющих) значения $X и Y$ в последовательные моменты времени $t$ соответствуют состояниям системы в процессе изменения согласно уравнениям (2).

Множество фазовых траекторий, построенных из различных начальных условий (задаваемых точкой с координатами $(Xо,Yо)$, образует так называемый фазовый портрет системы и позволяет анализировать характер изменений в системе без знания аналитических выражений решений системы уравнений (2).
Особый интерес представляют стационарные состояния (точки покоя $X_{ст}Y_{ст}$.) системы, в которых производные переменных по времени равны нулю.

$$dX/dt \mid X_{ст}Y_{ст} =0 ;dY/dt \mid X_{ст}Y_{ст} =0$$

Точки покоя системы (2) могут быть устойчивыми и неустойчивыми по Ляпунову. Как известно, исследование устойчивости любого, а значит, и постоянного решения можно свести к исследованию устойчивости нулевого решения. Поэтому далее будем считать, что система (2) имеет нулевое решение $X(t)=Y(t)=0$ , т. е.$P(0,0=Q(0,0)=0$ , и точка покоя совпадает с началом координат фазового пространства. Таким образом, устойчивость нулевого решения системы (2) означает устойчивость начала координат фазового пространства системы (2), и наоборот.
Рассмотрим типы стационарных состоянии и их устойчивости для нормальной системы двух линейных дифференциальных уравнений с постоянными коэффициентами:

$$dX/dt= P_1X+P_2Y \ldots(3)$$
$$dY/dt= P_3X+P_4Y$$

Точка $(0,0)$ является точкой покоя системы (3). Исследуем расположение траек-тории системы (3) в окрестности этой точки. Из теории известно, что решение можно представить в следующем виде:

$$X=\alpha_1 e^{k_t},Y=\alpha_2 e^{k_t}$$

Для определения $k$ получаем характеристическое уравнение:

$$\begin{vmatrix}
  P_1-k& P_2\\
  P_3& P_4-k
\end{vmatrix}=0\ldots(4)$$

Рассмотрим возможные случаи.
$I.$ Корни характеристического уравнения действительные и различные.
1) $k_1 < 0, k_2 < 0.$ Точка покоя асимптотически устойчива (устойчивый узел).
2) $k_1 > 0, k_2 > 0.$ Точка покоя неустойчива (неустойчивый узел).
3) $k_1 > 0, k_2 < 0.$ Точка покоя неустойчива (седло).
4) $k_1 = 0, k_2 > 0.$ Точка покоя неустойчива.
5) $k_1 = 0, k_2 < 0.$ Точка покоя устойчива, но не асимптотически

$II.$ Корни характеристического уравнения комплексные:
$k_1 = p + q_i, k_2 = p - q_i$
1) $p < 0 , q ≠ 0.$ Точка покоя асимптотически устойчива (устойчивый фокус).
2) $p > 0 , q ≠ 0.$ Точка покоя неустойчива (неустойчивый фокус).
3) $p = 0, q ≠ 0.$ Точка покоя устойчива (центр). Асимптотической устойчивости нет.

$III.$ Корни кратные: $k1 = k2$
1) $k_1 = k_2 < 0.$ Точка покоя асимптотически устойчива (устойчивый узел).
2) $k_1 = k_2 > 0.$ Точка покоя неустойчива (неустойчивый узел).
3) $k_1 = k_2 = 0.$ Точка покоя неустойчива. Возможен исключительный случай, когда все точки плоскости являются устойчивыми точками покоя.
Для системы (3) двух линейных уравнений с постоянными действительными коэффициентами характеристическое уравнение (4) приводится к виду $k_2 + a_1 k + a_2 = 0.$

1) Если $a_1 > 0 , a_2 > 0$, то нулевое решение системы (3) асимптотически устойчиво.
2) Если $а_1 > 0 , a_2 = 0$, или $a_1 = 0 , a_2 > 0$ , то нулевое решение устойчиво, но не асимптотически.
3) Во всех остальных случаях нулевое решение неустойчиво; однако при $a_1 = a_2 = 0$ возможен исключительный случай, когда нулевое решение устойчиво, но не асимптотически.
Разновидности особых точек в фазовом пространстве для системы (3) представлены на рисунке 2.

![Рисунок 2 Виды особых точек в фазовом пространстве](ris_02.png)


