# Lineárna diferenciálna rovnica n-tého
rádu s konštantnými koeficientmi

## Homogénna rovnica

Všeobecné riešenie homogénnej rovnice určíme pomocou koreňov
charakteristickej rovnice. Tú zostavíme tak, že derivácie hľadanej
funkcie nahradíme mocninami stupňa, ktorý je rovný rádu derivácie.
Funkciu samotnú chápeme ako nultú deriváciu

### Charakteristická rovnica má jednoduché reálne korene

Postup ilustrujeme na úlohe:

Riešime rovnicu

y''-5y'+6y=0 

Ako vždy, si samozrejme najskôr definujeme závislosť premenných
a zostavíme rovnicu

depends(y, x)$
DR: diff(y,x,2) - 5*diff(y,x) + 6*y = 0;

Zostrojíme charakteristickú rovnicu. To dosiahneme trojnásobnou
substitúciou za y'',y' a y.
Najskôr definujeme tri substitúcie mocnín za derivácie, potom ich postupne zrealizujeme.
zrealizujeme. Pri tom sa šikovne využije hodnota %, ktorá predstavuje výsledok posledného
príkazu. Pripomínam, že symbolom $ potláčame výstupy medzikrokov.

SUB1: p^2 = diff(y,x,2)$
SUB2: p = diff(y,x)$
SUB3: 1 = y$
subst(lhs(SUB1), rhs(SUB1), DR)$
subst(lhs(SUB2), rhs(SUB2), %)$
CHAR: subst(lhs(SUB3), rhs(SUB3), %);

Z koreňov charakteristicke rovnice vyskladáme všeobecné riešenie.
Všimnite si, že riešenie rovnice druhého rádu je závislé od dvoch konštánt.
Všeobecne paltí, že aký je rád rovnice, od toľkých konštánt závisí jej riešenie.

CHAR_KOREN1: solve(CHAR, p)[1]$
CHAR_KOREN2: solve(CHAR, p)[2]$
VSEOB_R: y = %k1*exp(rhs(CHAR_KOREN1)*x) + %k2*exp(rhs(CHAR_KOREN2)*x);

Získané riešenie VSEOB_R môžeme porovnať s výsledkom, ktorý získame
priamym použitím príkazu ode2 na zadanú rovnicu

ode2(DR, y, x);
method;

Ako vidíme, riešenie je zhodné s naším výsledkom, príkaz method
prezrádza postup, ktorý bol pri riešení použitý.

### Jednoduché komplexné korene

Riešime rovnicu

y''+6y'+13y=0 

Definujeme rovnicu (závislosť už máme z predchádzajúcej úlohy, 
ak ste ju preskočili, je potrebné pridať príkaz
depends(y, x)$)

DR: diff(y,x,2) + 6*diff(y,x) + 13*y = 0;

Opäť zrealizujeme substitúcie mocnín za derivácie
aby sme zostavili charakteristickú rovnicu.

SUB1: p^2 = diff(y,x,2)$
SUB2: p = diff(y,x)$
SUB3: 1 = y$
subst(lhs(SUB1), rhs(SUB1), DR)$
subst(lhs(SUB2), rhs(SUB2), %)$
CHAR: subst(lhs(SUB3), rhs(SUB3), %);

Pretože máme reálne koeficienty v rovnici, má táto komplexne
združené korene. V danom prípade nám preto stačí určiť iba jeden.
Z neho potom oddelením reálnej a imaginárnej časti zostavíme
všeobecné riešenie závislé od dvoch konštánt

CHAR_KOREN1: solve(CHAR, p)[1]$
VSEOB_R: y = %k1*exp(realpart(rhs(CHAR_KOREN1))*x)*cos(abs(imagpart(rhs(CHAR_KOREN1)))*x) + %k2*exp(realpart(rhs(CHAR_KOREN1))*x)*sin(abs(imagpart(rhs(CHAR_KOREN1)))*x);

Opäť môžeme porovnať zhodnosť výsledku s priamym použitím príkazu ode2()

