L/newtopt.mws

Newtonin menetelmä optimointiin 20.3.02 HA

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

> alias(Tr=Transpose):

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

> read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl");
#read("/p/edu/mat-1.414/maple/v202.mpl");

Newtonin menetelmä

> f:=(x1,x2)->100*(x2-x1^2)^2+(1-x1)^2; # Rosenbrockin banaanifunktio

> plot3d(f(x,y),x=-2.5..1.5,y=-8..6,view=[-2.5..1.5,-8..6,0..100],style=patchcontour);

> gr:=Transpose(Vector(grad(f(x1,x2),[x1,x2]))); # Transponoidaan tässä, niin ei tule jatkossa ongelmia pysty/vaakavektorian kanssa.

> H:=Matrix(hessian(f(x1,x2),[x1,x2]));

> p0:=<-2,5>;

> p:=p0:

> Hp:=subs(x1=p[1],x2=p[2],H);

> Eigenvalues(Hp);evalf(%);

Jotta menetelmän toiminen olisi taattua, olisi Hp:n oltava pos. def. Tässä ei aivan ole, mutta positiivinen ominaisarvo

dominoi voimakkaasti, joten "toimii todennäköisesti".

> grp:=subs(x1=p[1],x2=p[2],gr);

> h:=LinearSolve(Hp,-grp);

> p1:=p0+h;

> p1:=evalf(p1);

> p:=p1:

> Hp:=subs(x1=p[1],x2=p[2],H);

> det(Hp);

> Eigenvalues(Hp);

Tässä on jo oikeasti pos. def.

> X:='X':X[0]:=evalf(p0);

> N:=5:for i to N do
p:=X[i-1]:
Hp:=subs(x1=p[1],x2=p[2],H);
HH[i-1]:=Eigenvalues(Hp):
grp:=subs(x1=p[1],x2=p[2],gr);

> h:=LinearSolve(Hp,-grp);

> X[i]:=evalm(X[i-1]+h);
od:

>

>

> seq(evalm(X[k]),k=1..5);

>