Harj. 8 LV

22.3.2002 HA ja JP

Alustukset

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra):with(linalg):with(plottools): with(plots):setoptions3d(axes=boxed,orientation=[-30,50]):

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

> V2L:=vek->convert(vek,list):

1.

Tämä on puhtaasti vanhan kertausta. Tämä nyt lähinnä siksi, ettei kaikki ole kahden muuttujan tapauksia. Samalla harjaannutaan Maple-työskentelyssä tutun teeman parissa.

> f:=(x,y,z)->4*x*y*z-x^4-y^4-z^4;

f := proc (x, y, z) options operator, arrow; 4*x*y*...

> g:=grad(f(x,y,z),[x,y,z]);

g := vector([4*y*z-4*x^3, 4*x*z-4*y^3, 4*x*y-4*z^3]...

> solve({g[1]=0,g[2]=0,g[3]=0},{x,y,z});

{y = 0, z = 0, x = 0}, {x = 1, y = 1, z = 1}, {x = ...
{y = 0, z = 0, x = 0}, {x = 1, y = 1, z = 1}, {x = ...
{y = 0, z = 0, x = 0}, {x = 1, y = 1, z = 1}, {x = ...
{y = 0, z = 0, x = 0}, {x = 1, y = 1, z = 1}, {x = ...

> H:=Matrix(hessian(f(x,y,z),[x,y,z]));

H := _rtable[136382060]

> H111:=subs(x=1,y=1,z=1,H);

H111 := _rtable[136431448]

> Eigenvalues(H111);

_rtable[136492332]

Ominaisarvot ovat kaikki negatiivisia, joten kyseessä on negatiivisesti definiitti matriisi ja (1,1,1) on maksimipiste.

Huom! Kirjallisuudessa esiintyy diagonaalialideterminanttien avulla lausuttavia ehtoja.

Niitä ei tarvitse osata. (Jos joku soveltaa niitä oikein, niin ok, mutta en kehoita opettelemaan.)

Ominaisarvojen merkit kertovat syyn asioiden tilaan, siksi niitä suosin vahvasti.

Koetehtävissä ei tule olemaan mitään laskennallista etua determinanttiehtojen osaamisesta.

2.

> q:=3*x^2+2*x*y + 3*y^2;

q := 3*x^2+2*x*y+3*y^2

> A:=<<3,1>|<1,3>>;

A := _rtable[135111092]

> det(A);

8

Ominaisarvot ovat samanmerkkiset, merkki on sama kuin päälävistäjän alkioiden (summan) merkki. Siispä yhtälö on muotoa

> lambda[1]*y[1]^2+lambda[2]*y[2]^2=16;

lambda[1]*y[1]^2+lambda[2]*y[2]^2 = 16

missä lambda[1] ja lambda[2] ovat positiivisia, joten ellipsihän siitä tulee.

b)

> g:=q-16; # rajoitefunktio

g := 3*x^2+2*x*y+3*y^2-16

> f:=x^2+y^2; # kohdefunktio

f := x^2+y^2

> L:=f+lambda*g; # Lagrangen funktio

L := x^2+y^2+lambda*(3*x^2+2*x*y+3*y^2-16)

> lagyht:={diff(L,x)=0,diff(L,y)=0,g=0};

lagyht := {3*x^2+2*x*y+3*y^2-16 = 0, 2*x+lambda*(6*...

> solve(lagyht,{x,y,lambda});

{x = RootOf(_Z^2-2), lambda = -1/4, y = RootOf(_Z^2...

> ratk:=map(allvalues,[%]);

ratk := [{y = sqrt(2), lambda = -1/4, x = sqrt(2)},...

> p1:=subs(ratk[1],[x,y]);p2:=subs(ratk[2],[x,y]);p3:=subs(ratk[3],[x,y]);p4:=subs(ratk[4],[x,y]);

p1 := [sqrt(2), sqrt(2)]

p2 := [-sqrt(2), -sqrt(2)]

p3 := [-2, 2]

p4 := [2, -2]

> ell:=implicitplot(g=0,x=-3..3,y=-3..3,scaling=constrained,color=blue):

> akseli:=p->line([0,0],p);

akseli := proc (p) options operator, arrow; line([0...

> display(ell,akseli(p1),akseli(p2),akseli(p3),akseli(p4));

[Maple Plot]

> sqrt(subs(ratk[1],f));sqrt(subs(ratk[3],f));

2

2*sqrt(2)

Tässä näkyvät pikku- ja isoaksleien puolikkaiden pituudet.

Tarkistus:

> Eigenvectors(A);

_rtable[136966004], _rtable[137094300]

> 4*y1^2+2*y2^2=16;

4*y1^2+2*y2^2 = 16

> lhs(%)/16=1;

1/4*y1^2+1/8*y2^2 = 1

> iso:=sqrt(8); pieni:=2;

iso := 2*sqrt(2)

pieni := 2

3.

> restart: with(linalg):

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

> f:=(x,y,z)->(x/4)^2+(y/4)^2+(z/4)^2;

f := proc (x, y, z) options operator, arrow; 1/16*x...

> g:=x+y+z-L;

g := x+y+z-L

> Lag:=f(x,y,z)+lambda*g;

Lag := 1/16*x^2+1/16*y^2+1/16*z^2+lambda*(x+y+z-L)

> gr:=grad(Lag,[x,y,z,lambda]);

gr := vector([1/8*x+lambda, 1/8*y+lambda, 1/8*z+lam...

> solve({gr[1]=0,gr[2]=0,gr[3]=0,g=0},{x,y,z,lambda});

{z = 1/3*L, lambda = -1/24*L, y = 1/3*L, x = 1/3*L}...

Siis kaikki pätkät yhtä suuria : x=y=z=L/3.

> f(L/3,L/3,L/3); # Pinta-alojen summa:

1/48*L^2

Reunat: Tällöin jokin muuttuja = 0. Kaksi keskenään erilaista mahdollisuutta (symmetriasyistä):

1) yksi muuttuja = 0 ja 2) kaksi muutujaa = 0

ensin yksi, olkoon z:

> f(x,y,0);

1/16*x^2+1/16*y^2

> g:=x+y+0-L;

g := x+y-L

> Lag:=f(x,y,0)+lambda*g;

Lag := 1/16*x^2+1/16*y^2+lambda*(x+y-L)

> gr:=grad(Lag,[x,y,lambda]);

gr := vector([1/8*x+lambda, 1/8*y+lambda, x+y-L])

> solve({gr[1]=0,gr[2]=0,gr[3]=0},{x,y,lambda});

{x = 1/2*L, lambda = -1/16*L, y = 1/2*L}

> f(L/2,L/2,0);

1/32*L^2

Yksi pätkä (ei leikkausta):

> f(L,0,0);

1/16*L^2

Tässä on kaikki mahdolliisuudet. maksimi on tämä, jossa ei leikata ollenkaan, minimi se, jossa leikataan kolmeen samanpituiseen osaan.

4.

> restart:with(plots):with(linalg):with(LinearAlgebra):V2L:=vek->convert(vek,list):

>

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the assigned name GramSchmidt now has a global binding

Aloitetaan määrittelemällä kohdefunktio f ja gradienttifunktio g. Koska muuttujia on mukava välillä käsitellä jonona tai listana ja välillä taas vektorina, päädyin ehdottamaan kahta märitelmää: f ja fv, jossa jälkimmäisen argumenttina on vektori.

Vastaavasti tehdään g:n ja gv:n tapauksessa.

Kun nämä pikku kömpelyydet kestetään, niin loppu onkin sitten jo mukavaa.

> f:=(x,y)->2*(x^2+y^2)+x*y-5*(x+y);

f := proc (x, y) options operator, arrow; 2*x^2+2*y...

> fv:=X->2*(X[1]^2+X[2]^2)+X[1]*X[2]-5*(X[1]+X[2]);

fv := proc (X) options operator, arrow; 2*X[1]^2+2*...

> g:=[D[1](f),D[2](f)];

g := [proc (x, y) options operator, arrow; 4*x+y-5 ...

> gv:=v->g(v[1],v[2]);

gv := proc (v) options operator, arrow; g(v[1],v[2]...

Tehdään indeksoitu muttuja, johon kerätään SD:n tuottamia vektoreita. Ensin alkupiste:

> x[0]:=<1,-2>;

x[0] := _rtable[136129004]

> N:=5:for i to N do u:=-Normalize(Vector(gv(x[i-1])));
phi:=t->simplify(fv(evalm(x[i-1]+t*u)));
tmin:=solve(diff(phi(t)=0,t));
x[i]:=(evalm(x[i-1]+tmin*u));
od:

> X:=(Matrix(map(V2L,[seq(x[ii],ii=0..5)])));reitti:=convert(X,listlist);
plot(reitti,color=blue,scaling=constrained);
contourplot(f(x,y),x=-1..2,y=-1..3);
display(%,%%);

X := _rtable[136878080]

reitti := [[1, -2], [127/76, 13/19], [1, 259/304], ...

[Maple Plot]

[Maple Plot]

[Maple Plot]

>

5.

> restart:with(plots):with(linalg):with(LinearAlgebra):V2L:=vek->convert(vek,list):
f:=(x,y)->exp(x)*(4*x^2+2*y^2+4*x*y+2*y+1);
fv:=X->exp(X[1])*(4*X[1]^2+2*X[2]^2+4*X[1]*X[2]+2*X[2]+1);

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the assigned name GramSchmidt now has a global binding

f := proc (x, y) options operator, arrow; exp(x)*(4...

fv := proc (X) options operator, arrow; exp(X[1])*(...

> g:=[D[1](f),D[2](f)];gv:=v->g(v[1],v[2]);

g := [proc (x, y) options operator, arrow; exp(x)*(...

gv := proc (v) options operator, arrow; g(v[1],v[2]...

Aluksi vilkaisu:

> plot3d(f(x,y),x=-1..1,y=-1..1);

[Maple Plot]

> contourplot(f(x,y),x=-2..1,y=-2..0);

[Maple Plot]

>

Minimi näyttäisi olevan jossain ellipsimäisen käyrän sisällä. Kokeillaan alkupisteenä O:a.

> x[0]:=<0,0>;

x[0] := _rtable[135348308]

> u:=-Normalize(Vector(gv(x[0])));

u := _rtable[135488116]

> z:=evalm(x[0]+t*u);

z := vector([-1/2*t, -t])

> phi:=t->simplify(fv(evalm(x[0]+t*u)));

phi := proc (t) options operator, arrow; simplify(f...

> phi(t);

exp(-1/2*t)*(5*t^2-2*t+1)

> solve(diff(phi(t)=0,t));

11/5+4/5*sqrt(6), 11/5-4/5*sqrt(6)

> tmin:=evalf(min(%));

tmin := .240408206

> X[1]:=evalf(x[0]+tmin*u);

X[1] := _rtable[136159272]

>

>

> grad(f(x,y),[x,y]);

vector([exp(x)*(4*x^2+2*y^2+4*y*x+2*y+1)+exp(x)*(8*...

> solve({%[1]=0,%[2]=0},{x,y});

{y = 1, x = -3/2}, {y = -1, x = 1/2}

solve löytää KRP:t kunnioitettavan hyvin.

> plot3d(f(x,y),x=0..1,y=-2..0,style=patchcontour);

[Maple Plot]

> contourplot(f(x,y),x=0..1,y=-2..0,contours=20);

[Maple Plot]

Tämä ei ole ihan helppo tapaus. Viivahaun minimin määrittäminen ei sujukaan ihan mekaanisesti.

Tässä täytyisi edetä askel kerrallaan tai sitten implementoida interpolaatioon perustuva algoritmi.

No, jätetään tähän, kokeissahan ei sellaiset komplikaatiot tule kysymykseen.

6.

Kertauksen vuoksi jyrkimmän laskun menetelmän proseduuri:
Olkoon f(x) funktio, jonka minimiä haemme. valitaan alkupiste x(0).
1. minimoidaan f suoralla : x[n]-t*grad(f(x[n])) (t > 0)
2. päivitetään: x[n+1]:=x[n] + tmin*u, missä u = -grad(f(x[n])).

> restart:f:=(x,y)->x^2+c*y^2;
fv:=X->X[1]^2+c*X[2]^2;
g:=[D[1](f),D[2](f)];
gv:=v->g(v[1],v[2]);
#gv(x[0]); g(1,-2);
X0:=evalm(a*<c,1>);
evalm(-gv(X0));
u:=evalm(-gv(X0)/(2*a*c));
phi:=t->simplify(fv(evalm(X0+t*u)));
phi(t);
(c-t)^2+c*(1-t)^2: expand(%): collect(%,t^2);
collect(phi(t),t^2);
tmin:=solve(diff(phi(t),t)=0,t);
X1:=map(normal,evalm(X0+tmin*u));

Warning, the name changecoords has been redefined

f := proc (x, y) options operator, arrow; x^2+c*y^2...

fv := proc (X) options operator, arrow; X[1]^2+c*X[...

g := [proc (x, y) options operator, arrow; 2*x end ...

gv := proc (v) options operator, arrow; g(v[1],v[2]...

X0 := vector([a*c, a])

vector([-2*a*c, -2*a*c])

u := vector([-1, -1])

phi := proc (t) options operator, arrow; simplify(f...

a^2*c^2-4*a*c*t+t^2+c*a^2+c*t^2

(1+c)*t^2-4*c*t+c^2+c

(1+c)*t^2-4*a*c*t+a^2*c^2+c*a^2

tmin := 2*a*c/(1+c)

X1 := vector([a*c*(-1+c)/(1+c), -a*(-1+c)/(1+c)])

> c:=20.:seq(c*(c-1)^n/(c+1)^n,n=10..20);
seq(c*(c-1)^n/(c+1)^n,n=21..40);
seq(c*(c-1)^n/(c+1)^n,n=100..110);
c:=40.:
seq(c*(c-1)^n/(c+1)^n,n=200..201);

7.351450850, 6.651312673, 6.017854322, 5.444725340,...
7.351450850, 6.651312673, 6.017854322, 5.444725340,...

2.444839908, 2.211998013, 2.001331535, 1.810728532,...
2.444839908, 2.211998013, 2.001331535, 1.810728532,...

.9004521050e-3, .8146947616e-3, .7371047843e-3, .66...
.9004521050e-3, .8146947616e-3, .7371047843e-3, .66...

.1812216385e-2, .1723815586e-2

>

Suppenee hitaaasti, mutta nyt en jaksa selostaa enempää.