## Systém obyčajných lineárnych diferenciálnych rovníc
1.rádu s konštantnými koeficientmi

### Homogénne systémy

#### Matica systému má len jednoduché
reálne vlastné hodnoty

Postup riešenia si ilustrujeme na systéme

y1'=  y1 + 2y2
y2'= 4y1 + 3y2

Definujeme y1 a y2 ako funkcie premennej x

depends(y1, x)$
depends(y2,x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=y1+2*y2;
DR2:diff(y2,x,1)=4*y1+3*y2;

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1 a y2
a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2)]);

Pomocou príkazu eigenvectors() určíme vlastné hodnoty a vlastné vektory
matice systému A.

Je potrebné si uvedomiť štruktúru výsledku, ktorým je zoznam, na jeho prvom mieste stojí
zoznam obsahujúci zoznam vlastných hodnôt a zoznam ich násobností. Na druhom mieste je zoznam
ktorý obsahuje ďalšie dva zoznamy. Z nich prvý obsahuje vlastný vektor prislúchajúci prvej
vlastnej hodnote (uvedený opäť vpodobe zoznamu) a druhý obsahuje vlastný vektor prislúchajúci
druhej vlastnej hodnote (taktiež v podobe zoznamu).

Túto pomerne zložitú štruktúru je potrebné si uvedomiť pri extrahovaní
požadovaných hodnôt z výsledku.

VL_VEK:eigenvectors(A);

Zo získaných vlastných hodnôt a vektorov zostrojíme všeobecné
riešenie v tvare vektoru

Y:VL_VEK[2][1][1]*%c1*exp(VL_VEK[1][1][1]*x)+VL_VEK[2][2][1]*%c2*exp(VL_VEK[1][1][2]*x);

Z vektora Y môžeme extrahovať jednotlivé zložky riešenia y1 a y2

y1=Y[1];
y2=Y[2];

#### Jendoduché komplexné vlastné hodnoty
matice systému

Postup riešenia si ilustrujeme na systéme

y1'=   y1 + 3y2
y2'= -3y1 +  y2

Definujeme y1 a y2 ako funkcie premennej x

depends(y1, x)$
depends(y2,x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=y1+3*y2;
DR2:diff(y2,x,1)=-3*y1+y2;

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1 a y2
a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2)]);

Pomocou príkazu eigenvectors() určíme vlastné hodnoty a vlastné vektory
matice systému A.

VL_VEK:eigenvectors(A);

Zo získaných vlastných hodnôt a vektorov postupne zostrojíme všeobecné
riešenie.

V danom prípade sú vlastné hodnoty komplexné a preto aj vlastné
vektory majú komplexné súradnice. Pri konštrukcii riešenia preto
budeme musieť oddeliť reálne a komplexné časti vlastných hodnôt
a vlastných vektorov a s ich pomocou riešenie zostaviť.

Pri konštrukcii je zvolené označenie u a v tak, aby bolo zhodné 
s označením na slajdoch o systémoch obyčajných diferenciálnych rovníc.

u:realpart(VL_VEK[2][1][1])*exp(realpart(VL_VEK[1][1][1])*x)*cos(abs(imagpart(VL_VEK[1][1][1]))*x)-imagpart(VL_VEK[2][2][1])*exp(realpart(VL_VEK[1][1][2])*x)*sin(abs(imagpart(VL_VEK[1][1][1]))*x);

v:realpart(VL_VEK[2][1][1])*exp(realpart(VL_VEK[1][1][1])*x)*sin(abs(imagpart(VL_VEK[1][1][1]))*x)+imagpart(VL_VEK[2][2][1])*exp(realpart(VL_VEK[1][1][2])*x)*cos(abs(imagpart(VL_VEK[1][1][1]))*x);

Zo získaného fundamentálneho systému riešení zostrojíme všeobecné
riešenie systému v tvare vektora Y

Y:%c1*u+%c2*v;

Z vektora Y môžeme odseparovať riešenia po zložkách

y1=Y[1];
y2=Y[2];

#### Násobné vlastné hodnoty, vlastné vektory 
lineárne nezávislé

