Octave: Lotka-Voltera

Př. 1.

a) Namodelujte situaci ryby - rybáři, kdy se vzrůstajícím počtem rybářů klesá počet ryb. Ovšem, jakmile je rybářů příliš a ryb málo, začne klesat počet rybářů a tedy opět vzrůstat počet ryb.

b) Zkuste najít ideální počet rybářů a ryb, který vytvoří rovnováhu (rybáři nemusí přicházet/odcházet a ani ryby se nepřemnoží/nevyhynou).

ad a)Využijeme model Lotka-Voltera

dx/dt = alfa . x(t) - beta . x(t) . y(t) (rovnice pro kořist)

alfa - exponenciální nárůst kořisti (uvažujeme neomezené krmivo, místo k žití atd.)

beta - úbytek vlivem predátora

dy/dt = delta . x(t) . y(t) - lambda . y(t) (rovnice pro predátora)

Uvažujme natalitu (porodnost - tj. přirozený přírůstek) ryb 1,5 a úmrtnost (mortalitu) 0,5.

U rybářů srovnatelně míra příchozích - 0,4 a míra odcházejících 0,8. Pak:

alfa=1.5;
beta=0.5;
gama=0.8;
delta=0.4;

Uvažujme na začátku 8 ryb a 8 rybářů v časovém intervalu 10 časových jednotek (dní):

x0=[8,8]; tn=10;

x označím jako x(1) a y jako x(2) a dosadím lotka-volterrovy rovnice do anonymní fce:

f=@(t,x)[alfa*x(1)-beta*x(1)*x(2); delta*x(1)*x(2)-gama*x(2)];

Vyřeším rovnice pomocí ODE45. To jak dlouho rovnice trvá vyřešit, zjistím pomocí příkazů tic a toc (měří se doba vykonání kódu mezi nimi, která se pak automaticky vypíše):

tic, [t,x] = ode45(f,[0,tn],x0), toc;

Nyní si zobrazím graf:

plot(x(:,1),x(:,2)), xlabel ryby, ylabel rybáři

Celý kód:

alfa=1.5;
beta=0.5;
gama=0.8;
delta=0.4;
x0=[8,8]; tn=10;
f=@(t,x)[alfa*x(1)-beta*x(1)*x(2); delta*x(1)*x(2)-gama*x(2)];
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2)), xlabel ryby, ylabel rybáři

Graf:

image.png

Pokud bych chtěl vidět, jak se postupně mění počet ryb/rybářů, místo příkzu plot využiji comet:

comet(x(:,1),x(:,2));

Pozn.: V ukázce si můžeme všimnout, že počet ryb je v určitém momentu blízko k nule (graf se je daleko více blíže k y-ové nežli k x-ové ose). V praxi by tato oblast znamenala vyhynutí ryb a v důsledku i vymizení rybářů (tzv. atto-fox problem - podle situace, kdy zbude 10^(-18) lišky - v praxi by opět druh vyhynul).

image.png

ad b) Hledám optimální situaci.

Podržím si graf pomocí hold on a do toho stejného zadám jinou počáteční hodnotu. Zjistil jsem, že pokud zadám hodnotu [2,3], vychází bod => vytvořila se rovnováha (znamená to tedy, že v moři budou plavat dvě ryby a chytat je 3 rybáři):

hold on 

x0=[2,3.01]; 
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2))

Celý kód (včetně hledání dalších možností):

alfa=1.5;
beta=0.5;
gama=0.8;
delta=0.4;
x0=[8,8]; tn=10;
f=@(t,x)[alfa*x(1)-beta*x(1)*x(2); delta*x(1)*x(2)-gama*x(2)];
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2)), xlabel ryby, ylabel rybáři
hold on 

x0=[2,3.01]; 
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2))
hold on

x0=[2,2]; 
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2))
hold on

x0=[3,3]; 
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2))
hold on

x0=[5,2]; 
tic, [t,x] = ode45(f,[0,tn],x0), toc;
plot(x(:,1),x(:,2))
hold on

Graf:

image.png