10.labor

Contents

Egyismeretlenes elsőrendű közönséges differenciálegyenlet

Malthus-féle növekedési modell

Jelöljük a t időpillanatban egy faj egyedeinek számát N(t)-vel. Ekkor N'(t) jelöli a változás sebességet. Ha korlátlan mennyiségben áll rendelkezésre tápanyag, akkor N'(t) egyenesen arányos az adott pillanatban jelenlévő egyedek számával, azaz a differenciálegyenlet: N'(t)=rN(t), ahol r>0 arányossági tényező. Ezen differenciálegyenletet szeretnénk megoldani Matlabban.

function dy=malthus(t,y)
r=0.3;
dy=r*y;

% Az általános megoldó az ode45, ennek használatához szükségünk van a külön
% fájlban lévő differenciálegyenletre, egy kezdeti feltételre, és egy
% időintervallumra, amin meg szeretnénk oldani a differenciálegyenletet.

[T,Y]=ode45(@malthus,[0,20],1);

% Ekkor ábrázolva a megoldást:
plot(T,Y)
% Látjuk, hogy egy egyedből 20 időegység alatt több mint 400 lett. Ez
% történik például sáskajáráskor.

Paraméterek megadása

% Ha a paramétert nem a differenciálegyenlet belsejében szeretnénk megadni
% (mert például különböző paraméterértékekre is meg szeretnénk oldani az
% egyenletet), akkor a function handle-t tudjuk használni. Most úgy írjuk
% fel a függvényt, hogy a paramétertől is függ:

function dy=malthus2(t,y,r)
dy=r*y;

for i=1:5
    r=i/100;
    [T,Y]=ode45(@(t,y)malthus2(t,y,r),[0,20],1);
    plot(T,Y)
    hold on
end
hold off

Verlust-féle logisztikus modell

A Malthus-féle modell nem veszi figyelembe a környezet fenntartó képességét. Ezt például úgy vehetjük figyelembe, a növekedés üteme nem konstans lesz, hanem ha a környezet fenntartó képességét K-val jelöljük, akkor a növekedés üteme rN(t)(1-N(t)/K)-val lesz arányos (azaz ha N(t)>K, akkor negatív, így csökken az egyedszám). Ekkor:

function dy=verlust(t,y)
r=0.3;
K=10;
dy=r*(1-y/K)*y;

% A megoldás:
[T,Y]=ode45(@verlust,[0,20],1);
plot(T,Y)
% Ez is érdekes:
[T2,Y2]=ode45(@verlust,[0,20],21);
hold on
plot(T2,Y2);
hold off

Többismeretlenes elsőrendű közönséges differenciálegyenletrendszer

Kétfajos Lotka-Volterra ragadozó-zsákmány modell

Tegyük fel, hogy két faj él együtt: egy növényevő és egy ragadozó. Szeretnénk az ő kapcsolatukat modellezni. Növényevő önnmagában egy Verlust-típusú modell szerint növelkedne, azaz N(t)-vel jelölve az egyedszámot: N'(t)=r(1-N(t)/K)N(t) viszont ha jelen van egy ragadozó faj is (jelölje P(t) az ő elemszámát): N'(t)=r(1-N(t)/K)N(t)-a*N(t)*P(t) itt alfa az méri, hogy mennyire hatékony ragadozó P. Az általa elejtett zsákmány mennyisége arányos továbbá azzal, hogy mennyi zsákmányállat van jelen, illetve azzal, hogy mennyi ragadozó van jelen. A ragadozó száma zsákmányállat nélkül csökkenne (mert éhenhalnának), de az elejtett zsákmányállatok egy részét gyarapodásra fordítják, így az ő egyenletük: P'(t)=-bP(t)+cN(t)P(t)

% Ezen egyenletrendszer szeretnénk megoldani Matlabban. Ehhez:

function dy = lotvol(t,y)
r=0.5;
K=2;
a=0.1;
b=0.3;
c=0.5;

dy = [ r*(1-y(1)/K)*y(1)-a*y(1)*y(2);...
    -b*y(2)+c*y(1)*y(2)];

% Most a megoldáshoz 2 kezdeti feltételre van szükség:
[T,Y]=ode45(@lotvol,[0,60],[1,1]);
plot(T,Y)

Fázisportré

