## Smerové pole diferenciálnej rovnice

Ak si pripomenieme, že geometrický význam derivácie funkcie je smernica
dotyčnice v danom bode. Ak teda máme rovnicu v tvare y'=f(x,y),
tak po dosadní súradníc nejakého prípustného bodu v rovine
do pravej strany rovnice dostaneme ako výsledok smernicu dotyčnice
ku grafu funkcie, ktorá je riešením danej rovnice a prechádza týmto
konkrétnym bodom. Ak v každom bode vyznačíme túto dotyčnicu krátkym
vektorom v smere dotyčnice, dostaneme smeerové pole, ktoré poskytuje
určitú geometrickú predstavu o riešeniach rovnice.

V systéme wxmaxima zobrazíme smerové pole pomocou funkcie plotdf().
(V závislosti od verzie maximy môže byť potrebné načátať balíček plotdf
príkazom load(plotdf).

Zobrazme smerové pole rovnice y'=x*y.

plotdf(x*y);

Ako argument funkcie plotdf() sme dosadili pravú stranu
rovnice. Ak klikneme myšou do niektorej pozície v smerovom poli,
zobrazí sa integrálna krivka (graf riešenia), ktorá prechádza daným bodom.

Všimnime si, že maxima nevykresľuje všetky vektory smerového poľa
rovnakou dĺžkou. Tá sa zväčšuje v závislosti na absolútnej hodnote
smernice.

Je možné zobraziť v smerovom poli integrálnu krivku,
ktorá prechádza zadaným bodom (a zodpovedá napr. riešeniu
začiatočnej úlohy y'=f(x,y), y(x0)=y0.
    
Zobrazme smerové pole rovnice y'=x*y  spolu s trajektóriou
riešenia spĺňajúceho y(1)=1. 

plotdf(x*y,[trajectory_at,1,1]);

Ďalšou voľbou môžeme rozhodnúť, či má daná trajektória vo zvolenom bode
začínať (forward), končiť (backward) alebo prechádzať (both).
Ukážky

plotdf(x*y,[trajectory_at,1,1],[direction,forward]);

plotdf(x*y,[trajectory_at,1,1],[direction,backward]);

Alternatívne môžeme smerové pole zobraziž pomocou drawdf(), resp wxdrawdf(),
musíme však načítať balíček drawdf.

load(drawdf);

drawdf(x*y,[x,-5,5],[y,-5,5]);

Tu však nie je možné kliknutím zobraziť trajektórie. Tie je potrebné
explicitne zadať ako argument funkcie.

drawdf(x*y,[x,-5,5],[y,-5,5],[trajectory_at,1,1]);

Môžeme zadať aj viac trajektórií do jedného obrázku, napr.

drawdf(x*y,[x,-5,5],[y,-5,5],[trajectory_at,1,1],[trajectory_at,-1,0.5],[trajectory_at,-2,-1] );

## Rovnice so separovanými premennými

### Ručné riešenie

logabs:true;

Riešme rovnicu y-y^2+xy'=0.

Najskôr musíme deklarovať skutočnosť, že y je funkciou nezávisle premennej x:

depends(y, x)$

Definujeme rovnicu

DR: y - y^2 + x*diff(y,x) = 0;

Riešením vzhľadom na y' osamostatníme deriváciu na jednej starne rovnice

DR2: solve(DR, diff(y,x))[1];

V zhode s označením v slajdoch (strana č. 10) z predchádzajúcej rovnice
separujeme premenné a dostaneme riešenie

q: 1/denom(rhs(DR2));
p: 1/num(rhs(DR2));
integrate(q, x) = integrate(p, y) + %c;

Pomocou Simplify->Contract Logarithms v menu zjednodušíme:

logcontract(%);

Ak si uvedomíme, že integračná konštantu na pravej strane
môžeme považovať za logaritmus nejakej inej (kladnej) konštanty k,
máme na pravje strane súčet dvoch logaritmov, čo je vlastne
logaritmus súčinu, teda log(x) = log(k*(y-1)/y. Ak rovnicu vyriešime
vzhľadom na y, dostávame riešenie:

solve(x=k*(y-1)/y,y);

Poznámka: Možno ste si všimli, že vo výsledku integrácie 1/x nám chýba
absolútna hodnota. jej zobrazenie je možné docieliť pomocou nastavenia
premennej logabs na hodnotu true,teda

logabs:true;

Potom dostávame riešenie v takomto tvare

q: 1/denom(rhs(DR2));
p: 1/num(rhs(DR2));
integrate(q, x) = integrate(p, y) + %c;

logcontract(%);

Opäť, integračná konštanta vpravo je logaritmom nejakej inej,
kladnej konštanty k a súčet logaritmov je logaritmus súčinu.
Ak si uvedomíme, že absolútna hodnota má za následok nanajvýš
zmenu znamienka, môžeme rozšíriť obor konštanty k aj na záporné
hodnoty a absolútne hodnoty tak odstrániť. 
Dostaneme rovnaký výsledok

solve(x=k*(y-1)/y,y);

Zobrazovanie absolútnej hodnoty vo výsledkoch integrácie vypneme
nastavením logabs na hodnotu false:

logabs:false;

### Riešenie pomocou ode2

Rovnaký výsledok získame aj použitím príkazu ode2 na pôvodnú rovnicu

ode2(DR, y, x);

logcontract(%);

exp(lhs(%)) = exp(rhs(%));

exp(lhs(%)) = exp(rhs(%));

VR:solve(%,y)[1];

Príkazom method sa dozvieme, že maxima riešila rovnicu metódou
separácie premenných.

method;

Získali sme tak všeobecné riešenie rovnice. Pokiaľ by sme chceli riešiť
začiatočnú úlohu, použijeme príkaz ic1(), ktorého argumentmi sú všeobecné riešenie
a hodnoty x=x0 a y=y0 dané začiatočnou podmienkou.
Nájdime riešenie rovnice, ktoré spĺňa začiatočnú podmienku y(2)=3.

ic1(VR,x=2,y=3);

## Homogénna rovnica

Homogénna rovnica je rovnica, ktorú môžeme písať v tvare

y'=f(y/x).

Inými slovami, vieme ju vždy upraviť na tvar,
v ktorom závisle premenná a nezávisle premenná vystupujú
len v podobe vzájomného pomeru y/x.

### Ručné riešenie

Substitúciou u=y/x, odkiaľ y=u*x a derivovaním (pozor na deriváciu súčinu)
dostaneme y'=u'*x+u. Tak rovnicu prevedieme na rovnicu so separovanými
premennými.

Riešme rovnicu y'=-(x+y)/x

depends(y,x);

DR:diff(y,x,1)=-(x+y)/x;

Osamostatníme na ľavej strane y'

DR2: expand(solve(DR, diff(y,x))[1]);

Vykonáme substitúciu u=y/x

depends(u,x)$
SUB:u=y/x;

solve(SUB,y)[1]$
SUB2:diff(%,x);

Nasleduje samotná substitúcia do pôvodnej rovnice a jej úprava

subst(rhs(SUB2), lhs(SUB2), DR2)$
subst(lhs(SUB), rhs(SUB), %)$
DR_SEP: solve(%, diff(u,x))[1];

Dostali sme rovnicu so separovanými premennými, môžeme ju vyriešiť
pomocou ode2.

ode2(DR_SEP,u,x);

Vykonáme spätnú substitúciu a vyriešime vzhľadom na y

subst(rhs(SUB), lhs(SUB), %);
solve(%, y)[1];

Rovnaký výsledok (odlišuje sa len opticky) získame aj použitím
príkazu ode2 priamo na pôvodnú rovnicu.

ode2(DR, y, x);

Riešme rovnicu 2x^3*y'=y(2x^2-y^2)

DR:2*x^3*diff(y,x,1)=y*(2*x^2-y^2);

Osamostatníme na ľavej strane y'

DR2: expand(solve(DR, diff(y,x))[1]);

Vykonáme substitúciu u=y/x (bola definovaná vyššie,
tu už len skutočne dosadíme do našej rovnice)

subst(rhs(SUB2), lhs(SUB2), DR2)$
subst(lhs(SUB), rhs(SUB), %)$
DR_SEP: solve(%, diff(u,x))[1];

Pravú stranu musíme upraviť ručne. Vidíme, že na pravej
strane je vlastne -1/2x*(y/x)^3=-1/2x*u^3

DR_SEP:subst(-u^3/(2*x),rhs(%),DR_SEP);

Riešime rovnicu so separovanými premennými

ode2(DR_SEP,u,x);

Vykonáme spätnú substitúciu

subst(rhs(SUB), lhs(SUB), %);
solve(%, y)[1];

ode2(DR, y, x);

method;

## Lineárna dif. rovnica prvého rádu

Ide o rovnicu v tvare

y'=a(x)*y+b(x)

O jej riešení bližšie slajdy, str. 20
Tu si ilustrujeme postup pri priešení metódou variácie konštanty.

Riešme rovnicu

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

Zadefinujeme si rovnicu

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

Osamostatníme y'

DR2: expand(solve(DR, diff(y,x))[1]);

Ako prvý krok riešime homogénnu rovnicu bez b(x) na pravej strane.
Z predchádzajúcej rovnice teda vyberieme prvé dva členy na pravej strane,
ktoré obsahujú y. Využijeme na to príkaz part()

HOM: lhs(DR2) = part(rhs(DR2),1) + part(rhs(DR2),2);

Nájdeme riešenie homogenizovanej rovnice

HOM_RIES: ode2(HOM, y, x);

Variácia konštanty, predpokladáme, že c je funkciou závislou od x
a hľadáme jej vhodný tvar.
Definujeme závislosť C na x a dosadíme do riešenia homogenizovanej rovnice

depends(C, x)$
Y1: subst(C, %c, HOM_RIES);

Do rovnice označenej DR2 dosadíme riešenie v tvare Y1 a riešime pre
neznámu funkciu C premennej x:

DR2, diff(Y1, x), Y1$
solve(%, diff(C, x))[1]$
C_VYS: C = rhs(integrate(%, x));

Na záver dosadíme funkciu C (uloženú v C_VYS) do riešenia Y1
homogenizovanej rovnice HOM

Y1,C_VYS;

Opäť, pri priamom použití príkazu ode2() dostávame rovnaký výsledok,
príkaz method nám prezradí, že bolo nájdené ako riešenie lineárnej rovnice

ode2(DR, y, x);
method;

## Bernoulliho rovnica

Bernoulliho rovnicou rádu r rozumieme rovnicu v tvare

y' = a(x)y + b(x)y^r ; r≠0 a r≠1

Substitúciou u=y^(1-r) ju prevedieme na lineárnu rovnicu, ktorú
vyriešime metódou variácioe konštanty.

Vyriešme rovnicu y'+y*x=x^3*y^2

Definujeme y ako funkciu premennej x a zadefinujeme zadanú rovnicu

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

Vykonáme riešiacu substitúciu u=y^(1-r), v našom prípade je r=2

depends(u, x)$
SUB: u = y^(1-2);
SUB2: diff(SUB, x);

Výsledky vyplývajúce zo substitúcie dosadíme do pôvodnej rovnice DR

expand(solve(DR, diff(y,x))[1]/y^2)$
subst(lhs(SUB), rhs(SUB), %)$
DR2: ratsubst(lhs(SUB2), rhs(SUB2), %);

Získali sme lineárnu rovnicu, ktorú už vieme riešiť

ode2(DR2, u, x), SUB$
VR: solve(%, y)[1];
method;

Rovnaký výsledok získame priamym použitím príkazu ode2() na pôvodnú
rovnicu, V premennej odeindex je uložená hodnota rádu r Bernoulliho rovnice.

expand(ode2(DR, y, x));
method;
odeindex;