assignment 2, Ordinära differentialekvationer, MMG511 / MVE162
Author
Elias Kamyab Orvar
License
Creative Commons CC BY 4.0
Abstract
Göteborg universitet / Chalmers
\documentclass{article}
\usepackage[T1]{fontenc} % Svenska bokstäver
\usepackage[utf8]{inputenc} % Teckenkodning i filen
\usepackage[swedish]{babel} % Svenska rubriker i innehållsförteckning etc
\usepackage{amsfonts} % Innehåller t.ex. fonterna \mathbb
\usepackage{amsmath} % Innehåller mycket användara kommandon för matte
\usepackage{amsthm} % Ger omgivningar för att skapa satser och bevis
\usepackage{graphicx} % för att infoga bilder
% Ofta vill man definiera egna kommandon för sådana saker som man skriver ofta
\usepackage{sidecap}
\usepackage{listings}
\newcommand{\tb}{\textbackslash} % Gör att den late kan skriva bara \tb
\newcommand{\R}{\mathbb{R}} % Ger reella talen
\newcommand{\Z}{\mathbb{Z}} % Ger heltalssymbolen
% Skapa olika typer av satser som man vill använda i sin fil (kräver amsthm)
% Det första är namnet på omgivningen som använder när man skapar den
% och det andra är det som skrivs ut.
\newtheorem{sats}{Sats}
\newtheorem{prop}{Proposition}
\newtheorem{lemma}{Lemma}
% Kommentarer inleds med procenttecken och visas inte
% Här kommer en titel och författare att visa i början av dokumnetet
% och i sidhuvudena (om man vill)
\title{assignment 2, Ordinära differentialekvationer, MMG511 / MVE162}
\author{Elias Kamyab Orvar}
\date{25 maj 2016} % \today ger dagens datum
\begin{document}
%Skapa titelinfo
\maketitle
%Lägga till rubriker för sektioner etc.
\section{Inledning till rapporten}
I denna projekt ska vi undersöka beteende av två stycken Predator-Prey system (där bara två arter är inblandade) baserat på två ekologiska modeller dvs Lokta-Volterra modellen och Holling Tanner respons model. Vidare skall vi studera hur stabilitet, instabilitet eller asymptotiska egenskapen av modeller viariera beroende på vilka värde för parametrarna vi väljer. För stöd till vår anaalytiskt analys använder vi även numeriska och grafiska beräkningar.
\section{Analys}
\subsection{Lokta-volterra model}
Den givna systemet in Lokta-Volterra modellen är: \newline $ F(x_1(t), x_2(t))=
\begin{cases}
\dot{x_1} = (a - bx_2)x_1 \implies \frac{\dot{dx_1}}{dt}= ax_1- bx_1x_2 \ \ \ \ \ \ \ \ \ \ (1) \\
\dot{x_2} = (-c + dx_1)x_2 \implies \frac{\dot{dx_2}}{dt}= -cx_2+dx_2x_1 \ \ \ \ \ (2)
\end{cases}\newline $
På systemet ovan $x_1$ och $x_2$ representerar storleken av de dem två populationerna med avseende på tiden. \newline
Förändringsandel av prey per capita är $ \frac{\dot{x_1}}{x_1}= a- bx_2 $
och förändringsandel per capita av predator $ \frac{\dot{x_2}}{x_2}= -c + dx_1 $.
Parametrarna $ a, b, c$ och $d$ är alltid positiva i systemet samt $x_1(0) \neq 0, x_2(0) \neq 0 $; vi kan ej börja analysen med noll prey och predator.\newline
På (1) dvs. på prey-ekvationen ser vi att $\dot{x_1}$ beror på termerna $ax_1$ och $-bx_1x_2$. I saknad av predator $(x_2)$ kommer $\dot{x_1}$ bero bara på termen $ax_1$ därmed blir derivatan (tillväxten) alltid positiv och population av prey $(x_1)$ växer till oändligheten när tiden går mot oändligheten. I denna model antas att prey har obegränsade tillgång till föda och det finns inga andra faktorer som påverkar systemet utifrån dvs. den enda naturliga fiende/begränsningen som prey har, är predator. \newline
Förhållandevist större/ganska stor population av predator gentemot prey-populationen kan medföra en negativt tillväxt av prey och leder till minskning av prey-population ty den andra termen då dominerar.
På (2) dvs. predator-ekvationen beror tillväxten $ \dot{x_2}$ på termerna $-cx_2$ och $ dx_2x_1 $.
Termen $-cx_2$ representerar naturlig dödsfall eller emigration av predator, alltså ju mindre konstanten $c$ är desto långsammare minskar rovdjurbefolknigen (när tiden ökar) samt stor c påverkar tillväxten (derivatan) av predator negativt. \newline Vi ser i (2) att en positiv tillväxt av predator beror på termen $dx_1x_2$. I brist på föda dvs. små $x_1$ eller väldigt små d (födelsetal av predator), dominerar $-cx_2$ över $dx_1x_2$ och då blir derivatan negativt vilket gör att antalet predator minskar med tiden till en tid då prey-populationen återhämtar sig som senare i sin tur leder till att predator-populationen ökar. I avsaknad av prey $(x_1 = 0)$ dör ut hela predator populationen.
\newline
För att hitta stationära punkter sätter vi $ \begin{pmatrix} \dot{x_1} \\ \dot{x_2}
\end{pmatrix}=\begin{pmatrix}
0 \\ 0
\end{pmatrix} $ dvs. $
\begin{cases}
\dot{x_1} = ax_1- bx_1x_2 = 0 \\
\dot{x_2} = -cx_2+dx_2x_1 = 0
\end{cases}\newline $ Då får vi flera stationära punkter, en av dem är origin
$ \begin{pmatrix} x_1 \\ x_2
\end{pmatrix} = \begin{pmatrix}
0 \\ 0 \end{pmatrix}. \newline $
Systemets variational matris $D_f$ som vi döper till $A(x_1, x_2)$ är: \newline
$A(x_1, x_2)=
\begin{bmatrix}
a-bx_2 & -bx_1 \\ dx_2 & dx_1-c
\end{bmatrix}$
$ A(0, 0)=
\begin{bmatrix}
a & 0 \\ 0 & -c
\end{bmatrix} \implies Det(A) < 0$, $\begin{cases}
trace(A)<0 $ om $c>a \\ trace(A)>0$ om $ c<a \\
trace(A)=0$ om $ c=a
\end{cases} \newline \lambda_1=a$ och $\lambda_2 = -c$, alltså origin är en sadel punkt (instabilt).
Den här stationära punkten är ointressant ty det är omöjligt och orealistiskt att vi startar med noll bytesdjur och rovdjur i naturen. \bigbreak
Sätter vi $\dot{x_2}=0$ får vi kordinatlinjen
$ \begin{pmatrix}
\frac{c}{d} \\ \alpha \end{pmatrix}$ där $ \alpha \in R$ godtyckligt, stoppar vi $\alpha$ i $\dot{x_1}$ får vi null-cline (en linje) som är pararellt med y-axeln. \newline $\begin{pmatrix}
x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix}
\frac{c}{d} \\ \alpha \end{pmatrix} \implies \dot{x_1} = (a-b\alpha)\frac{c}{d}=c(\frac{a-b\alpha}{d})$. se figure 1
\begin{figure}[!ht]
\caption{nuclines där antigen x1' eller x2'= 0}
\centering
\includegraphics[width=1\textwidth]{L-Vnuclines.jpg}
\end{figure}
\bigbreak
Sätter vi $\dot{x_1}=0$ får vi kordinatlinjen
$ \begin{pmatrix}
\alpha \\ \frac{a}{b} \end{pmatrix}$ där $ \alpha \in R$ godtyckligt, om sedan stoppar vi det i $\dot{x_2}$ får vi (med hjälp av matlab) null-cline (en linje) som är pararellt med x-axeln. $\begin{pmatrix}
x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix}
\alpha \\ \frac{a}{b} \end{pmatrix} \implies \dot{x_2} = (d\alpha-c)\frac{a}{b}=a(\frac{d\alpha-c}{b}). se figure 1 \newline$
De båda nullclines skär varandra i första kvadranten (första kvadranten är vår invariant set) på punkten $\begin{pmatrix}
x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix}
\frac{c}{d} \\ \frac{a}{b} \end{pmatrix} $ som i sin tur är en stationär punkt.
Nu ska vi kolla egenskapen om nyligen-hittade stationära punkt med hjälp av linjärisering:
$\newline A(\frac{c}{d}, \frac{a}{b}) = \begin{bmatrix}
0 & -\frac{bc}{d} \\ \frac{ad}{b} & 0
\end{bmatrix} \implies trace(A)=0$ och $Det(A)>0. \newline$
Vi får komplexa egenvärde med $Re\lambda_{1,2}(A)=0$ med $\lambda_{1,2} = \pm i\sqrt[]{ac}$, Eftersom realdelen av egenvärdena är noll så enligt definition 4.2.1 [Sze-Bi-Hsu] är vår stationära punkt ej hyperbolisk därmed kan vi inte använda oss av satserna (om linjärisering) på kapital 4 i boken för att bedömma egenskaperna av förnämnda punkten. Denna punkt är center och för att veta mer om dess stabilitet eller ostabilitet samt asymptotiskt egenskap behöver vi använda en annan metod:\newline
Den givna systemet har en första integral (Lyaponuv funktion) som är given på uppgiften: \bigbreak
$V(x_1,x_2) = g(x_1)h(x_2) = (x_1)^{c}exp(-dx_1)(x_2)^{a}exp(-bx_2)$ där \newline $V: \Omega \subset \mathbb{R}^{2} \to \mathbb{R}, V \in C^{1}, \begin{pmatrix}
\frac{c}{d} \\ \frac{a}{b}
\end{pmatrix} = x^{*} \in \Omega, \Omega$ är en öppen mängd i $ \mathbb{R}^{2} $ med: \newline
$ \begin{cases}
g(x_1) = (x_1)^{c}exp(-dx_1) \\ h(x_2) = (x_2)^{a}exp(-bx_2)
\end{cases} \newline$
$ \dot{V} = \nabla V(x_1,x_2)f(x_1, x_2) = \frac{\partial V}{\partial x_1}\dot{x_1}+\frac{\partial V}{\partial x_2}\dot{x_2} =
\begin{bmatrix}
e^{-dx_1}(cx_1^{c-1}-dx_1^{c})x_2^{a}e^{-bx_2} \\ e^{-bx_2}(ax_2^{a-1}-bx_2^{a})x_1^{c}e^{-dx_1}
\end{bmatrix}
\begin{bmatrix}
(a-bx_2)x_1 \\ (-c+dx_1)x_2
\end{bmatrix} = \newline =
\begin{bmatrix}
e^{-dx_1-bx_2}(cx_1^{c-1}x_2^{a}-dx_1^{c}x_2^{a})(ax_1-bx_2x_1) \\ e^{-dx_1-bx_2}(ax_2^{a-1}x_1^{c}-bx_2^{a}x_1^{c})(-cx_2+dx_1x_2)
\end{bmatrix} = \newline = e^{-dx_1-bx_2}
\begin{bmatrix}
acx_1^{c}x_2^{a}-cbx_1^{c}x_2^{a+1}-adx_1^{c+1}x_2^{a}+dbx_1^{c+1}x_2^{a+1} \\
-acx_1^{c}x_2^{a}+cbx_1^{c}x_2^{a+1}+adx_1^{c+1}x_2^{a}-dbx_1^{c+1}x_2^{a+1}
\end{bmatrix} = 0$ (konstant) \newline
\newline
Alltså $h'(x_2)x_2(-c+dx_1) = g'(x_1)x_1(a-bx_2) $ dvs. $\dot{V}(x) \leq 0 \ \forall x \in \Omega$ och enligt sats 5.2.1 $x^{*}$ är stabilt men inte asymptotiskt stabil ty ($\dot{V} < 0 \ \forall x \in \Omega \backslash \{0\} $) gäller inte här. \newline
$\dot{V}(x_1,x_2) \equiv 0$, varje lösning är en periodisk lösning och vår system (1), (2) har en familj av neutrala stabila periodiska orbits. \newline
Nu vet vi att alla våra lösningar är stabila. Om vi ritar nivå-kurvan av V i en 3-D bild ser vi att maximum-värde av V kommer att ligga på $\begin{pmatrix} \frac{c}{d} \\ \frac{a}{b}
\end{pmatrix} $; allt tyder på en periodisk beteende av lösningar så länge vi inte startar exakt i $x^{*}$ eller på axlarna (då ena populationen är noll). Ett sånt periodiskt beteende ses väl på figure 2 och 3 där på figure 3 plottar vi $x_{1}$ och $x_{2}$ som funktioner av tiden, på figure 2 valde vi parametrarna som $(a, b, c, d)=(1, \ 0.75, \ 0.5, \ 0.25)$
\begin{figure}[!ht]
\caption{persiodiska lösningar, L-V model, oilka startvärde av prey och predator}
\centering
\includegraphics[width=1.2\textwidth,height=0.7\textwidth]{L-Vfirstbild.jpg}
\end{figure}
\begin{figure}[!ht]
\caption{Lokta-Volterra prey-predator model över tiden}
\centering
\includegraphics[width=1\textwidth,height=0.25\textwidth]{L-Vx_t_y_t_.jpg}
\end{figure}
\bigbreak
\subsection{Holling Tanner model}
Givet är två differenteial ekvationer av första ordning för prey respektive predator med alla tänkbara parametrar r, s, K, D, $\omega$, J > 0. \newline
$ F(x_1(t), x_2(t))=
\begin{cases}
\dot{x_1} = r(1-\frac{x_1}{K})x_1-\frac{x_1x_2\omega}{(D+x_1)}\ \ \ \ \ \ (3) \\
\dot{x_2} = s(1-(\frac{Jx_2}{x_1}))x_2 \ \ \ \ \ \ \ \ \ \ \ \ (4)
\end{cases}\newline $ På prey-ekvationen (3) ser vi att termen [$r(1-\frac{x_1}{K})x_1=rx_1-\frac{rx_1^{2}}{K}$] består av två del-termer. Termen beror på r och K, med K som carriying kapacity vilket ger en gräns på resurser såsom föda med mera. Detta gör att antalet prey kan inte växa till oändligehten i frånvaro/minskning av antalet predator. Om $x_1 > K$ då populationen av prey kan upplånas/minskar drastiskt i med att allt mat har redan konsumerats. \newline
parametern r är growth-rate av prey som finns på båda deltermenra. Vid stora K dominerar den första del-termen och prey-populationen har lätt för sig att föröka sig; växa i antal. Vid små K blir det tvärtom. \newline
Nästa term som påverkar tillväxten av prey är $-\frac{x_1x_2\omega}{(D+x_1)}$, här ser vi att om vi inte har någon predator i systemet kommer denna term försvinna och prey population växer till en viss nivå K stycken med positivt tillväxt för att sedan minskar igen. så över tiden kommer antalet prey stabilisera sig kring K. På andra sidan vid stora populationer av predator förhållandevis till prey-populationen kommer denna term dominerar och leder till minskning av prey-populationen \newline
På predator ekvationen (4) har vi termerna $sx_2$ och $-\frac{Jsx_2^{2}}{x_1}$. Här ser vi med blotta ögat på första termen att en positivt tillväxt av predator beror på s (growth rate/kan vara födelsetal) och antalet prey $x_1$ (mat), ju mer prey desto snabbare tillväxt av predator och vice versa. Andra termen påverkar tillväxten negativt och beror på antalet prey i nämnaren dvs. för små population av prey kommer denna term blir stort och därmed dominerar över första termen och vi får en negativ tillväxt och vice versa. \newline
Holling Tanner metoden är mer komplicerad men mer realistiskt i ekologiska miljöer ty det handlar om maktspel mellan prey och predator samt naturresurser att enas död blir den andres bröd. \newline
För att hitta stationära punkterna sätter vi båda $\dot{x_1}=\dot{x_2}=0$. Till skillnad från Lokta-Volterra ser vi att systemet är inte definierat i origin. En Stationär punkt som uppfyller båda kraven är $\begin{pmatrix} K \\ 0 \end{pmatrix}$.\newline
För att kolla punktens egenskaper utför vi linjärisering av systemet med systemets varitonsmatris:
$A(x_1, x_2)=
\begin{bmatrix}
r-\frac{2x_1r}{K}-\frac{x_2D\omega}{(D+x_1)^{2}} & -\frac{x_1\omega}{(D+x_1)} \\ \frac{sJx_2^{2}}{x_1^{2}} & s-\frac{2sJx_2}{x_1}
\end{bmatrix}$ \newline
vilket ger oss
$A(K, 0)=
\begin{bmatrix} -r & -\frac{K\omega}{(D+K)} \\ 0 & s \end{bmatrix} med \begin{cases}
Det(A)<0 \ \ alltid \\ trace(A)<0 \ om \ r>s \\ trace(A) > 0 \ om \ r<s \\ trace(A) = 0 \ om \ r=s
\end{cases} \newline \implies $ alltså en sadel punkt (ostabilt) oavsett vad trace(A) är för något. Denna punkt är inte av intresse att analysera ty vi har inga predator i systemet. \newline
Ett sätt att lättare hitta en eventuellt stationär punkt är att vi hitta våra nullclines och utnyttjar deras skärnings punkt som en ny stationär punkt. \newline
För att hitta nulclines sätter vi $\dot{x_2} = 0$ och håller $x_1$ som konstant då får vi en linje med lutning $\frac{1}{J}$. Sätter vi $\dot{x_1} = 0$ och håller $x_2$ som konstant $\implies$
$0=r(1-\frac{x_1}{K})x_1-\frac{x_1^{2}\omega}{J(D+x_1)} $ utvecklar vi den några steg får vi slutligen en andra gradsekvation i form av $x_1^{2}+x_1(D-K+\frac{K\omega}{Jr})-KD=0$ alltså en kurva (paraboll). \newline
Löser vi den lämpliga $x_1$ som ligger på första kvadranten ur andragrads ekvationen får vi då: \newline
$x_1=\frac{1}{2} [ -(D+K(\frac{K\omega}{Jr}-1))+\sqrt[]{(D+K(\frac{K\omega}{Jr}-1))^{2}+4DK} \ ] = x_1^{*}$
Alltså den nya stationära punkten är skärningspunkten mellan linjen och kurvan dvs. \newline $\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} x_1^{*} \\ \frac{x_1^{*}}{J} \end{pmatrix}$, vidare stoppar vi den i systemetes variational matrisen, får vi då som resultat: \newline
$A(x_1^{*}, \frac{x_1^{*}}{J})=
\begin{bmatrix}
r-\frac{2x_1^{*}r}{K}-\frac{x_1^{*}D\omega}{J(D+x_1^{*})^{2}} & -\frac{x_1^{*}\omega}{(D+x_1^{*})} \\ \frac{s}{J} & -s
\end{bmatrix}$ med följande egenskaper: \newline
$trace(A)=r-s-\frac{2x_1^{*}r}{K}-\frac{x_1^{*}D\omega}{J(D+x_1^{*})^{2}}$ \newline
$Det(A)=-rs+\frac{2x_1^{*}rs}{K}+\frac{x_1^{*}Ds\omega}{J(D+x_1^{*})^{2}}+\frac{sx_1^{*}\omega}{J(D+x_1)^{*}} $, (division med rs) ger $\newline \implies -1 + \frac{2x_1^{*}}{K}+\frac{x_1^{*}D\omega}{rJ(D+x_1^{*})^{2}}+\frac{x_1^{*}\omega}{r(D+x_1^{*})} \implies \newline
\begin{cases}
sadelpunkt \ (instabilt) \ om \ \frac{2x_1^{*}}{K}+\frac{x_1^{*}D\omega}{rJ(D+x_1^{*})^{2}}+\frac{x_1^{*}\omega}{r(D+x_1^{*})}<1 \ \ ty \ det(A)<0\\ icke \ sadelpunkt \ om \ \frac{2x_1^{*}}{K}+\frac{x_1^{*}D\omega}{rJ(D+x_1^{*})^{2}}+\frac{x_1^{*}\omega}{r(D+x_1^{*})}>1
\end{cases} \newline$
Ett sätt att undersöka systemetes beteende är att utgå från trace(A) för att analysera stabiliteten.
trace(A) kan bli både positivt och negativt beroende på val av parametrarna. Enligt den numeriska beräkningen om $\frac{K-D}{2}<x_1^{*}<K$ då kommer vi alltid ha en asymptitiskt stabil stationär punkt ty realdelen av egenvärdena blir negativt (se figure 4)
\begin{figure}[!ht]
\caption{asymptotiskt stabilitet}
\centering
\includegraphics[width=1\textwidth,height=0.4\textwidth]{hollingass.jpg}
\end{figure}
på bild 4 har vi valt parametrarna som $(\omega, r, s, K, D, J)=(1, \ 2, \ 2, \ 10, \ 2, \ 1.5) \implies det(A)=1.0768, tra(A)=-3.0629$ med statiänra punkten $ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 7.3776 \\ 4.9184 \end{pmatrix}$
\newline
Däremot om $x_1^{*}<\frac{K-D}{2}$ då visas att vi inte har asymptotiskt stabilt punkt, (se figure5 med tra(A)=0.0078 och Det(A)=-24.8769)\newline
På bild-5 får vi en ring runtomkring statiänra punkten sådant att inga banor (trajectories) inträder. Denna ring kallas $\omega-$limit set eller periodisk orbit. Trapping region eller positivt invariant set (stor boll eller rektangel där alla trajectories dras innår) ligger på första kvadranten. gränserna till denna boll hittas genom att se vad som händer om vi startar längs $x_1-$axeln; om vi starta på $x_1<K$ så ökar vi fram till K sedan blir det noll (stannar vi), startar vi på $x_1>K$ går vi mot K och sedan stannar vi. Alltså K är gränsen på $x_1-$led. Höjden på rektangen eller bollen kan uppskattas som $\frac{K}{J}$ (lite osäker om denna höjd).
Eftersom tra(A)=0.0078 så vi får både egenvärdena med realdelan större än noll vilket leder till att denna staionärpunkt är repeller, så om vi startar tillräckligt nära punkten kommer alla banor lämnar punkten. se figure-6. allt på bilderna 5-6 är av sådant som sägs på Poincare-Bendixson sats [6.1.1 Hsu].
\begin{figure}[!ht]
\caption{periodisk lösning med w-limit set, innehåller enbart en stationär punkt}
\centering
\includegraphics[width=1.1\textwidth,height=0.6\textwidth]{Hollejasy.jpg}
\end{figure}
\begin{figure}[!ht]
\caption{periodisk lösning, med start nära stionärpunkten, repeller}
\centering
\includegraphics[width=1.1\textwidth,height=0.5\textwidth]{repeller.jpg}
\end{figure}
\bigbreak
\bigbreak
\bigbreak
\bigbreak
\bigbreak
\bigbreak
\bigbreak
\section{reslutat}
Som skillnad mellan modellerna bortsett från små detaljer kan man betona att på H-T modellen behövde vi inte använda oss av något test-funktion och linjärisering av systemet räckte till för ett fullständigt analys. men L-V modellen där vi fick realdelerna av egenvärdena som noll vilket gjorde att vi blev tvungna att utföra testfunktion (first integral) för att ta reda på mer om systemets beteende. \newline
H-T modellen är känsligt för val av parmetrarna och lösningar skiftar mellan asysmtotiskt stabilt och icke-asymptotiskt (periodiska lösningar) beroende på val av parametrar men på L-V modellen har vi dock alltid periodiska lösnigar. \bigbreak
H-T modellen upplevs mer realsitiskt än V-L ty vi använder oss av mer parametrar (faktorer) dessutom tas hänsyn till naturresurser.\newline
På L-V modellen beror de två arterna starkt på varandra och populationer växlar fram och tillbaka mellan dem två arterna medan på H-T modellen verkar att populationerna går mot någonstans av jämnvikt vilket man tycker logiskt ty djur anpassar sig väldigt bra till sin omgivning genom evolutionen.\bigbreak
På båda modellerna gäller gamla forna citaten att "enas död blir den andres bröd " mer eller mindre beroende på modell. ökaning av prey leder till ökning av predator, ökning av predator leder till minskning av prey som i sin tur leder till minskning av predaor som i sin tur i senare skede leder till återhämtning av prey samt i ännu senare skede återhämtning av predator.
\bigbreak
\bigbreak
%Avsluta med en innehållsförteckning
\tableofcontents
\begin{verbatim}
% Första modellen (Lotka-Volterra)
temp=[1 0.75 0.5 0.25];
a=temp(1);
b=temp(2);
c=temp(3);
d=temp(4);
f = @(t,Y) [Y(1)*(a-b*Y(2)); Y(2)*(-c+d*Y(1))];
y1 = linspace(-10,20,20);
y2 = linspace(-10,15,20);
[x,y] = meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(x));
t=0;
for i = 1:numel(x)
Yprim = f(t,[x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x,y,u,v,'black');
figure(gcf)
title('Phase Space Plot (rate plot) start på [10, 20] bytesdjur och [10, 20] rovdjur')
xlabel('Antal bytesdjur i någon enhet')
ylabel('Antal rovdjur i anågon enhet')
axis([0, 10, 0, 8]);
hold on
figure(1) % vi ritar streckorna
nullclinex=[];
nullcliney2 = [];
for k=0:55
nullclinex=[nullclinex ((-c*a)/b+(d*k*a)/b)];
end
for j = 0:55
nullcliney2 = [nullcliney2 ((a*c)/d - (b*c*j)/d)];
end
nullcliney=ones(1,56);
nullcliney=nullcliney.*(a/b);
nullclinex2 = ones(1, 56);
nullclinex2 = nullclinex2.*(c/d);
plot(nullclinex,nullcliney, 'red')
hold on
plot(nullclinex2, nullcliney2, 'Blue');
plot(c/d,a/b, 'k*');
for Prey = 1:0.5:4
for Pred = 1:0.5:4
[ts,X] = ode45(f,linspace(0,20,400),[Prey; Pred]);
plot(X(:,1), X(:,2))
%plot(X(1,1), X(1,2),'bo') % starting point
%plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold on
figure(2)
plot(ts, X);
title('Frequency plot: simple Lotka-Volterra predator/prey');
xlabel('time');
ylabel('population');
legend('prey', 'predators');
%hold on
%figure(3)
%(X(:,1)^.c)*exp(-d*X(:,1)).*(X(:, 2).^a)*exp(-b.*X(:,2));
%% andra modellen (Holling Tanner response model)
clear all; clc;
temp=[1 2 2 10 2 1.5];
w=temp(1);
r=temp(2); % growth rate of prey
s=temp(3); %death rate of predator
K=temp(4); % kapacity
D=temp(5); % half saturation
J=temp(6); % maximal birth of predator
f = @(t,Y) [Y(1)*r*(1-Y(1)/K)-(Y(1)*Y(2)*w)/(D+Y(1)); s*Y(2)*(1-(J*Y(2))/Y(1))];
y1 = linspace(0.1,20,300);
y2 = linspace(0.1,20,300);
[x,y] = meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(x));
t=0;
for i = 1:numel(x)
Yprim = f(t,[x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x,y,u,v,'black');
figure(gcf)
title('Phase Plot start på [1:20] bytesdjur och [1:20] rovdjur')
xlabel('Antal bytesdjur i någon enhet')
ylabel('Antal rovdjur i någon enhet')
axis([0, 20, 0, 20]);
hold on
figure(1)
for Prey = 1:1:20
for Pred = 1:1:20
%(c/d, a/b)=(20, 14) är en stationär punkt
[ts,X] = ode45(f,linspace(0,300,3000),[Prey; Pred]);
plot(X(:,1), X(:,2))
%plot(X(1,1), X(1,2),'bo') % starting point
%plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold on
x1=linspace(-3,20);
plot(x1,(-r/(w*K))*(x1.^2+x1.*(D-K)-D*K),'black');
hold on
plot(x1,x1./J,'black')
%
A=((K*w)/(J*r)+D-K);
statpkt=[(-1/2)*(A-sqrt(A^2+4*D*K));(-1/(2*J))*(A-sqrt(A^2+4*D*K))];
%
Jacobi=[(r-((2*r*statpkt(1))/K)-(statpkt(2)*w*D)/(D+statpkt(1))^2) (statpkt(1)*w)/(D+statpkt(1));
J*s*statpkt(2)^2/statpkt(1)^2 (s-(2*J*s*statpkt(2))/statpkt(1))];
tr=trace(Jacobi)
dt = det(Jacobi)
plot(statpkt(1),statpkt(2),'*')
%%
clear all; clc;
temp=[3 6 3 40 2 0.4];
%temp=[3 3 2 50 15 0.8];
w=temp(1);
r=temp(2); % growth rate of prey
s=temp(3); %death rate of predator
K=temp(4); % kapacity
D=temp(5); % half saturation
J=temp(6); % maximal birth of predator
f = @(t,Y) [Y(1)*r*(1-Y(1)/K)-(Y(1)*Y(2)*w)/(D+Y(1)); s*Y(2)*(1-(J*Y(2))/Y(1))];
y1 = linspace(0.1,50,300);
y2 = linspace(0.1,50,300);
[x,y] = meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(x));
t=0;
for i = 1:numel(x)
Yprim = f(t,[x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x,y,u,v,'black');
figure(gcf)
title('Phase Plot start på [1:5:20] bytesdjur och [1:5:20] rovdjur')
xlabel('Antal bytesdjur i någon enhet')
ylabel('Antal rovdjur i någon enhet')
axis([0, 50, 0, 50]);
hold on
figure(1)
for Prey = 4.7
for Pred = 11.9
%(c/d, a/b)=(20, 14) är en stationär punkt
[ts,X] = ode45(f,linspace(0,30,1500),[Prey; Pred]);
plot(X(:,1), X(:,2))
plot(X(1,1), X(1,2),'bo') % starting point
%plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold on
x1=linspace(-2,100,1000);
plot(x1,(-r/(w*K))*(x1.^2+x1.*(D-K)-D*K),'black');
hold on
plot(x1,x1./J,'black')
%
A=((K*w)/(J*r)+D-K);
statpkt=[(-1/2)*(A-sqrt(A^2+4*D*K));(-1/(2*J))*(A-sqrt(A^2+4*D*K))];
%
Jacobi=[(r-((2*r*statpkt(1))/K)-(statpkt(2)*w*D)/(D+statpkt(1))^2) (statpkt(1)*w)/(D+statpkt(1));
J*s*statpkt(2)^2/statpkt(1)^2 (s-(2*J*s*statpkt(2))/statpkt(1))];
tr=trace(Jacobi)
dt = det(Jacobi)
plot(statpkt(1),statpkt(2),'*')
\end{verbatim}
\end{document}