Postupujeme rovnako, ako v prípade jednoduchých vlastných hodnůt.

Postup riešenia si ilustrujeme na systéme

y1'=      y2 + y3
y2'= y1      + y3
y3'= y1 + y2

Definujeme y1, y2 a y3 ako funkcie premennej x

depends(y1, x)$
depends(y2, x)$
depends(y3, x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=y2+y3;
DR2:diff(y2,x,1)=y1+y3;
DR3:diff(y3,x,1)=y1+y2;

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1, y2 
a y3 a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2),coeff(rhs(DR1),y3)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2),coeff(rhs(DR2),y3)],[coeff(rhs(DR3),y1),coeff(rhs(DR3),y2),coeff(rhs(DR3),y3)]);

Pomocou príkazu eigenvectors() určíme vlastné hodnoty a vlastné vektory
matice systému A.

VL_VEK:eigenvectors(A);

Ako vidíme, druhá vlastná hodnota 1 má násobnosť 2, avšak vlastné vektory
príslúchajúce tejto vlastnej hodnote [1,0,-1],[0,1,-1] sú lineárne nezávislé.

Všeobecné riešenie teda zoztavíme rovnakým spôsobom ako v predchádzajúcich prípadoch.

Y:VL_VEK[2][1][1]*%c1*exp(VL_VEK[1][1][1]*x)+VL_VEK[2][2][1]*%c2*exp(VL_VEK[1][1][2]*x)+VL_VEK[2][2][2]*%c3*exp(VL_VEK[1][1][2]*x);

Z vektora Y môžeme odseparovať riešenia po zložkách

y1=Y[1];
y2=Y[2];
y3=Y[3];

#### Násobné vlastné hodnoty, vlastné vektory 
lineárne závislé

Postupnosť lineárne závislých vlastných vektorov musíme nahradiť
reťazcom  zovšeobecnených vlastných vektorov.
(Vysvetlenie pojmu viď. slajdy).

Postup riešenia si ilustrujeme na systéme
y1'= -y1 + y2
y2'=     - y2 + 4y3
y3'= y1       - 4y3

Definujeme y1, y2 a y3 ako funkcie premennej x

depends(y1, x)$
depends(y2, x)$
depends(y3, x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=-y1+y2;
DR2:diff(y2,x,1)=-y2+4*y3;
DR3:diff(y3,x,1)=y1-4*y3;

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1, y2 
a y3 a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2),coeff(rhs(DR1),y3)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2),coeff(rhs(DR2),y3)],[coeff(rhs(DR3),y1),coeff(rhs(DR3),y2),coeff(rhs(DR3),y3)]);

Pomocou príkazu eigenvectors() určíme vlastné hodnoty a vlastné vektory
matice systému A.

VL_VEK:eigenvectors(A);

Z výsledkov vidíme, že dvijnásobnej vlastnej hodnote -3 zodpovedá len jeden vlastný vektor [1,-2,1]. To znamená, že vlastné vektory sú lineáre závislé a musíme nájsť reťazec zovšeobecnených vlatných vektorov. Ten získame a ko riešenie rovnice (A-λE)*v=[1,-2,1].

Najskôr zostavíme maticu B=(A-λE)

E:diagmatrix(3,1)$
B:A-VL_VEK[1][1][1]*E;

Vynásobením matice B stĺpcovým vektorom v=(v1,v2,v3) získame ľavé strany rovníc pre prvý (a v tomto prípade aj posledný) člen reťazca zovšeobecnených vlastných vektorov.

L_STRANY:B.transpose([v1,v2,v3]);

Následne zostavíme systém rovníc pre súradnice vektora v=(v1,v2,v3)

pom_rov1:L_STRANY[1][1]=VL_VEK[2][1][1][1];
pom_rov2:L_STRANY[2][1]=VL_VEK[2][1][1][2];
pom_rov3:L_STRANY[3][1]=VL_VEK[2][1][1][3];

Riešenie systému rovníc pre zovšeobecnený vlastný vektor, existuje nekonečne veľa riešení, ktoré sú závislé od jedného parametra (voľná neznáma)