Kirajzolhatjuk a fázisportrét is, azaz, hogy hogyan alakul az állatok száma egymáshoz képest:

plot(Y(:,1),Y(:,2))

Másodrendű közönséges differenciálegyenlet

Kőhajítás

Tegyük fel, hogy eldobunk egy követ függőlegesen felfelé, jelöljük a kő magasságát t időpontban y(t)-vel. Ekkor hat rá a gravitáció, ami gyorsulást okoz lefelé irányba. Ez azt jelenti, hogy függőleges irányban a gyorsulása: -g (ami nagyjából 9.8 m/s^2). Ekkor a differenciálegyenlet y''(t)=-9.8 lesz. A Matlab azonban csak elsőrendű differenciálegyenleteket tud kezelni, ezért át kell alakítanunk egy elsőrendű rendszerré. Legyen v(t) a kő sebessége a t időpontban. Ekkor y'(t)=v(t), illetve v'(t)=-9.8 (mert a sebesség deriváltja a gyorsulás); azaz vektorosan összefoglalva a helyzetet (első koordináta) és a sebességet (második koordináta):

 function dyv = ko(t,yv)
 %yv(1), a magasság: y(t)
 %yv(2) a függőleges sebesség: v(t)
 dyv = [yv(2); -9.8]; %fontos, hogy oszlopvektorként adjuk meg

% Legyen a kezdősebesség 20 m/s, és a földről indul, azaz 0 magasságról, így
% a megoldó az első 10 másodpercre

[T3,Y3]=ode45(@ko, [0, 10],[0,20]);
plot(T3,Y3(:,1), 'rO', T3,Y3(:,2),'b')

Esemény detektálás

Láthatjuk, hogy a magasság negatív a végén. Ennek nincs értelme, hiszen a kő elvben nem mehet a föld alá, ezért meg kell állítani a megoldót, amikor eléri a 0 magasságot. Ehhez esemény-detektálásra van szükség, azaz kell készíteni egy esemény függvényt:

function [ertek,megalle,irany] = esemeny(t,y)
% Ez a függvény veszi észre, hogy elérjük a 0 magasságot, 
% és megállítja a megoldót ekkor.
ertek = y(1)-0;     % a 0 magasságot keresem, de lehetne más érték is
megalle = 1;   % állítsa-e meg a megoldót
irany = 0;   % csak ha csökken: -1, bármikor amikor eléri: 0, csak ha nő: 1

% Ennek használatához a megoldónak át kell adni az eseményt detektáló
% függvényt, amelyet az odeset parancson keresztül tehetünk meg.

beallitasok = odeset('Events',@esemeny);
[T4,Y4] = ode45(@ko,[0 10],[0,20],beallitasok);
plot(T4,Y4(:,1), 'rO', T4,Y4(:,2),'b')

% Itt amikor 0 lett a magasság, nem számol tovább. A T4 vektorban összegyűjti,
% hogy mely pontban számolta ki (közelítette) a megoldást. T4 utolsó eleme
% az, amikor megállt a megoldó

T4(end)

% ebből láthatjuk, hogy kb 4 másodperc alatt esett vissza a földre a kő.

ans =

    4.0816

Feladatok

1.feladat

A kőhajításnál vegyük figyelembe a légellenállást is, ami jó közelítéssel a sebességgel egyenesen arányos (valamely k arányossági tényezővel), azaz most y''(t)=-9.8-k*y'(t). Ezt átírjuk rendszerré, hasonlóan az előbbiekhez a sebesség lesz az új változó, azaz: y'=v; v'=-9.8-k*v; A k értéke legyen 0.1. Derítsük ki, hogy mennyi idő alatt esik le a kő a földre.

2.feladat

Tegyük fel, hogy két halfaj él egy területen: egy növényevő és egy ragadozó. Ekkor a halászás hatását úgy írhatjuk le, hogy mindkét faj egy hányadát elveszítjük egy adott idő alatt. Ha Lotka-Volterra modellt használunk, akkor mindkét fajnál megjelenik negatív tag; a növényevőknél egy -e*N(t) illetve ragadozóknál egy -e*P(t) tag. Oldjuk meg ezen új egyenletet Matlabban, ha e=0.05. Ábrázoljuk az új rendszer fázisportréját a régivel egy képen.

Published with MATLAB® R2016a