Normaaliryhmä¶

Normaaliryhmän ratkaiseminen numeerisesti Matlabilla¶

Merkitään \(\mathbf{x}=(x_1,\ldots,x_n)\in\mathbb R^n\). Tällöin normaaliryhmän funktioita voidaan merkitä lyhyesti \(f_i(t,x_1,\ldots,x_n)=f_i(t,\mathbf{x})\). Muodostetaan reaaliarvoisista funktioista \(f_i\) vektoriarvoinen funktio

Myös alkuehdot \(x_1(t_0)=b_1,\ldots,x_n(t_0)=b_n\) voidaan antaa vektorimuodossa, kuten

Tutkitaan esimerkiksi alkuarvotehtävää

Funktio \(\mathbf{f}\) on

Matlabissa tämä voidaan syöttää seuraavasti.

f=@(t,x)[x(1)+t*x(2);x(1)+5*x(2)+exp(t)]

Tässä x on vektori, jonka komponentteihin viitataan x(1) ja x(2). Jos halutaan välillä \([t_0,t_1]\) numeerinen arvio ratkaisulle, joka toteuttaa alkuehdon \(\mathbf{x}(t_0)=(b_1,b_2)\), voidaan käyttää Matlabin ode45-funktiota.

[t,x] = ode45(f,[t0 t1],[b1 b2]);

Esimerkissämme alkuehto on \(\mathbf{x}(-1)=(1{,}5,\ -0{,}5)\), joten

[t,x] = ode45(f,[-1 0],[1.5 -0.5]);

antaa approksimaation ratkaisulle välillä \([-1,0]\). Nyt vektorissa t on järjestyksessä muuttujan \(t\) arvoja väliltä \([-1,0]\). Matriisin x ensimmäisessä sarakkeessa x(:,1) on vastaavilla paikoilla arviot arvoille \(x_1(t)\) ja toisessa sarakkeessa x(:,2) arviot arvoille \(x_2(t)\). Funktioiden \(x_1(t)\) ja \(x_2(t)\) kuvaajat voidaan siten piirtää.

plot(t,x(:,1),t,x(:,2))

Matlabin ohjeissa funktiolle \(\mathbf{f}\) käytetään nimeä odefun.

Esimerkki.

Ratkaistaan aiempi esimerkki Matlabilla. Symbolista ratkaisua ei voida muodostaa.

x=dsolve('D3x=-(Dx)^2-3*x+exp(t),x(0)=1,Dx(0)=0,D2x(0)=0','t')
Warning: Explicit solution could not be found.
x = [ empty sym ]

Ratkaistaan vastaava normaaliryhmä numeerisesti välillä \(t\in[0,3]\) funktiolla ode45 ja tulostetaan ratkaisun \(x(t)=x_1(t)\) kuvaaja:

% Normaaliryhmän oikea puoli t:n ja x:n funktiona.
% x on vektori, jonka komponentteihin viitataan x(1), x(2), x(3)
odefun=@(t,x)[x(2);x(3);-x(2)^2-3*x(1)+exp(t)]

t0 = 0;      % alkuajanhetki
t1 = 3;      % loppuajanhetki
x0 = [1 0 0] % alkuehdot
[t,x] = ode45(odefun,[t0 t1],x0);
% t:n arvot ovat vektorissa t
% xi:n likiarvot ovat matriisin x i:nnessä sarakkeessa
plot(t,x(:,1)) % x1:n kuvaaja

Esimerkki.

Ratkaistaan aiemman esimerkin kahden jousen systeemi Matlabilla tapauksessa \(m_1=m_2=1\), \(k_1=1\) ja \(k_2=2\) välillä \(t\in[0,20]\) funktiolla ode45 ja tulostetaan ratkaisujen \(x_1(t)\) ja \(x_2(t)\) kuvaajat.

% Normaaliryhmän oikea puoli t:n ja x:n funktiona.
% x on vektori, jonka komponentteihin viitataan
% x(1), x(2), x(3) ja x(4).
odefun=@(t,x)[x(3);x(4);-x(1)+2*(x(2)-x(1));-2*(x(2)-x(1))]

t0 = 0;        % alkuajanhetki
t1 = 20;       % loppuajanhetki
x0 = [1 0 0 0] % alkuehdot
[t,x] = ode45(odefun,[t0 t1],x0);
% t:n arvot ovat vektorissa t
% xi:n likiarvot ovat matriisin x i:nnessä sarakkeessa
plot(t,x(:,1),t,x(:,2),'--') % x1:n ja x2:n kuvaajat

Tämän tehtävän Matlab osaa ratkaista myös symbolisesti.

s=dsolve('D2x=-x+2*(y-x),D2y=-2*(y-x)',...
'x(0)=1,y(0)=0,Dx(0)=0,Dy(0)=0','t')
s.x
s.y