pom_ries:solve([pom_rov1,pom_rov2,pom_rov3],[v1,v2,v3]);

Dosadením nejakej konkrétnej hodnoty za parameter %r1 (v tomto prípade %r1=1) získame jedno partikulárne riešenie.

POZN. Môže sa stať, že parameter %r1 bol použitý pri riešení iného systému, je potrebné sa pozrieť, aký parameter bol skutočne použitý systémom vo výsledku pom_ries a v prípade potreby ho v nasledujúcej substituúcii hodnoty prepísať.

vl_vek:subst(1,%r1,pom_ries[1]);

Z tohto výsledku extrakciou pravých strán získame prvok zovšeobecneného reťazca vlastných vektorov

ZV_VL_VEK:[rhs(vl_vek[1]),rhs(vl_vek[2]),rhs(vl_vek[3])];

Ten použijeme na zostavenie všeobecného riešenia systému vo vektorovom tvare

Y:4*VL_VEK[2][2][1]*%c1*exp(VL_VEK[1][1][2]*x)+VL_VEK[2][1][1]*%c2*exp(VL_VEK[1][1][1]*x)+VL_VEK[2][1][1]*x*%c3*exp(VL_VEK[1][1][1]*x)+%c3*ZV_VL_VEK*exp(VL_VEK[1][1][1]*x);

Alebo po zložkách

y1=Y[1],expand;
y2=Y[2],expand;
y3=Y[3],expand;

### Nehomogénne systémy

Podobne ako pri rovniciach vyššieho rádu s konštantnými koeficientmi hľadáme riešenie v podobe súčtu všeobecného riešenia homogénneho systému a jedného partikulárneho riešenia systému nehomogénneho.

#### Metóda variácie konštánt

Postup riešenia si ilustrujeme na systéme

y1'= 2y1 - 3y2
y2'=   y1 - 2y2 + 2 sin x

Definujeme y1 a y2 ako funkcie premennej x

depends(y1, x)$
depends(y2,x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=2*y1-3*y2;
DR2:diff(y2,x,1)=y1-2*y2+2*sin(x);

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1 a y2
a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2)]);

Pomocou príkazu eigenvectors() určíme vlastné hodnoty a vlastné vektory
matice systému A.

VL_VEK:eigenvectors(A);

Zo získaných vlastných hodnôt a vektorov zostrojíme všeobecné
riešenie v tvare vektoru

Y:VL_VEK[2][1][1]*%c1*exp(VL_VEK[1][1][1]*x)+3*VL_VEK[2][2][1]*%c2*exp(VL_VEK[1][1][2]*x);

Z vektora Y môžeme extrahovať jednotlivé zložky všeobecného riešenia y1v a y2v

y1v=Y[1];
y2v=Y[2];

Zo všeobecného riešenia zostavíme fundamentálnu maticu systému, ktorú získame z koeficientov pri konštantách %c1 a %c2

U(x):=matrix([coeff(Y[1],%c2),coeff(Y[1],%c1)],[coeff(Y[2],%c2),coeff(Y[2],%c1)])$
U(x);

Pre určenie partikulárneho riešenia metódou variácie konštánt budeme potrebovať inverznú maticu ku fundamentálnej matici. Označíme ju U_inv(x)

U_inv(x):=invert(U(x))$
U_inv(x);

Definujeme vektor f(x) ako nehomogenitu systému

f_vek(x):=matrix([0],[2*sin(x)]);

Partikulárne riešenie určíme podľa vzorca variácie konštánt ako U(x)*∫U_inv(x)*f(x)dx

Y_part:U(x).integrate(U_inv(x).f_vek(x),x),ratsimp;

Riešenie nehomogénneho systému potom určíme ako súčet všeobecného riešenia Y a partikulárneho riešenia Y_part

Y_nehom:Y+Y_part;

Alebo riešenie extrahujeme z vektorového tvaru po zložkách

y1=Y_nehom[1];
y2=Y_nehom[2];

#### Špeciálny tvar pravej strany

