load(draw)$

# Numeričko traženje ekstrema

Traženje ekstrema s teorijske strane je riješen problem. Međutim, u praksi se mogu javiti komplicirani sustavi koje nije moguće egzaktno riješiti kako bismo dobili stacionarne točke. U tom slučaju moramo upotrijebiti neku numeričku metodu. Ovdje nećemo ulaziti u detalje tih metoda niti ćemo ih spominjati jer to prelazi okvire ovog predmeta. Na jednom jednostavnom primjeru pokazat ćemo situaciju u kojoj možemo naići na određene poteškoće koje su više vezane uz samu Maximu, iako se ručnim rješavanjem ekstremi mogu egzaktno pronaći.

## Zadatak

Odredite ekstreme funkcije f(x,y)=4x+3y+7 uz uvjet x^4.24*y^3.91=18.

### Rješenje (1. način)

Iz uvjeta možemo izraziti jednu varijablu pomoću druge pa problem svodimo na traženje ekstrema funkcije jedne varijable

f(x,y):=4*x+3*y+7;

Kako je opća potencija (koja se javlja u uvjetu) definirana za pozitivne brojeve, možemo pretpostaviti da su varijable x i y pozitivne.

assume(x>0,y>0);

Maxima zna iz zadanog uvjeta jednu varijablu izraziti pomoću druge, no problem je što ona jednadžbu rješava u skupu kompleksnih brojeva pa na izlazu daje 391 moguće rješenje. Rješenja nisu namjerno ispisana zbog uštede prostora (zato je stavljen znak dolara na kraju naredbe).

rj:solve(x^4.24*y^3.91=18,y)$

length(rj);

U ovom slučaju ovaj problem možemo lagano izbjeći tako da jednadžbu rješavamo numeričkim putem i u tom slučaju dobivamo samo jedno realno rješenje kao što bismo dobili i ručnim rješavanjem jednadžbe. 

rj:solve(x^4.24*y^3.91=18,y),numer;

rj[1];

rhs(rj[1]);

y0:rhs(rj[1]);

Dobiveno rješenje uvrstimo u početnu funkciju f i dobivamo funkciju jedne varijable.

fun:f(x,y0);

Tražimo stacionarne točke dobivene funkcije jedne varijable. Jednadžbu rješavamo ponovo numeričkim putem kako bismo izbjegli ranije navedene probleme.

der:diff(fun,x);

stac:solve(der=0,x),numer;

x0:rhs(stac[1]);

U slučaju kompliciranijih jednadžbi možemo koristiti naredbu find_root koja zahtijeva da pogodimo neki dovoljno mali interval koji sadrži nultočku. U većini situacija to možemo pogoditi pomoću slike tako da funkciju crtamo unutar nekog intervala tako dugo dok graf ne presječe x-os.

wxplot2d(der,[x,0.8,5]);

Vidimo sa slike da je nultočka sigurno npr. unutar segmenta [1,3]

x0:find_root(der,1,3);

Kako je druga derivacija u stacionarnoj točki pozitivna, u pripadnoj stacionarnoj točki x0=1.29111035643873 funkcija postiže minimum koji je u ovom slučaju ujedno i globalni.

der2:diff(fun,x,2);

subst(x0,x,der2);

Iz uvjeta y0 izračunamo vrijednost druge varijable.

y0:subst(x0,x,y0);

Funkcija f(x,y)=4x+3y+7 uz uvjet x^4.24*y^3.91=18 postiže uvjetni globalni minimum u točki (x0,y0)=(1.29111035643873,1.587497325055168) i taj minimum iznosi 16.92693340092043.

f(x0,y0);

### Rješenje (2. način)

pomoću Lagrangeove funkcije

L(x,y,t):=f(x,y)+t*(x^4.24*y^3.91-18);

L(x,y,t);

Lx:diff(L(x,y,t),x);

Ly:diff(L(x,y,t),y);

Lt:diff(L(x,y,t),t);

Maxima ne zna riješiti dobiveni sustav, iako se on može ručno lagano riješiti. Stoga trebamo maximi malo pomoći.

solve([Lx=0,Ly=0,Lt=0],[x,y,t]);