ode2(DR, y, x);
method;

### Viacnásobné reálne korene

Riešime rovnicu

y''''-4y'''-2y''+12y'+9y=0 

Definujeme závislosť y od x (ak už ju nemáme definovanú z predchádzajúcich úloh)

depends(y,x);

Zostavíme dif. rovnicu podľa zadania

DR: diff(y,x,4) -4*diff(y,x,3)- 2*diff(y,x,2) +12*diff(y,x,1)+ 9*y = 0;

Vykonáme substitúcie mocnín namiesto derivácií a zostavíme
charakteristickú rovnicu

SUB1: p=diff(y,x)$
SUB2: p^2 = diff(y,x,2)$
SUB3: p^3=diff(y,x,3)$
SUB4: p^4 = diff(y,x,4)$
SUB5: 1 = y$
subst(lhs(SUB1), rhs(SUB1), DR)$
subst(lhs(SUB2), rhs(SUB2), %)$
subst(lhs(SUB3), rhs(SUB3), %)$
subst(lhs(SUB4), rhs(SUB4), %)$
CHAR: subst(lhs(SUB5), rhs(SUB5), %);

Nájdeme jej korene a pomocou príkazu multiplicities
určíme aj ich násobnosti

CHAR_KORENE:solve(CHAR,p);
nasobnosti:multiplicities;

Vidíme, že charakteristická rovnica mí dva dvojnásobné korene,
ku každému z nich teda prislúchajú dve lineárne nezávislé riešenia,
ich lineárnouo kombináciou zostavíme všeobecné riešenie.
Tu je vyskladané na dvakrát, najskôr pre každý z dvojnásobných
koreňov, následne zlúčené do celého všeobecného riešenia

VSEOB_R1: %k1*exp(rhs(CHAR_KORENE[1])*x) + %k2*x*exp(rhs(CHAR_KORENE[1])*x)$
VSEOB_R2: %k3*exp(rhs(CHAR_KORENE[2])*x) + %k4*x*exp(rhs(CHAR_KORENE[2])*x)$
VSEOB_R: y = VSEOB_R1+VSEOB_R2;

Príkaz ode2 už nám neposkytne riešenie, rieši iba rovnice
prvého a druhého rádu. 

ode2(DR,y,x);
method;

BONUS: Pokúste sa zostavenie všeobecného riešenia zautomatizovať
v podobe nejakej funkcie, ktorá na základe zoznamu koreňov 
a zoznamu  ich násobností vyskladá všeobecné riešenie.

## Nehomogénne rovnice

Riešenie nehomogénnej rovnice hľadáme v tvare súčtu všeobecného riešenia
príslušnej homogenizovanej rovnice a nejakého partikulárneho riešenia
rovnice nehomogénnej.
Riešenie nehomogénnej rovnice hľadáme metódou variácie konštánt.
V niektorých špeciálnych prípadoch môžeme postupovať tiež metódou
neurčitých koeficientov.

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

Metódu ilustrujeme na riešení rovnice

y''-y'=e^x/(e^x+1)

Najskôr definujeme závislosť y od nezávisle premennej x (ak sme tak už neurobili
v predchádzajúcich úlohách) a zostavíme diferenciálnu rovnicu.

depends(y, x)$
DR: diff(y,x,2)-diff(y,x)= exp(x)/(exp(x)+1);

Zo zadania vytvoríme homogenizovanú rovnicu DR_HOMOG s nulovou pravou
stranou. Jej všeobecné riešenie VR_HOMOG získame pomocou príkazu ode2().
Pri rovniciach rádu vyššieho než 2 opakuejme ručný postup z časti
o homogénnej rovnici. 

DR_HOMOG: lhs(DR) = 0$
VR_HOMOG: ode2(DR_HOMOG, y, x);

