| (%i1) | 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.
1 Zadatak
Odredite ekstreme funkcije f(x,y)=4x+3y+7 uz uvjet x^4.24*y^3.91=18.
1.1 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
| (%i2) | 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.
| (%i3) | assume(x>0,y>0); |
![(%o3) [x>0,y>0]](ekstremi_numerika_htmlimg/ekstremi_numerika_1.png)
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).
| (%i4) | rj:solve(x^4.24·y^3.91=18,y)$ |
| (%i5) | 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.
| (%i6) | rj:solve(x^4.24·y^3.91=18,y),numer; |
![(rj) [y=2.094312869232622/x^1.084398976982097]](ekstremi_numerika_htmlimg/ekstremi_numerika_5.png)
| (%i7) | rj[1]; |
| (%i8) | rhs(rj[1]); |
| (%i9) | y0:rhs(rj[1]); |
Dobiveno rješenje uvrstimo u početnu funkciju f i dobivamo funkciju jedne varijable.
| (%i10) | 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.
| (%i11) | der:diff(fun,x); |
| (%i12) | stac:solve(der=0,x),numer; |
![(stac) [x=1.291110356438728]](ekstremi_numerika_htmlimg/ekstremi_numerika_12.png)
| (%i13) | 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.
| (%i14) | wxplot2d(der,[x,0.8,5]); |
Vidimo sa slike da je nultočka sigurno npr. unutar segmenta [1,3]
| (%i15) | 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.
| (%i16) | der2:diff(fun,x,2); |
| (%i17) | subst(x0,x,der2); |
Iz uvjeta y0 izračunamo vrijednost druge varijable.
| (%i18) | 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.
| (%i19) | f(x0,y0); |
1.2 Rješenje (2. način)
pomoću Lagrangeove funkcije
| (%i20) | L(x,y,t):=f(x,y)+t·(x^4.24·y^3.91−18); |
| (%i21) | L(x,y,t); |
| (%i22) | Lx:diff(L(x,y,t),x); |
| (%i23) | Ly:diff(L(x,y,t),y); |
| (%i24) | 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.
| (%i25) | solve([Lx=0,Ly=0,Lt=0],[x,y,t]); |
![(%o25) []](ekstremi_numerika_htmlimg/ekstremi_numerika_28.png)
Iz prve dvije jednadžbe možemo lagano eliminirati varijablu t (dodatni parametar koji smo uveli u Lagrangeovoj funkciji).
| (%i26) | el:eliminate([Lx=0,Ly=0],[t]); |
![(el) [100*(318*x^3.24*y^3.91-391*x^4.24*y^2.91)]](ekstremi_numerika_htmlimg/ekstremi_numerika_30.png)
| (%i27) | 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].
| (%i28) | rj:solve(Lt=0,y),numer; |
![(rj) [y=2.094312869232622/x^1.084398976982097]](ekstremi_numerika_htmlimg/ekstremi_numerika_33.png)
| (%i29) | rj[1]; |
Dobivamo jednadžbu s jednom nepoznanicom koju rješavamo numeričkim putem.
| (%i30) | 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.
| (%i31) | solve(jed=0,x),numer; |
![(%o31) [x^1.084398976982097=1.703303049657219/x]](ekstremi_numerika_htmlimg/ekstremi_numerika_37.png)
Možemo jednadžbu pomnožiti s x pa dobivamo novu ekvivalentnu jednadžbu koju maxima zna riješiti.
| (%i32) | jed2:expand(jed·x); |
| (%i33) | x0:solve(jed2=0,x),numer; |
![(x0) [x=1.29111035643873]](ekstremi_numerika_htmlimg/ekstremi_numerika_40.png)
| (%i34) | x0[1]; |
| (%i35) | 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.
| (%i36) | wxplot2d(jed,[x,0.8,5]); |
| (%i37) | x0:find_root(jed,1,3); |
Izračunamo drugu varijablu iz uvjeta rj[1].
| (%i38) | rj[1]; |
| (%i39) | 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.
| (%i40) |
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.
| (%i41) | Lx; |
| (%i42) | jed3:subst([x=x0,y=y0],Lx); |
Vrijednost parametra t
| (%i43) | t0:solve(jed3=0,t); |
![(t0) [t=-10997204/162516435]](ekstremi_numerika_htmlimg/ekstremi_numerika_56.png)
| (%i44) | t0:rhs(t0[1]); |
Sada odredimo sve parcijalne derivacije drugog reda Lagrangeove funkcije.
| (%i45) | hes:hessian(L(x,y,t),[x,y,t]); |
![(hes) matrix(<BR>
[13.7376*t*x^2.24*y^3.91, 16.5784*t*x^3.24*y^2.91, 4.24*x^3.24*y^3.91],<BR>
[16.5784*t*x^3.24*y^2.91, 11.3781*t*x^4.24*y^1.91, 3.91*x^4.24*y^2.91],<BR>
[4.24*x^3.24*y^3.91, 3.91*x^4.24*y^2.91, 0]<BR>
)](ekstremi_numerika_htmlimg/ekstremi_numerika_58.png)
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.
| (%i46) | hes0:subst([x=x0,y=y0,t=t0],hes); |
![(hes0) matrix(<BR>
[-10.03787161598452, -9.851985104577398, 59.11191062746484],<BR>
[-9.851985104577398, -5.499221864639432, 44.33393297059865],<BR>
[59.11191062746484, 44.33393297059865, 0]<BR>
)](ekstremi_numerika_htmlimg/ekstremi_numerika_59.png)
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.
| (%i47) | determinant(hes0); |
| (%i48) | f(x0,y0); |
2 3D slika
| (%i49) |
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]; |
