Normaaliryhmä

Seuraavaksi luodaan katsaus useammasta differentiaaliyhtälöstä muodostuvaan differentiaaliyhtälöryhmään. Esimerkkien avulla pohditaan, kuinka käytännön ongelmasta muokataan normaalimuotoinen differentiaaliyhtälöryhmä ja kuinka se voidaan ratkaista numeerisesti Matlabilla. Varsinainen lineaaristen differentiaaliyhtälöryhmien teoria jätetään tulevalle opintojaksolle.

Määritelmä.

Olkoot \(x_1,x_2,\ldots,x_n\) tuntemattomia reaalifunktioita. Ensimmäisen kertaluvun differentiaaliyhtälöistä koostuvaa yhtälöryhmää

\[\begin{split}\begin{cases} x_1'(t)=f_1(t,x_1(t),x_2(t),\ldots,x_n(t)),\\ x_2'(t)=f_2(t,x_1(t),x_2(t),\ldots,x_n(t)),\\ \mspace{47mu}\vdots\\ x_n'(t)=f_n(t,x_1(t),x_2(t),\ldots,x_n(t)), \end{cases}\end{split}\]

missä funktiot \(f_i : I\times\mathbb R^n\to\mathbb R\) ovat \(n+1\) muuttujan reaaliarvoisia funktioita, kutsutaan normaaliryhmäksi (first-order system). Välillä \(I\) derivoituvat funktiot \(x_1,x_2,\ldots,x_n\) muodostavat normaaliryhmän ratkaisun välillä \(I\), jos yhtälöryhmän yhtälöt ovat voimassa aina, kun \(t\in I\).

Normaaliryhmään liittyvällä alkuarvotehtävällä tarkoitetaan sellaisen ratkaisun hakemista, joka toteuttaa ennalta asetetut ehdot

\[x_1(t_0)=b_1,\qquad x_2(t_0)=b_2,\qquad\ldots\qquad x_n(t_0)=b_n\]

välin \(I\) pisteessä \(t_0\).

Esimerkki.

Differentiaaliyhtälöryhmä

\[\begin{split}\begin{cases} x'=2x+y,\\ y'=x+2y-e^{2t}, \end{cases}\end{split}\]

on funktioiden \(x\) ja \(y\) normaaliryhmä, jonka määrittelevät funktiot \(f_1(t,x,y)=2x+y\) ja \(f_2(t,x,y)=x+2y-e^{2t}\). Ryhmän eräs ratkaisu on

\[\begin{split}\begin{cases} x(t)=e^{2t}-e^t,\\ y(t)=e^t, \end{cases}\end{split}\]

sillä näille funktioille on voimassa

\[\begin{split}\begin{cases} x'(t)=2e^{2t}-e^t=2\left(e^{2t}-e^t\right)+e^t=2x(t)+y(t),\\ y'(t)=e^t=\left(e^{2t}-e^t\right)+2e^t-e^{2t}=x(t)+2y(t)-e^{2t}. \end{cases}\end{split}\]

Pieni normaaliryhmä (\(n=2\) tai \(3\)) voidaan yrittää ratkaista palauttamalla se yhdeksi korkeamman kertaluvun yhtälöksi. Suurempien ryhmien tapauksessa tämä menetelmä ei useinkaan ole käyttökelpoinen (katso huomatus normaaliryhmistä).

Esimerkki.

Ratkaise alkuarvotehtävä

\[\begin{split}\begin{cases} x_1'=4x_1-3x_2,&(1)\\ x_2'=6x_1-7x_2,&(2) \end{cases} \qquad \begin{cases} x_1(0)=2\\ x_2(0)=-1. \end{cases}\end{split}\]

Ratkaisu.

Menetellään seuraavasti.

1. Ratkaistaan \(x_1\) yhtälöstä (2) muuttujien \(x_2'\) ja \(x_2\) funktiona ja lasketaan sen derivaatta muuttujien \(x_2''\) ja \(x_2'\) funktiona.

   \[\begin{split}\begin{cases} x_1 = \frac{1}{6}x_2' + \frac{7}{6}x_2, & (3) \\ x_1' = \frac{1}{6}x_2'' + \frac{7}{6}x_2', & (4) \end{cases}\end{split}\]

2. Sijoitetaan (3) ja (4) yhtälöön (1). Saadaan toisen kertaluvun yhtälö funktiolle \(x_2\), josta \(x_2\) ratkaistaan.

   \[\frac16x_2''+\frac76x_2'=4\left(\frac16x_2'+\frac76x_2\right)-3x_2 \Leftrightarrow x_2''+3x_2'-10x_2=0\]

   Tämä on toisen kertaluvun lineaarinen homogeeninen yhtälö funktiolle \(x_2\), jonka karakteristisen yhtälön \(\lambda^2+3\lambda-10=0\) juuret ovat \(\lambda=2\) ja \(\lambda=-5\). Yleinen ratkaisu on siis

   \[x_2=c_1e^{2t}+c_2e^{-5t}.\qquad(5)\]

3. Nyt \(x_1\) saadaan sijoittamalla (5) yhtälöön (3).

   \[x_1=\frac16\left(2c_1e^{2t}-5c_2e^{-5t}\right)+\frac76\left(c_1e^{2t}+c_2e^{-5t}\right)=\frac{3}{2}c_1e^{2t}+\frac{1}{3}c_2e^{-5t}\qquad(6)\]

Yhtälöt (5) ja (6) ovat nyt alkuperäisen normaaliryhmän yleinen ratkaisu. Muistetaan alkuehdot

\[\begin{split}\begin{cases} x_1(0)=\frac{3c_1}{2}+\frac{c_2}{3}=2\\ x_2(0)=c_1+c_2=-2, \end{cases}\end{split}\]

joista ratkaistaan \(c_1=2\) ja \(c_2=-3\). Siispä alkuarvotehtävän ratkaisu on

\[\begin{split}\begin{cases} x_1(t)=3e^{2t}-e^{-5t}\\ x_2(t)=2e^{2t}-3e^{-5t}. \end{cases}\end{split}\]

Määritelmä.

Muotoa

\[x^{(n)}(t)=f\left(t,x(t),x'(t),x''(t),\ldots,x^{(n-1)}(t)\right)\]

oleva differentiaaliyhtälö on funktion \(x=x(t)\) \(n\). kertaluvun normaalimuotoinen differentiaaliyhtälö.

Normaalimuotoinen yhtälö voidaan palauttaa normaaliryhmäksi asettamalla

\[x_1=x,\qquad x_2=x',\qquad x_3=x'',\qquad \ldots\qquad x_n=x^{(n-1)},\]

sillä tällöin

\[x_1'=x'=x_2,\qquad x_2'=x''=x_3,\qquad \ldots\qquad x_n'=x^{(n)}.\]

Nyt normaalimuotoisen differentiaaliyhtälön ratkaisemiseksi riittää ratkaista normaaliryhmä

\[\begin{split}\begin{cases} x_1' = x_2 \\ x_2' = x_3 \\ \phantom{x_3'}\vdots \\ x_{n-1}' = x_n \\ x_n' = f(t, x_1, x_2, \ldots, x_n), \end{cases}\end{split}\]

ja tämän ratkaisusta poimitaan \(x(t)=x_1(t)\).

Esimerkki.

Muokkaa alkuarvotehtävä

\[x'''=-(x')^2-3x+e^t,\qquad x(0)=1,\ x'(0)=0,\ x''(0)=0\]

normaaliryhmäksi.

Ratkaisu.

Asetetaan

\[x_1=x,\qquad x_2=x'\qquad\text{ja}\qquad x_3=x'',\]

jolloin normaaliryhmäksi tulee

\[\begin{split}\begin{cases} x_1' = x_2 \\ x_2' = x_3 \\ x_3' = -x_2^2 - 3x_1 + e^t \end{cases}\qquad \begin{cases} x_1(0)=1\\ x_2(0)=0\\ x_3(0)=0. \end{cases}\end{split}\]

Myös useamman normaalimuotoisen differentiaaliyhtälön ryhmä voidaan palauttaa normaaliryhmäksi.

Esimerkki.

Tarkastellaan oheisen kuvan mukaista kahden jousen systeemiä.

Paikkojen \(x_1(t)\) ja \(x_2(t)\) nollakohdat on kiinnitetty siten, että tasapainotilassa \(x_1=x_2=0\). Massaan \(m_1\) vaikuttaa jousen 1 aiheuttaman jousivoiman \(-k_1x_1\) lisäksi jousen 2 aiheuttama voima \(k_2(x_2-x_1)\), missä \(x_2-x_1\) on jousen 2 venymä. Jousen 2 aiheuttama voima osoittaa siis oikealle, jos jousi 2 on venynyt eli \(x_2-x_1>0\). Massaan \(m_2\) vaikuttaa tämän vastavoima eli \(-k_2(x_2-x_1)\). Systeemin käyttäytymistä kuvaa siten differentiaaliyhtälöryhmä (vertaa vaimennettuun vapaaseen värähtelijään)

\[\begin{split}\begin{cases} m_1x_1''=-k_1x_1+k_2(x_2-x_1)\\ m_2x_2''=-k_2(x_2-x_1). \end{cases}\end{split}\]

Muunnetaan tämä normaaliryhmäksi merkitsemällä \(x_3=x_1'\) ja \(x_4=x_2'\), jolloin saadaan

\[\begin{split}\begin{cases} x_1'=x_3\\ x_2'=x_4\\ x_3'=(-k_1x_1+k_2(x_2-x_1))/m_1\\ x_4'=(-k_2(x_2-x_1))/m_2. \end{cases}\end{split}\]

Huomautus.

Miksi normaaliryhmä? Ensimmäisen kertaluvun differentiaaliyhtälöiden ja normaaliryhmien ratkaisemiseksi on tehokkaita numeerisia menetelmiä. Käytännössä differentiaaliyhtälö tai yhtälöryhmä usein palautetaankin normaaliryhmäksi ja ratkaistaan numeerisesti. Muuttujien ja yhtälöiden lukumäärä voi sovelluksissa hyvin olla tuhansia.

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

\[\mathbf{f}(t,\mathbf{x})=(f_1(t,\mathbf{x}),\ldots,f_n(t,\mathbf{x})).\]

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

\[\mathbf{x}(t_0)=(b_1,\ldots,b_n).\]

Tutkitaan esimerkiksi alkuarvotehtävää

\[\begin{split}\begin{cases} x_1'=x_1+tx_2,\\ x_2'=x_1+5x_2+e^t, \end{cases} \qquad \begin{cases} x_1(-1)=1{,}5,\\ x_2(-1)=-0{,}5. \end{cases}\end{split}\]

Funktio \(\mathbf{f}\) on

\[\mathbf{f}(t,\mathbf{x})=\left(x_1+tx_2,\,x_1+5x_2+e^t\right).\]

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.

Ratkaisu.

1. Säiliöön 1 tulee suolaa

   \[\dfrac{y\text{ kg}}{200\text{ l}}\cdot10\ \dfrac{\text{l}}{\text{min}}=0{,}05y\ \frac{\text{kg}}{\text{min}}\]

   ja sieltä poistuu suolaa

   \[\dfrac{x\text{ kg}}{100\text{ l}}\cdot30\ \dfrac{\text{l}}{\text{min}}=0{,}3x\ \frac{\text{kg}}{\text{min}}.\]

   Koska \(x'(t)\) on suolan määrän muutosnopeus (kg/min) säiliössä 1 hetkellä \(t\), niin saadaan yhtälö \(x'=-0{,}3x+0{,}05y\). Vastaavalla tavoin säiliötä 2 tarkastelemalla päätellään \(y'=0{,}3x-0{,}15y\).

2. On ratkaistava alkuarvotehtävä

   \[\begin{split}\begin{cases} x'&=-0{,}3x+0{,}05y\\ y'&=0{,}3x-0{,}15y \end{cases} \qquad \begin{cases} x(0)=0{,}2\\ y(0)=0{,}1. \end{cases}\end{split}\]

   Ratkaistaan systeemi numeerisesti Matlabin ode45-funktiolla ja piirretään kuva.

   odefun=@(t,x)[-0.3*x(1)+0.05*x(2); 0.3*x(1)-0.15*x(2)];
   [t,x]=ode45(odefun,[0 80],[0.2 0.1]);
   plot(t,x(:,1),t,x(:,2),'--')
   set(gca,'fontsize',14)
   legend('x(t)','y(t)')
   xlabel('t')
   ylabel('x,y')
   

   
   align:center

   Säiliössä 1 on aluksi 200 g suolaa ja määrä vähenee monotonisesti kohti nollaa. Säiliön 2 suolapitoisuus aluksi nousee ja vähenee sitten kohti nollaa. Vaikuttaako ratkaisu järkevältä?

   Yhtälöpari voidaan ratkaista myös symbolisesti dsolve-komennolla.

   s=dsolve('Dx=-0.3*x+0.05*y, Dy=0.3*x-0.15*y, x(0)=0.2, y(0)=0.1')
   s.x
   s.y
   

width:50.0% width:50.0% width:60.0% width:60.0%

Posting submission...