КОМПЬЮТЕРНЫЙ ПРАКТИКУМ ПО КУРСУ

“Системный анализ и основы компьютерного моделирования экосистем”

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

В пособии приведены методические указания и порядок выполнения работ, обеспечивающих компьютерную поддержку авторского курса А.М. Варюхина “Системный анализ и основы моделирования в экосистем”. В каждой работе сформулированы цели ее выполнения, кратко изложены теоретические основы, на базе которых она построена, приведены контрольные вопросы и рекомендуемая литература. Пособие предназначено для студентов специальности 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

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

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

## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone '"Europe/Saratov"'

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

#========ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ 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 ====

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

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

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)

##   time       X1
## 1  0.0 100.0000
## 2  0.1 101.8145
## 3  0.2 103.6581
## 4  0.3 105.5312
## 5  0.4 107.4340
## 6  0.5 109.3669

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

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

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

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

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)$. Соответствие знаков этих последних коэффициентов различным типам взаимодействий приведено в таблице:

Таблица типов взаимодействий популяций

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

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

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

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

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

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

Рассмотрим модель взаимодействия хищников и их добычи, когда между особями одного вида нет соперничества. Пусть $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

#----классическая модель Лотки и Вольтера----
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

#----модель Лотки и Вольтера с поправкой------
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$ (хищников) во времени и фазовом пространстве.

Задайте в первом скрипте интервал моделирования в диапазоне $(0-20)$, значения вектора коэффициентов правых частей классической системы $(0< \alpha, \beta, \gamma, \delta, <10)$,начальные значения численности $x_1 ,x_2$ $(0 < x_1 < 5, 0 < x_2< 5)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.
Задайте во втором скрипте дополнительно параметр $\mu$ в диапазоне $(0.1-1.0)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.
Используя пакет $shiny$, самостоятельно создайте интерактивную модель Вольтера Лотки с логистической поправкой с web инфейсом.
Создайте Rmd документ, который содержит:

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

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

При каких значениях параметра кривые системы (2) и (3) подобны?
Как значения параметра $\mu$ для системы с поправкой влияют на форму кривых?
Как начальные значения $x_1$ и $x_2$ влияют на форму кривых?
Дайте определение константам $\alpha,\beta , \gamma , \delta$, как их значения влияют на формы кривых классической системы и системы с поправкой.

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

Литература

Гиляров А.М. Популяционная экология, М.: Из-во МГУ, 1990.
Демидович Б.П. Лекции по математической теории устойчивости. М.: Наука, 1967
Зарядов И.С. Введение в статистический пакет R. Учебно-методическое пособие, РУДН, 2010.
Кабаков Р. R в действии. Анализ и визуализация данных в программе R/пер.с англ. М.: ДМК Пресс, 2016.-588 с.: ил.
Мастицкий С. Э.Визуализация данных с помощью ggplot2.М.: ДМК Пресс, 2017.
Пискунов Н.С. Дифференциальное и интегральное исчисления. Для втузов. Учеб-ник. М.: Наука, 1985.
Полуэктов PA., Пых ЮА., Швытов ИА. Динамические модели экологи-ческих систем. Л.: Гидрометеоиздат. 1980.
Ризниченко Г.Ю., Рубин А.Ь. Математические модели биологических продукционных процессов: Учебное пособие. М.: Из-во МГУ, 1993.
Смит Дж. М, Модели в экологии, М.: Мир, 1976.
Федоров В.Д., Гильманов Г. Г. Экология. М.: Из-во МГУ. 1980.
Yihui Xie, J.J. Allaire, Garrett Grolemund (2019).R Markdown The Definitive Guide. Published July 31, 2018 by Chapman and Hall/CRC 304 Pages.ссылка на книгу
Learn Shiny ссылка на сайт

Приложение

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

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

\[\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$ $ X_0X(t;t_0,X_0 +X_0)-X(t) ,t → + ∞ $. То решение $\bar X(t)$ системы (1) называется асимптотически устойчивым в положительном направлении (или асимптотически устойчивым). Аналогично определяются различные типы устойчивости решения в отрицательном направлении.

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

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

Геометрически устойчивость по Ляпунову решения $\bar X(t)$ можно интерпретировать следующим образом (рисунок 1): все решения ) $\bar X(t;t_0,\bar X_0 +\bar \Delta X_0)$, близкие в начальный момент $t0$ к решению $\bar X(t)$ (т.е. начинающиеся в пределах $δ$ - трубки ), не выходят за пределы $ε$ - трубки при всех значениях $t ≥ t_0$ .
Асимптотическая устойчивость есть устойчивость с дополнительным условием. Геометрически это означает, что любое решение $\bar X_1(t)$ , начинающееся в момент $t_0$ в $Δ$ - трубке, с течением времени неограниченно приближается к решению $\bar X(t)$. Трубка радиуса $Δ$ называется областью притяжения решения $\bar X(t)$.
Неустойчивость по Ляпунову геометрически означает, что среди решений, близких в начальный момент $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 Виды решений по Ляпунову

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

Система обыкновенных дифференциальных уравнений называется автономной (или стационарной, или консервативной, или динамической), если независимая переменная не входит явно в систему уравнений. Нормальную автономную систему 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.$

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

Рисунок 2 Виды особых точек в фазовом пространстве

Модели популяций

Varykhin A.

22 11 2020

КОМПЬЮТЕРНЫЙ ПРАКТИКУМ ПО КУРСУ

“Системный анализ и основы компьютерного моделирования экосистем”

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

Саратов 2020

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

“Исследование моделей роста популяций”

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

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

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

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

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

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

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

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

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

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

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

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

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

Литература

Приложение

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

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