H/harj9ohje.mws

3.4.02 HA

Ohjeita PNS-sovitukseen on LSQ.mws :sä ja optimointiin SD.mws :ssä sekä newtopt.mws :ssä. Näitä on päivitetty,

mws-versiot ovat paremmin ajantasalla kuin html:t.

Apufunktioita

Optimointiyöarkeilla esiintyviä komentojonoja on kirjoitettu proseduureiksi. Ne saa käyttöön lukemalla:

> #read("/p/edu/mat-1.414/maple/v202.mpl"):

> #read("c:\\maple\\v202.mpl"): # Windows:ssa

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

Tekstitiedoston v202.mpl voi myös ottaa www-sivulta ../ maple/,

jolloin sen voi lukea siitä hakemistosta, jonne sen

tallettaa. (Joku häiriö editoinnissa aina välillä.)

Otetaan funktiot kuitenkin tähän mukaan, jotta kaikki varmasti ne saavat.

> SD2:=proc(f,p0,niter)
local fv,xx,g,gv,Lmax,X,i,u,phi,lambda,p,dp,tmin;
Lmax:=1;
fv:=unapply(f(xx[1],xx[2]),xx);
g:=[D[1](f),D[2](f)];
gv:=v->g(v[1],v[2]);
X[0]:=Vector(p0);
for i to niter do
u:=-Normalize(Vector(gv(X[i-1])));
phi:=t->simplify(fv(evalm(X[i-1]+t*u)));
lambda:=Lmax;
while(phi(evalf(lambda))>phi(0.0)) do
lambda:=lambda/2. od;
p:=interp([0,lambda/2,lambda],map(phi,[0,lambda/2,lambda]),t);
dp:=diff(p,t);
tmin:=solve(dp=0,t);
X[i]:=(evalm(X[i-1]+tmin*u));
od;
[seq(convert(X[i],list),i=0..niter)];
end:

> Newtopt2:=proc(f,p0,niter)
local fv,p,X,x1,x2,h,H,gr,i,grp,Hp;
# fv:=unapply(f(X[1],X[2]),X);
gr:=Vector(linalg[grad](f(x1,x2),[x1,x2]));
H:=Matrix(linalg[hessian](f(x1,x2),[x1,x2]));
p:=Vector(evalf(p0)); X:=<p>;
for i to niter
do
grp:=subs(x1=p[1],x2=p[2],gr);
Hp:=subs(x1=p[1],x2=p[2],H);
h:=LinearSolve(Hp,-grp);
p:=p+h;
X:=<X|p>;
od;
convert(Transpose(X),listlist);
end:



Newtonsys2:=proc(Fun,x0,niter)
# Fun on 2. muuttujan funktio
# x0 alkupiste, niter iteraatioiden lkm.
# Esim:
# F:=(x,y)->[y*exp(x)-2,x^2+y-4];
# Newtonsys2(F,<-1,0.5>,4);
local F,xc,X,i,Jc,Fc,h,J,x,y;
F:=Fun(x,y);
xc:=evalf(x0); X:=<<x0>>;
J:=Matrix(linalg[jacobian](F,[x,y]));
for i from 0 to niter do
Jc:=subs(x=xc[1],y=xc[2],J);
Fc:=Vector(subs(x=xc[1],y=xc[2],F));
h:=LinearAlgebra[LinearSolve](Jc,-Fc);
xc:=xc+h;
X:=<X | xc>;
od:
convert(Transpose(X),listlist);
end:

> piirrareitti:=proc(X,xrange,yrange)
local reitti,reitinkuva,reittinumerot,rnumkuva;
#reitti:=convert(Transpose(X),listlist);
reitti:=X;
reitinkuva:=plot(reitti,style=line,symbol=circle,color=brown):
reittinumerot:=seq([op(reitti[i]),` `||`i`],i=1..nops(reitti));
rnumkuva:=textplot([reittinumerot],align=RIGHT):
#display([reitinkuva,rnumkuva]);
#display([reitinkuva,rnumkuva],view=[xrange,yrange]);
display([reitinkuva],view=[xrange,yrange]);
end:

# Joskus on mukava numeroida reittipisteet, joskus taas ei.
# Molempiin oma funktionsa, numeroversio päättyköön n:ään.

piirrareittin:=proc(X,xrange,yrange)
local reitti,reitinkuva,reittinumerot,rnumkuva;
#reitti:=convert(Transpose(X),listlist);
reitti:=X;
reitinkuva:=plot(reitti,style=line,symbol=circle,color=brown):
reittinumerot:=seq([op(reitti[i]),` `||`i`],i=1..nops(reitti));
rnumkuva:=textplot([reittinumerot],align=RIGHT):
#display([reitinkuva,rnumkuva]);
display([reitinkuva,rnumkuva],view=[xrange,yrange]);
#display([reitinkuva],view=[xrange,yrange]);
end:

# Esim:
# xa:=<1,1>:
# X:=Newtonsys2(F,xa,5);
# piirrareitti(F,X,xa[1]-1..xa[1]+1,xa[2]-1..xa[2]+1);
# display(%%,view=[1.9..2.1,0.2..0.7]);

>

1 ja 2

Kts. LSQ.mws/html

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra):alias(Tr=Transpose,Van=VandermondeMatrix):

> with(plots):

Sopivaa on varmaankin tehdä teht. 1 vaiheittain siihen tapaan kuin LSQ-työarkilla,

Tehtävässä 2 voisi olla kätevää määritellä funktio:

> PNS:=(C,ydata)->convert(LinearSolve(Transpose(C).C,Transpose(C).Vector(yd)),list);

> # C:=Van(xd,nops(xd),deg+1); kertoimet:=PNS(C,yd); #poly:=sum(kertoimet[i]*x^(i-1),i=1..deg+1);

Tässä on syytä suorittaa:

> xd:=evalf([1940, 1950, 1960, 1970, 1980, 1990]);

Homma hoituu siististi enempää skaalaamattakin, mutta järkevää tällaiaissa tapauksissa on skaalata, esim. seuraava on luoneteva skaalausfunktio:

> T:=x->(x-1965)/25;

> td:=map(T,xd);td:=evalf(%);n:=nops(td):

Laskun jälkeen muunnetaan polynomi takaisin käänteisfunktiolla alkuperäiseen x-muuttujaan, jos halutaan. Toisaalta polynomimallia käytetään arvojen laskemiseen. Riittää siis antaa x ja soveltaa siihen T-funktiota sekä laskea tällä arvolla

sovituspolynomi.

Olkoon tämä osa kuitenkin "vapaaehtoista" lisäpiirrettä, joka kannattaa muistaa käytännön laskennan tullessa eteen.

3.

Ehkä sovitaan, että otetaan ensinmainittu funktio ensisijaiseksi. Minimin suhteen kannattaa ensin katsoa pelkkää polynomia,

siitä pääsee vaikkapa suoraan neliöiksi täydentämällä johtopäätökseen jopa globaalista minimistä, jota exp(x):llä kertominen ei hetkauta. Maksimin suhteen taitaa olla niin ja näin, globaalia nyt ei mitenkään ainakaan.

Vaikka tehtävää mainostetaan Newtonina, kannattaa ottaa SD/Newton-yhdistelmä käyttään. SD.mws:ssä on ihan valmista esimerkkiainesta. Newtonissahan on se ongelma, että alkupisteen löytäminen ei ola ihan helppoa. SD etenee hitaasti, mutta varmasti pitkin reittiä, jossa funktion arvot koko ajan pienenevät. Jossain kohdassa kannattaa vaihtaa Newtoniin.

4.

Tässä sopii yllä oleva Newtonsys2. Tarkoitus on, että osaat myös selittää, mitä ohjelmakoodissa tapahtuu.

Alkupisteitä täytyy hakea piirtämällä, ratkaisupisteitä voi olla paljonkin. Newtonia ei tarvitse soveltaa kaikkiin kuvasta löytyviin, kunhan nyt johonkin.

5.

Ehdotin vaihdettavaksi: d1=6 ja d2=5, mutta voi kokeilla kumminpäin vaan.

Oma AV-teht. 4:n ratkaisuni oli tämän tyylinen. (Kts. seuraavaa kuvaa.) Silloin: B:n koordinaatit saadaan alla olevan mukaisesti.

> restart:

Warning, the name changecoords has been redefined

> with(LinearAlgebra):with(plots): with(plottools):

Warning, the name arrow has been redefined

> OO:=[0,0]: A:=[d1*cos(alpha),d1*sin(alpha)]; Av:=Vector(A);

> Bv:=Av+<d2*cos(alpha+beta),d2*sin(alpha+beta)>; B:=convert(Bv,list);

Piirretään kuva konkreettisilla numeroarvoilla.

> with(plottools):V2L:=V->convert(V,list);

> d1:=6;d2:=5;

> alpha:=Pi/6; beta:=2*Pi/3;

> display(plot([OO,A,B],x=0..6),plot(A[2],5..6,linestyle=5,color=grey),textplot([[.5,.2,'al'],[5.2,3.5,bet],[5.4,3.1,al],[5.1,2.8,'A'],[1,5.7,'B',align=TOP]],align=RIGHT),line(OO,V2L(2*Av)),plot([B],style=point,symbol=circle,symbolsize=20));

Oma AV-teht. 4:n ratkaisuni oli tämän tyylinen. Silloin: B:n koordinaatit ovat siis:

Tässä on samalla piirtämismallia. (Komentolitania on pitkä, koska rupesin tekstejä taiteilemaan.)

>

>