Fundamentálny systém riešení (pojem viď. slajdy) uložíme ako Y1 a Y2.
Následne definujeme závislosť "konštánt" k1 a k2 ako funkcií závislých
od premennej x (variácia konštánt) a zostavíme systém rovníc pre neznáme
funkcie K1 a K2 (štruktúra systému viď. slajdy)

UPOZORNENIE:
Pri pokuse o opätovné použitie preddefinovaných príkazov na inú rovnicu
môže dôjsť ku kolízii definície závislosti K1 a K2 od x vzhľadom na už skôr
získané hodnoty K1 a K2 (zo skôr riešeného príkladu). Je preto potrebné
tieto hodnoty najskôr zrušiť napr. príkazom kill(K1,K2) alebo 
kill(all) zrušiť všetky už definované hodnoty.

Y1: coeff(rhs(VR_HOMOG), %k1)$
Y2: coeff(rhs(VR_HOMOG), %k2)$
depends(K1, x)$
depends(K2, x)$
PODM1: diff(K1,x)*Y1 + diff(K2,x)*Y2 = 0;
PODM2: diff(K1,x)*diff(Y1,x) + diff(K2,x)*diff(Y2,x)
= rhs(DR);

Riešením systému rovníc PODM1 a PODM2 získame derivácie funkcií K1 a K2

RIES: solve([PODM1,PODM2], [diff(K1,x),diff(K2,x)]);

Integráciou výsledkov (pozor - riešenie systému rovníc má
dátovú štruktúru "zoznam zoznamov", dôležité pre vyberanie
jednotlivých položiek) získame hľadané funkcie K1 a K2.
S použitím K1 a K2 zostrojíme partikulárne riešenie YP.

K1: integrate(rhs(RIES[1][1]),x)$
K2: integrate(rhs(RIES[1][2]),x)$
YP: ratsimp(K1*Y1 + K2*Y2);

Všeobecné riešenie pôvodnej nehomogénnej rovnice teraz získame
ako súčet všeobecného a partikulárneho riešenia.

VSEOB_R: lhs(VR_HOMOG) = rhs(VR_HOMOG) + YP;

kill(all);

### Špeciálny tvar pravej strany

#### Polynóm

Vyriešime rovnicu

y''-4*y'+4*y=x^2

Definujeme závislosť y od x a tvar diferenciálnej rovnice

depends(y, x)$
DR: diff(y,x,2)-4*diff(y,x)+4*y= x^2;

Zo zadania vytvoríme homogenizovanú rovnicu DR_HOMOG s nulovou pravou
stranou. Jej všeobecné riešenie VR_HOMOG získame pomocou príkazu ode2().
Pri rovniciach rádu vyššieho než 2 opakuejme ručný postup z časti
o homogénnej rovnici.

DR_HOMOG: lhs(DR) = 0$
VR_HOMOG: ode2(DR_HOMOG, y, x);

Na pravej strane pôvodnej rovnice je polynóm

rhs(DR-DR_HOMOG);

Zo všeobecného riešenia VR_HOMOG vidíme, že 0 nie je koreňom
charakteristickej rovnice. Preto partikulárne riešenie predpokladáme
v tvare polynómu rovnakého stupňa s neznámymi koeficientami a,b,c

YP:a*x^2+b*x+c;

Pripravíme si derivácie partikulárneho riešenia,
aby sme ich mohli dosadiť do ľavej strany diferenciálnej
rovnice

SUB: diff(YP,x,2) = diff(y,x,2)$
SUB1:diff(YP,x,1)=diff(y,x,1)$
SUB2: YP=y$

Vykonáme pripravené substitúcie - derivácie partikulárneho riešenia

subst(lhs(SUB), rhs(SUB), DR)$
subst(lhs(SUB1), rhs(SUB1), %)$
POM_ROV:expand(subst(lhs(SUB2), rhs(SUB2), %));

Z poslednej rovnice POM_ROV získame porovnaním koeficientov
pri rovnakých mocninách systém rovníc pre koeficienty a,b,c.
Riešením tohto systému získame konkrétne hodnoty

