Harj. 9 AV

27.3.02 HA

Alustukset

Kts. myös L/LSQ.mws ja L/newtopt.mws

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra): with(plots): alias(Tr=Transpose,Lsolve=LinearSolve):

> with(linalg):with(plottools):

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

> setoptions3d(axes=boxed,orientation=[-30,50]):

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

> read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl");

Teht. 1

> xd:=[2,3,4,5]: yd:=[5,9,15,2]:

> CT:=Matrix([[1,1,1,1],xd]); C:=Tr(CT);

> A:=CT.C;

>

> sum(xd[i]^2,i=1..4);

> b:=CT.Vector(yd);

> sum(yd[i],i=1..4);

> sum(xd[i]*yd[i],i=1..4);

> a:=Lsolve(A,b);

> datakuva:=plot([seq([xd[i],yd[i]],i=1..4)],style=point,symbol=circle,symbolsize=15):

> suora:=a[1]+a[2]*x;

> display(datakuva,plot(suora,x=1..6,color=blue));

>

2.

> f:=(x1,x2)->2*(x1^2+x2^2)+ x1*x2 -5*(x1+x2);

> gr:=Vector(grad(f(x1,x2),[x1,x2]));

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

> p0:=<1,-2>; p:=p0:

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

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

> p0+h;

> solve({gr[1]=0,gr[2]=0},{x1,x2});

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

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

> p0+h;

Toden totta, tehtäväpaperissa oleva ennustus näyttää toteutuvan!

No miksi?

Newtonissa otetaan funktion muutokselle Taylorin 2. asteen polynomiapproksimaatio, jonka tarkka minimi haetaan (ratkaisemalla yksikäsitteinen KRP), mimi siksi, että oletetaan pos. def.

Nyt, kun kohdefunktio on suoraan 2. asteen polynomi, se yhtyy tarkalleen omaan Taylorin 2. asteen polynomiinsa.

Newtonaskel laskee sen minimin tarkasti.

3.

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra): with(linalg):with(plots): alias(Tr=Transpose,Lsolve=LinearSolve):

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

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

> read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl");

> f1:=x^2-x*y+2*y^2-10; f2:=x^3*y^2-2;

> fv:=x-><x[1]^2-x[1]*x[2]+2*x[2]^2-10,x[1]^3*x[2]^2-2>;

> F:=(x,y)->[x^2-x*y+2*y^2-10,x^3*y^2-2];

> implicitplot({f1=0,f2=0},x=0..4,y=-3..3);

> J:=Matrix(jacobian([f1,f2],[x,y]));

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

> Jp:=subs(x=pp[1],y=pp[2],J);fp:=fv(pp);

> h:=evalf(LinearSolve(Jp,-fp));

> p[1]:=p[0]+h;

> pp:=p[1]:Jp:=subs(x=pp[1],y=pp[2],J);fp:=fv(pp);h:=evalf(LinearSolve(Jp,-fp));p[2]:=p[1]+h;

> for j to 5 do
pp:=p[j-1]:Jp:=subs(x=pp[1],y=pp[2],J);fp:=fv(pp);h:=-evalf(LinearSolve(Jp,fp));p[j]:=p[j-1]+h; od:

> p[5];

> Matrix(map(V2L,[seq(p[i],i=0..5)]));

>

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

> n:=20:for j to n do
pp:=p[j-1]:Jp:=subs(x=pp[1],y=pp[2],J);fp:=fv(pp);h:=-evalf(LinearSolve(Jp,fp));p[j]:=p[j-1]+h; od:

> #op(Matrix(map(V2L,[seq(p[i],i=0..n)]))):

> p[n];

Eipä ota supetakseen.

> p[0]:=<3,0>;

> n:=1:for j to n do
pp:=p[j-1]:Jp:=subs(x=pp[1],y=pp[2],J);fp:=fv(pp);h:=-evalf(LinearSolve(Jp,fp));p[j]:=p[j-1]+h; od:

Error, (in LinearAlgebra:-LA_Main:-LinearSolve) inconsistent system

> fp;Jp;

Jakobiaani on singulaarinen, Newton jysähtää.

Tässä olen koonnut proc:ksi 2-ulotteisen Newtonin, sitä on mukava käytellä:

> print(Newtonsys2);

> Newtonsys2(F,<1,1>,6);

> Newtonsys2(F,<3,1>,6);

> Newtonsys2(F,<3,-1>,5);

>

4.

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra): with(plots): alias(Tr=Transpose,Lsolve=LinearSolve):

> with(linalg):with(plottools):

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

> f1:=d1*cos(alpha)+d2*cos(alpha+beta)-x0;f2:=d1*sin(alpha)+d2*sin(alpha+beta)-y0;

> F:=unapply([f1,f2],x,y);

> J:=Matrix(jacobian([f1,f2],[alpha,beta]));

> JI:=MatrixInverse(J);

> p0:=<alpha,beta>;

> p1:=p0-JI.Vector(F(alpha,beta));

> det(J);

> simplify(%);

> expand(%);

> factor(%);

> map(factor@expand,JI);

> subs(sin(alpha)^2+cos(alpha)^2=1,%);

> map(combine,%);

> JI:=%;

> p1:=p0-JI.Vector(F(p0[1],p0[2]));

>

Aivan siedettävä iteraatiokaava rarkaistussa muodossakin tässä tapauksessa. (Ehkä parasta kuitenkin jättää matriisikertolaskut tehtäväksi numeerisesti.)