Iz prve dvije jednadžbe možemo lagano eliminirati varijablu t (dodatni parametar koji smo uveli u Lagrangeovoj funkciji).

el:eliminate([Lx=0,Ly=0],[t]);

el[1];

Iz treće jednadžbe možemo dobiti eksplicitnu vezu između nepoznanica x i y. Nakon toga tu vezu uvrstimo u gornju jednadžbu el[1].

rj:solve(Lt=0,y),numer;

rj[1];

Dobivamo jednadžbu s jednom nepoznanicom koju rješavamo numeričkim putem.

jed:subst(rj[1],el[1]);

Međutim solve naredba numeričkim putem ne daje rješenje jer joj je zapis jednadžbe u nezgodnom obliku.

solve(jed=0,x),numer;

Možemo jednadžbu pomnožiti s x pa dobivamo novu ekvivalentnu jednadžbu koju maxima zna riješiti.

jed2:expand(jed*x);

x0:solve(jed2=0,x),numer;

x0[1];

x0:rhs(x0[1]);

Ili možemo na početnoj jednadžbi koristiti naredbu find_root kojoj ne smeta početni zapis. Naravno, opet moramo pomoću slike otkriti neki interval unutar kojeg se nalazi nultočka pripadne funkcije.

wxplot2d(jed,[x,0.8,5]);

x0:find_root(jed,1,3);

Izračunamo drugu varijablu iz uvjeta rj[1].

rj[1];

y0:subst(x0,x,rhs(rj[1]));

Ponovo smo dobili istu stacionarnu točku x0=1.29111035643873, y0=1.587497325055168. No, u ovom slučaju ne znamo postiže li se u toj točki uopće uvjetni ekstrem i ako se postiže je li minimum ili maksimum. Uvjet x^4.24*y^3.91=18 nije kompaktni skup u ravnini, što možemo vidjeti na slici, a iz same jednadžbe je jasno da pripadna krivulja nije omeđena.

wxdraw2d(grid=true,color=red,line_width=2,xlabel="x",ylabel="y",
    ip_grid_in=[10,10],xrange=[0,5],yrange=[0,5],
implicit(x^4.24*y^3.91=18,x,0.4,5,y,0.3,5));

Stoga moramo provesti daljnja ispitivanja Lagrangeove funkcije. U tom slučaju trebamo i pripadnu vrijednost parametra t, kojeg možemo dobiti iz prve ili druge jednadžbe sustava jer vrijednosti od x i y sada znamo.

Lx;

jed3:subst([x=x0,y=y0],Lx);

Vrijednost parametra t

t0:solve(jed3=0,t);

t0:rhs(t0[1]);

Sada odredimo sve parcijalne derivacije drugog reda Lagrangeove funkcije.

hes:hessian(L(x,y,t),[x,y,t]);

Izračunamo Hesseovu determinantu u dobivenoj stacionarnoj točki (x0,y0,t0). Ako je ta determinanta strogo manja od nula, tada se radi o lokalnom minimumu, a ako je strogo veća od nule, tada se radi o lokalnom maksimumu.

hes0:subst([x=x0,y=y0,t=t0],hes);

Funkcija f(x,y)=4x+3y+7 uz uvjet x^4.24*y^3.91=18 postiže uvjetni lokalni minimum u točki (x0,y0)=(1.29111035643873,1.587497325055168) i taj minimum iznosi 16.92693340092043.
U ovom slučaju taj uvjetni minimum je ujedno i globalni ako pogledamo donju 3D sliku.

determinant(hes0);

f(x0,y0);

## 3D slika

wxdraw3d(xu_grid=30,yv_grid=30,surface_hide=false,
xlabel="x",ylabel="y",zlabel="f(x,y)",view=[70,30],
color=red,line_width=2,parametric(u,(18/u^4.24)^(1/3.91),0,u,0.3,5),
point_type=7,point_size=1.2,points([[x0,y0,0]]),
color=skyblue,line_width=1,explicit(f(x,y),x,0,5,y,0,7.7),
color=blue,line_width=2,parametric(u,(18/u^4.24)^(1/3.91),f(u,(18/u^4.24)^(1/3.91)),u,0.3,5),
points([[x0,y0,f(x0,y0)]]),user_preamble= "set xyplane at 0"),wxplot_size=[800,600];