ROV1:coeff(lhs(POM_ROV), x^2)=coeff(rhs(DR), x^2)$
ROV2:coeff(lhs(POM_ROV), x)=coeff(rhs(DR), x)$
ROV3:subst(0,x,lhs(POM_ROV))=subst(0,x,rhs(DR))$
KOEF:solve([ROV1,ROV2,ROV3],[a,b,c]);

S využitím získaných koeficientov zostavíme partikulárne riešenie,
ktoré pripočítame ku riešeniu všeobecnému.

YP: rhs(KOEF[1][1])*x^2+ rhs(KOEF[1][2])*x+ rhs(KOEF[1][3])$
VSEOB_R:lhs(VR_HOMOG) =rhs(VR_HOMOG) + YP;

#### Pravá strana v tvare P(x)*exp(αx), P(x) je polynóm

Vyriešime rovnicu

y''+y'-2*y=x*e^(-2x)

Definujeme závislosť y od x a tvar diferenciálnej rovnice

depends(y, x)$
DR: diff(y,x,2)+diff(y,x)-2*y= x*exp(-2*x);

Zo zadania vytvoríme homogenizovanú rovnicu DR_HOMOG s nulovou pravou
stranou. Jej všeobecné riešenie VR_HOMOG získame pomocou príkazu ode2().
Pri rovniciach rádu vyššieho než 2 opakuejme ručný postup z časti
o homogénnej rovnici.

DR_HOMOG: lhs(DR) = 0$
VR_HOMOG: ode2(DR_HOMOG, y, x);

Overíme tvar pravej strany pôvodnej rovnice

rhs(DR-DR_HOMOG);

Zo všeobecného riešenia VR_HOMOG vidíme, že -2 je jednoduchým koreňom
charakteristickej rovnice. Preto partikulárne riešenie predpokladáme
v tvare polynómu prvého stupňa s neznámymi koeficientami a,b vynásobeného x
a exp(-2x)

YP:x*(a*x+b)*exp(-2*x);

Pripravíme si substitúcie - derivácie partikulárneho riešenia

SUB: diff(YP,x,2) = diff(y,x,2)$
SUB1:diff(YP,x,1)=diff(y,x,1)$
SUB2: YP=y$

A teraz substitúcie vykonáme

subst(lhs(SUB), rhs(SUB), DR)$
subst(lhs(SUB1), rhs(SUB1), %)$
POM_ROV:expand(subst(lhs(SUB2), rhs(SUB2), %)*exp(2*x));

Aby sme mohli porovnávať pri rovnakých mocninách, vynásobili sme rovnicu
exp(2*x), čím sme sa zbavili členov v tavre exp(-2*x), ktoré znemožňujú
vyberanie koeficientov (resp. koeficienty sú vytvárané aj s mocninami x).

Z poslednej rovnice POM_ROV získame porovnaním koeficientov
pri rovnakých mocninách systém rovníc pre koeficienty a,b.

Pre získanie absolútnych členov rovnice (teda koeficientov pri x^0), sme
využili dosadenie x=0.

Riešením tohto systému získame konkrétne hodnoty

ROV1:coeff(lhs(POM_ROV), x)=coeff(rhs(POM_ROV),x)$
ROV2:subst(0,x,lhs(POM_ROV))=subst(0,x,rhs(POM_ROV))$
KOEF:solve([ROV1,ROV2],[a,b]);

S využitím získaných koeficientov zostavíme partikulárne riešenie,
ktoré pripočítame ku riešeniu všeobecnému.

YP: rhs(KOEF[1][1])*x^2*exp(-2*x) + rhs(KOEF[1][2])*x*exp(-2*x)$
VSEOB_R:lhs(VR_HOMOG) =rhs(VR_HOMOG) + YP;

#### Pravá strana v tvare  P(x)*exp(αx)*cos(βx) resp-
P(x)*exp(αx)*sin(βx), P(x) je polynóm

Vyriešime rovnicu

y''-y=x*cos(x)*e^(x)