Tarkastellaan peto-saalis-järjestelmää (predator-prey system), jossa on yksi petolaji ja sillä yksi saalislaji, esimerkiksi pöllöt ja myyrät. Kun saaliin määrä kasvaa, niin petojen lukumääräkin kasvaa. Toisaalta petojen määrän kasvu hidastaa saalispopulaation kasvua ja lopulta kääntää saaliin määrän laskuun aiheuttaen ennen pitkää myös petopopulaation pienenemisen. Analysoidaan tätä kilpajuoksua tarkemmin. Olkoon \(x(t)\) saaliseläinten lukumäärää ja \(y(t)\) petojen lukumäärää. Tarkastellaan ensin kahta ääritapausta 1 ja 2 ja sitten yleistä tapausta 3.

1. Jos petoja ei ole (\(y=0\)), niin saalispopulaatio kasvaa eksponentiaalisesti, kuten

   \[x'=ax,\qquad a>0.\]

   (katso esimerkki populaation koosta).

2. Jos saalista ei ole (\(x=0\)), niin petopopulaatio vähenee eksponentiaalisesti, kuten

   \[y'=-by, \qquad b>0.\]

3. Jos \(x,y>0\), niin kohdassa 1 petojen määrän kasvaminen vähentää saalispopulaation kasvunopeutta \(x'\). Oletetaan, että vähennys kasvunopeuteen on suoraan verrannollinen petojen ja saaliseläinten kohtaamistaajuuteen, ja oletetaan edelleen, että kohtaamistaajuus on suoraan verrannollinen tuloon \(xy\). Eli jos esimerkiksi petojen määrä kaksinkertaistuu, niin kohtaamistaajuus myös kaksinkertaistuu, tai jos esimerkiksi saaliin määrä puolittuu, niin kohtaamistaajuuskin puolittuu. Niinpä päädytään yhtälöön

   \[x'=ax-pxy\qquad a,p>0.\]

   Vastaavasti kohdassa 2 saaliin määrän lisääntyminen kasvattaa petopopulaation kasvunopeutta \(y'\) suoraan verrannollisesti kohtaamistaajuuteen. Niinpä

   \[y'=-by+qxy\qquad b,q>0.\]

   Nämä kaksi yhtälöä muodostavat nyt (yksinkertaistetun) Lotka-Volterran mallin peto-saalis-järjestelmästä.

Esimerkki.

Olkoon nummella hetkellä \(t=0\) (kuukausina) \(70\) kania ja \(40\) kettua. Olkoon ajan hetkellä \(t\) kanien lukumäärää \(x(t)\) ja kettujen lukumäärää \(y(t)\). Oletetaan, että lukumäärät noudattavat mallia

\[\begin{split}\begin{cases} x'=0{,}2x-0{,}005xy\\ y'=-0{,}5y+0{,}01xy. \end{cases}\end{split}\]

Ratkaistaan differentiaaliyhtälöpari numeerisesti ja piirretään kuva ratkaisuista \(x(t)\) ja \(y(t)\) Matlabilla.

odefun=@(t,x)[0.2*x(1)-0.005*x(1)*x(2); -0.5*x(2)+0.01*x(1)*x(2)]
[t,x]=ode45(odefun,[0:0.1:60],[70,40]);
figure
plot(t,x(:,1),t,x(:,2),'--')
set(gca,'fontsize',14)
hleg1=legend('x(t) kanit', 'y(t) ketut')
set(hleg1,'Location','SouthEast')
xlabel('t')
ylabel('x,y')
axis image

Nähdään, että populaatioiden koot \(x(t)\) ja \(y(t)\) ovat jaksollisia ja että syklin pituus (kahden huipun välinen etäisyys) on noin \(20\) kk ja kettujen maksimimäärä saavutetaan noin \(5\) kk kanien maksimimäärän jälkeen.

Esimerkki.

Tarkastellaan niin sanottua kahden säiliön ongelmaa. Säiliöissä 1 on 100 litraa ja säiliössä 2 on 200 litraa suolavettä. Olkoon \(x(t)\) suolan määrää säiliössä 1 ja \(y(t)\) suolan määrää (kg) säiliössä 2 ajan \(t\) (min) funktiona. Säiliöön 1 johdetaan suolatonta vettä 20 l/min ja säiliöstä 2 poistuu suolavettä 20 l/min. Lisäksi säiliöstä 1 johdetaan suolavettä säiliöön 2 30 l/min ja säiliöstä 2 säiliöön 1 10 l/min (eli säiliöissä 1 ja 2 on jatkuvasti 100 ja 200 litraa suolavettä). Oletetaan lisäksi, että kummassakin säiliössä suola on sekoittamisen ansiosta kunakin ajanhetkenä tasaisesti jakautunut.

1. Kirjoita differentiaaliyhtälöpari funktioille \(x(t)\) ja \(y(t)\).
2. Tarkastele tilannetta siinä tapauksessa, että alkuhetkellä (\(t=0\)) säiliössä 1 on 200 g suolaa ja säiliössä 2 on 100 g suolaa.