Špeciálny tvar pravej strany je podobný špeciálnym tvarom pravej strany pri rovniciach vyššieho rádu, teda môže byť v tvare P(x), e^(αx)*P(x) alebo e^(αx)*sin(βx)*P(x) resp. e^(αx)*cos(βx)*P(x), kde P(x) je polynóm.

POZN. Vzhľadom na to, že ide o systém rovníc, má pravá strana f(x) vektorový tvar a uvedené tvary funkcií treba chápať tak, že platia po zložkách, teda napr. ak hovoríme, že P(x) je polynóm, znamená to, že každá zložka vektora f(x) je polynóm a pod.

Postup riešenia si ilustrujeme na systéme
y1'= -2y1 + 2y2 + 5*sin(2*x)
y2'=   y1 - 2y2 + 2 sin x

Definujeme y1 a y2 ako funkcie premennej x

depends(y1, x)$
depends(y2,x)$

Zostavíme zadaný systém diferenciálnych rovníc

DR1:diff(y1,x,1)=-2*y1+2*y2 + 5*sin(2*x);
DR2:diff(y2,x,1)=2*y1-5*y2;

Pomocou príkazov coeff extrahujeme z pravých strán koeficienty pri y1 a y2
a zostavíme z nich maticu systému

A:matrix([coeff(rhs(DR1),y1),coeff(rhs(DR1),y2)],[coeff(rhs(DR2),y1),coeff(rhs(DR2),y2)]);

VL_VEK:eigenvectors(A);

Zo získaných vlastných hodnôt a vektorov zostrojíme všeobecné
riešenie v tvare vektoru

Y_vseob:VL_VEK[2][1][1]*%c1*exp(VL_VEK[1][1][1]*x)+2*VL_VEK[2][2][1]*%c2*exp(VL_VEK[1][1][2]*x);

Z vektora Y_vseob môžeme extrahovať jednotlivé zložky všeobecného riešenia y1v a y2v

y1v=Y_vseob[1];
y2v=Y_vseob[2];

Definujeme partikulárne riešenie v predpokladanom tvare.

POZN. Definujeme ako maticu 2x1 čo zodpovedá stĺpcovému vektoru, je však potrebné sa na zložky odvolávať dvojicou indexov (ako v matici). Z výhodou potom môžeme využívať maticové operácie.

YP(x):=matrix([a1*cos(2*x)+b1*sin(2*x)],[a2*cos(2*x)+b2*sin(2*x)]);

Pre účely porovnávania koeficientov určíme tvar ľavých strán systému po dosadení partikulárneho riešenia

LS:diff(YP(x),x,1);

Definujeme vektor pravých strán (nehomogenitu systému).

f(x):=matrix([5*sin(2*x)],[0]);

Pre účely porovnávania koeficientov určíme tvar pravých strán systému po dosadení partikulárneho riešenia

PS:A.YP(x)+f(x),expand;

Porovnaním koeficientov pri rovnakých členoch na ľavých a pravých stranách ayatému získame systém rovníc pre hľadané neurčité koeficienty

rov1:coeff(LS[1][1],cos(2*x))=coeff(PS[1][1],cos(2*x));
rov2:coeff(LS[2][1],cos(2*x))=coeff(PS[2][1],cos(2*x));
rov3:coeff(LS[1][1],sin(2*x))=coeff(PS[1][1],sin(2*x));
rov4:coeff(LS[2][1],sin(2*x))=coeff(PS[2][1],sin(2*x));

Vyriešime daný systém rovníc

koeficienty:solve([rov1,rov2,rov3,rov4],[a1,b1,a2,b2]);

Výsledky dosadíme za koeficienty do predpokladaného tvaru partikulárneho riešenia

Y_part:subst(rhs(koeficienty[1][1]),a1,subst(rhs(koeficienty[1][2]),b1,subst(rhs(koeficienty[1][3]),a2,subst(rhs(koeficienty[1][4]),b2,YP(x)))));

Riešenie nehomogénneho systému určíme ako súčet všeobecného a partikulárneho riešenia

Y:Y_vseob+Y_part;

Z vektora Y môžeme extrahovať jednotlivé zložky riešenia y1 a y2

y1=Y[1][1];
y2=Y[2][1];