Definujeme závislosť y od x a tvar diferenciálnej rovnice

depends(y, x)$
DR: diff(y,x,2)-y= x*cos(x)*exp(x);

Zo zadania vytvoríme homogenizovanú rovnicu DR_HOMOG s nulovou pravou
stranou. Jej všeobecné riešenie VR_HOMOG získame pomocou príkazu ode2().
Pri rovniciach rádu vyššieho než 2 opakuejme ručný postup z časti
o homogénnej rovnici.

DR_HOMOG: lhs(DR) = 0$
VR_HOMOG: ode2(DR_HOMOG, y, x);

Overíme tvar pravej strany pôvodnej rovnice

rhs(DR-DR_HOMOG);

Zo všeobecného riešenia VR_HOMOG vidíme, že 1+i nie je koreňom
charakteristickej rovnice. Preto partikulárne riešenie predpokladáme
v tvare 

yp=(ax+b)cos(x)e^x + (cx+d)sin(x)e^x

POZN.: Všimnite si, že do partikulárneho riešenia vstupujú obe zložky,
 tj. komponenty obsahujúce cos(x) ak sin(x).

YP:(a*x+b)*cos(x)*exp(x)+(c*x+d)*sin(x)*exp(x);

Vykonáme pripravené substitúcie - derivácie partikulárneho riešenia.

Pre účely porovnania koeficientov vynásobíme výslednú rovnicu
e^(-x), čím sa zbavíme exponenciálnej funkcie.

SUB: diff(YP,x,2) = diff(y,x,2)$
SUB1:diff(YP,x,1)=diff(y,x,1)$
SUB2: YP=y$
subst(lhs(SUB), rhs(SUB), DR)$
subst(lhs(SUB1), rhs(SUB1), %)$
POM_ROV:expand(subst(lhs(SUB2), rhs(SUB2), %)*exp(-x));

Aby sme mohli porovnávať koeficienty podľa mocnín x, je ešte potrebné
z rovnice POM_ROV "odstrániť" funkcie sin(x) resp. cos(x). To môžeme
zrealizovať tak, že postupne za sin(x) dosadíme 1 a za cos(x) nulu, 
resp. naopak za sin(x) nulu a za cos(x) číslo 1.
Získame tak dve pomocné rovnice, v ktorých už budeme porovnávať
koeficienty pri mocninách x.

POM_ROV_sin:subst(1,sin(x),subst(0,cos(x),POM_ROV));
POM_ROV_cos:subst(1,cos(x),subst(0,sin(x),POM_ROV));

Zo získaných rovníc POM_ROV_sin a POM_ROV_cos zostavíme porovnaním koeficientov
pri rovnakých mocninách systém rovníc pre hľadané koeficienty a,b,c,d.

Pre získanie absolútnych členov rovnice (teda koeficientov pri x^0), sme
využili dosadenie x=0.

Riešením tohto systému získame konkrétne hodnoty

ROV1:coeff(lhs(POM_ROV_sin), x)=coeff(rhs(POM_ROV_sin),x)$
ROV2:subst(0,x,lhs(POM_ROV_sin))=subst(0,x,rhs(POM_ROV_sin))$
ROV3:coeff(lhs(POM_ROV_cos), x)=coeff(rhs(POM_ROV_cos),x)$
ROV4:subst(0,x,lhs(POM_ROV_cos))=subst(0,x,rhs(POM_ROV_cos))$
KOEF:solve([ROV1,ROV2,ROV3,ROV4],[a,b,c,d]);

S využitím získaných koeficientov zostavíme partikulárne riešenie,
ktoré pripočítame ku riešeniu všeobecnému.

YP: rhs((KOEF[1][1])*x+rhs(KOEF[1][2]))*cos(x)*exp(x) + (rhs(KOEF[1][3])*x+rhs(KOEF[1][4]))*sin(x)*exp(x)$
VSEOB_R:lhs(VR_HOMOG) =rhs(VR_HOMOG) + YP;