### /home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl 3.1.2002 # read("c:\\opetus\\maple\\v202.mpl"); (PC-WIN:ssä) # read("/p/edu/mat-1.414/maple/v202.mpl"); # read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl"); HA:n priv. ############################################ ##Datan muokkausta, yleisiä apufunktioita: ############################################ V2L:=vek->convert(vek,list): # Sama kuin v2l L2V:=lis->convert(lis,Vector): # LinearAlgebra-vektori v2l:=vek->convert(vek,list): # HAM-assa painovirhe (vec) l2v:=lis->convert(lis,vector): #Liitetään kaksi vektoria peräkkäin: vappend:=(v1,v2)->vector([op(convert(v1,list)),op(convert(v2,list))]): #Matriisin jonouttaminen pitkäksi vektoriksi m2l:=M->map(op,convert(M,listlist)): # Vastaa Matlabin M(:), mutta tässä # riveittäin #Diskretointi vrt. Matlab-tuttavamme linspace linspace:=proc() local i,n,a,b; a:=args[1];b:=args[2]; if nargs=2 then n:=10 else n:=args[3] fi; [seq(a+i*(b-a)/(n-1),i=0..n-1) ] end: # Esim: #x10:=linspace(0,1); x20:=linspace(0,1,20); #Kts. 5.3 Diskretointi [HAM] s. 117. Vaihdoimme oletusarvon 10:een, #Matlab:n 100:n sijasta. ################### #Konformikuvauksia# ################### gympyra:=(z0,r)->complexplot(z0+r*exp(I*Theta),Theta=0..2*Pi): gympyrakuva:=(z0,r,f)->complexplot(f(z0+r*exp(I*Theta)),Theta=0..2*Pi): #display([seq(gympyra(I+0.1*k,1),k=0..20)],insequence=true); #display([seq(gympyra(0,0.1*k),k=0..10)]); #gsade:=(Theta,r1,r2)->complexplot(r*exp(I*Theta),r=r1..r2): #gsadekuva:=(Theta,r1,r2,f)->complexplot(f(r*exp(I*Theta)),r=r1..r2): #display([seq(gsade(k*2*Pi/10,0.5,3),k=0..9)]); #nymp:=7: nsateet:=15: rmax:=4: hr:=rmax/nymp: hTheta:=2*Pi/nsateet: #display([seq(gsade(k*hTheta,0,rmax),k=0..nsateet),seq(gympyra(0,k*hr), #k=0..nymp)]); ############## # Iteraatiot # ############## porras:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]): newtonkulma:=x->plot([[x,0],[x,f(x)],[N(x),0]],color=black): NNN:=x->x-f(x)/D(f)(x): # Pelkkä N on liian salakavala! # N:=NNN ################################################### #Newtonin iteraatio kompleksitasossa ja fraktaalit# ################################################## newt3:=proc(z) local zz,m,j1,j2,j3; j1:=1.; j2:=evalf(exp(I*2*Pi/3)); j3:=evalf(exp(-I*2*Pi/3)); zz:=evalf(z); for m from 0 to 50 while abs(zz^3-1)>=0.001 do zz:=zz-(zz^3-1)/(3*zz^2); od; if m > 45 then RETURN(10) fi; if abs(zz-j1) < 0.1 then RETURN(0) elif abs(zz-j2) < 0.1 then RETURN(1) elif abs(zz-j3) < 0.1 then RETURN(2) else RETURN(-1) fi; end: # j1:=1.:j2:=evalf(exp(I*2*Pi/3)):j3:=evalf(exp(-I*2*Pi/3)): z2xyw:=(z,w)->[Re(z),Im(z),w]: # "Kertakäyttöfkt." # juuret3d:=spacecurve([z2xyw(j1,0),z2xyw(j2,1),z2xyw(j3,2)],style=point,symbol=circle,color=black,symbolsize=20): #with(plots): # Piirretään "alkuarvokuva", jossa kutakin juurta kohti suppeneva alkupiste # nostetaan sille tasolle, joka on juuren numero. Väritetään korkeuden mukaan. # Värikoodattu versio saadaan katsomalla suoraan ylhäältä. #aakuva:=plot3d('newt3'(x+I*y),x=-2..2,y=-2..2,grid=[40,40],view=[-2..2,-2..2,0..3],axes=FRAME,style=PATCHNOGRID,orientation=[-130,50]): # Valitse Color-valikosta Z HUE ################# # Interpolaatio # ################# # Lagrangen kertojapolynomi: L:=proc(j,xd,x) local oso,nimi,i,j1; j1:=j+1; oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j1]); nimi:=subs(x=xd[j1],oso); oso/nimi; end: #Parametrit # j -- antaa L:n indeksin, eli kyseessä on Lj # j=0..n, missä n on xd-listan pituus - 1 # xd -- xdata, [x0,x1,...,xn], lista (n+1 pistettä) # x -- polynomin muuttuja #Lagrangen polynomijono muodostetaan nyt näin, kun data on listassa xd # seq(L(j,xd,x),j=0..nops(xd)-1); ######################################################################## # L1 -- Lagrangen kertojapolynomi: versio, jossa datapisteet [x1,..,xn] # ja Lagrangen kertojapolynomit L1[1],...,L1[n] (astetta n-1). # Hyödyllinen mm. Gaussin integroinnissa. L1:=proc(j,xd,x) local oso,nimi,i; oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j]); nimi:=subs(x=xd[j],oso); oso/nimi; end: ## end L1 ####################### ################################################# # Interpolaatiopolynomi Lagrangen menetelmällä: ################################################# lagint:=(xd,yd,x)->add(yd[j+1]*L(j,xd,x),j=0..nops(xd)-1): # lagint([a,b,c],[ya,yb,yc],x); # xd:=[a,b,c]: lagint(xd,map(f,xd),x); ########################################################################## # interplot -- piirretään datapisteisiin liittyvä interpolaatiopolynomi # ja datapisteet samaan kuvaan # Tämä on niin usein esiintyvä tehtävä, että kannatti kirjoittaa proc:ksi. # HA 7.2.01 ####################################################################### interplot:=proc(xd,yd) local xyd,p,x; xyd:=zip((x,y)->[x,y],xd,yd); p:=interp(xd,yd,x); plots[display]([plot(p,x=xd[1]..xd[nops(xd)],color=blue),plot(xyd,style=point,symbol=cross,symbolsize=20,color=black)]); end: # Esim: xd:=[...]: f:=x->...: yd:=map(f,xd): # interplot(xd,yd); # display([interplot(xd,yd),f(x),x=xd[1]..xd[nops(xd)]]); ###################################################################### # Pelkkä datan piirto, tätä tarvitaan usein. (Tässä on riisuttu interplotista # muut osat pois ) # Data annetaan kahtena listana; xd ja yd # dataplot:=proc(xd,yd) local xyd; xyd:=zip((x,y)->[x,y],xd,yd); plot(xyd,style=point,symbol=cross,symbolsize=20,color=black); end: ## Numeerinen integrointi### # ## 1) Gaussin menetelmä ################# gaussolmut:=k->fsolve(orthopoly[P](k,_xxx_)=0,_xxx_): gausspainot:=n->seq(int(L1(_i_,[gaussolmut(n)],_x_),_x_=-1..1),_i_=1..n): # Tällaisissa yhden rivin (->) määrittelyissä ei lokaaleja muuttujia # saa mukaan. Tällöin on syytä joko tehdä kömpelömpi proc-rakenne # tai tehdä lokaaleiksi tarkoitetuista "outoja", kuten _xxx_ # (pelkkä _x_ on jo aika kelvollisen outo) . gausstaulukko:=n->[[gaussolmut(n)],[gausspainot(n)]]: gaussint:=proc(n,f) # Suorita ensin gtaulu:=gausstaulukko(n); global gtaulu; local s,w,i; s:=gtaulu[1]; w:=gtaulu[2]; sum(w[i]*f(s[i]),i=1..n); end: # n:=3: gtaulu:=gausstaulukko(n): # seq([gaussint(3,x->x^k),evalf(int(x^k,x=-1..1))],k=1..2*n-1); # Olis järkevämpää: gaussint1:=proc(f) global gtaulu; local s,w,i; s:=gtaulu[1]; w:=gtaulu[2]; sum(w[i]*f(s[i]),i=1..nops(s)); end: # Kaikkein järkevintä on panna gtaulu parametrilistaan: # Olkoon tämä se lopullinen 16.4.02 Gaussint:=proc(taulukko,f) local s,w,i; s:=taulukko[1]; w:=taulukko[2]; sum(w[i]*f(s[i]),i=1..nops(s)); end: # ----------- end Gaussint ----------------- # ---------- Gaussint-esim. ----- # n:=3: gtaulu:=gausstaulukko(n): Gaussint(gtaulu,f); # # Voidaan haluttaessa kutsua myös suoraan: # Gaussint(gaustaulukko(n),f); # # ---------end Gaussint-esim. ----- # # Integrointi yli mieliv. välin [a,b] # gaussab:=(n,f,a,b)->gaussint(n,t->f(1/2*t*b-1/2*t*a+1/2*b+1/2*a)*(1/2*b-1/2*a)): # Usean muuttujan integrointi # Aluksi yli neliön: gaussintnelio:=proc(m,n,f) # Suorita ensin gtaulux:=gausstaulukko(m); # gtauluy:=gausstaulukko(n); global gtaulux,gtauluy; local sx,wx,sy,wy,i,j; sx:=gtaulux[1]; wx:=gtaulux[2]; sy:=gtauluy[1]; wy:=gtauluy[2]; sum(wx[i]*sum(wy[j]*f(sx[i],sy[j]),j=1..n),i=1..m); end: # gtaulux:=gausstaulukko(3); # gtauluy:=gtaulux: #gaussintnelio(3,3,(x,y)->x^4*y^4),int(int(x^4*y^4,x=-1..1),y=-1..1); #evalf(%); # Tästä on hyvä jatkaa. Gaussnelio:=proc(gtaulux,gtauluy,f) # Suorita ensin gtaulux:=gausstaulukko(m); # gtauluy:=gausstaulukko(n); local sx,wx,sy,wy,i,j,m,n; sx:=gtaulux[1]; wx:=gtaulux[2]; sy:=gtauluy[1]; wy:=gtauluy[2]; m:=nops(sx); n:=nops(sy); sum(wx[i]*sum(wy[j]*f(sx[i],sy[j]),j=1..n),i=1..m); end: # gtaulux:=gausstaulukko(3); # gtauluy:=gtaulux: #gaussintnelio(3,3,(x,y)->x^4*y^4),int(int(x^4*y^4,x=-1..1),y=-1..1); #evalf(%); # Tästä on hyvä jatkaa. ################ # Yksiulotteinen: Gaussab:=(taulu,f,a,b)->Gaussint(taulu,t->f(1/2*t*b-1/2*t*a+1/2*b+1/2*a)*(1/2*b-1/2*a)): # Kaksiulotteinen: Gaussxy:=(gtaulux,gtauluy,f,a1,b1,a2,b2)->Gaussab(gtaulux,x->Gaussab(gtauluy,y->f(x,y),a2(x),b2(x)),a1,b1): # Suorita ensin gtaulux:=gausstaulukko(m);gtauluy:=gausstaulukko(n); # x-projisoituva alue # a1 ja b1 ovat x-vakiorajat # a2 ja b2 ovat yhden muutt. funktioita # Alue: a1 <= x <= b1, a2(x) <= y <= b2(x) # Kolmiulotteinen: Gaussxyz:=(gtaulux,gtauluy,gtauluz,f,a1,b1,a2,b2,a3,b3)->Gaussxy(gtaulux,gtauluy,(x,y)->Gaussab(gtauluz,z->f(x,y,z),a3(x,y),b3(x,y)),a1,b1,a2,b2): # Suorita ensin gtaulux:=gausstaulukko(m);gtauluy:=gausstaulukko(n); # gtauluz:=gausstaulukko(p); # xy-projisoituva alue # a1 ja b1 ovat x-vakiorajat # a2 ja b2 ovat yhden muutt. funktioita # a3 ja b3 ovat 2 muutt. fkt. # Alue: a1 <= x <= b1, a2(x) <= y <= b2(x), a3(x,y) <= z <= b3(x,y) # Monte Carlo-menetelmä # # Tehdään ensin 2d-versio, jossa integroimisalue ympäröidään # suorakulmiolla T. # Ehto alueeseen kuulumiselle esitetään tyyliin # g(x,y) <= 1 # Tätä pitää tarpeen mukaan modifioida. # # MonteCarlo2d:=proc (f, g, T, N) local k, summa, lkm, xv, yv, x, y, lx, ly, A; xv := T[1]; yv := T[2]; lx := xv[2]-xv[1]; ly := yv[2]-yv[1]; A := lx*ly; summa := 0.; lkm := 0; for k to N do x := xv[1]+lx*satu(); y := yv[1]+ly*satu(); if evalf(g(x,y)) <= 1. then summa := summa+f(x,y); lkm := lkm+1 end if end do; A*summa/N end proc: # Argumentit: f -- integroitava funktio, esim: f:=(x,y)->x*y; # g -- ehto, esim: g:=(x,y)->x^2+y^2; # T -- ympäröivä suuntaissärmiö, esim: # T:=[[-1,1],[-1,1]]; # N -- Kuinka monta pistettä arvotaan. # Alustus: # satu:=rand(0..10^6)/10.^6; # # Tehdään 3d-versio, jossa integroimisalue V ympäröidään # # suuntaissärmiöllä T. # Ehto alueeseen kuulumiselle esitetään tyyliin # g(x,y,z) <= 1 # Näitä pitää tarpeen mukaan modifioida. # Vertailun vuoksi Riemannin summa: # Samat argumentit muuten, mutta Riemannin N vastaa Monte Carlon N^2 # Tässä ei mitään satunnaislukualustusta. Riemann2d:=proc(f,g,T,N) local k,summa,lkm,xv,yv,x,y,lx,ly,A,X,Y,i,j; xv:=T[1]; yv:=T[2]; lx:=xv[2]-xv[1]; ly:=yv[2]-yv[1]; A:=lx*ly; X:=linspace(op(xv),N); Y:=linspace(op(yv),N); # xy:=seq(seq([X[i],Y[j]],j=1..N),i=1..N); summa:=0.; lkm:=0; for i to N-1 do for j to N-1 do if evalf(g(X[i],Y[j])) <= 1. then summa:=summa+f(x,y); lkm:=lkm+1; end if; end do; end do; A*summa/N^2; end: # # MonteCarlo3d:=proc(f,g,T,N) local k,summa,lkm,xv,yv,zv,x,y,z,lx,ly,lz,V; xv:=T[1]; yv:=T[2]; zv:=T[3]; lx:=xv[2]-xv[1]; ly:=yv[2]-yv[1];lz:=zv[2]-zv[1]; V:=lx*ly*lz; summa:=0.; lkm:=0; for k to N do x:=xv[1]+lx*satu(); y:=yv[1]+ly*satu(); z:=zv[1]+lz*satu(); if evalf(g(x,y,z)) <= 1. then summa:=summa+f(x,y,z); lkm:=lkm+1; end if; end do; V*summa/N; end: # Argumentit: f -- integroitava funktio, esim: f:=(x,y,z)->x*y*z; # g -- ehto, esim: g:=(x,y,z)->z^2+(sqrt(x^2+y^2)-3)^2; # T -- ympäröivä suuntaissärmiö, esim: # T:=[[1,4],[-3,4],[-1,1]]; # N -- Kuinka monta pistettä arvotaan. # Alustus: # satu:=rand(0..10^6)/10.^6; ################################################################ ### Sarjat ### # # Termien järjestäminen # Järjestetään vuorottelevan harmonisen sarjan termit järjestykseen: # (x(1)+x(3)+x(2)) + (x(5)+x(7)+x(4)) + ... # Kaksi paritonta + yksi parillinen (kaksi pos + 1 neg.) # Palautetaan indeksivektori # jarj132:=proc(n) local ind,iplus,iminus,k; iplus:= seq(2*k-1,k=1..3*n); iminus:=seq(2*k,k=1..n); ind:=map(op,[seq([iplus[2*k-1],iplus[2*k],iminus[k]],k=1..n)]); end: # x:=k->(-1)^(k+1)/k; # ind:=jarj132(8);nops(ind); # add(x(i),i=ind):evalf(%); # add(x(i),i=1..nops(ind)):evalf(%); # jarj132(10);sort(%); ## ######################################################## # 3d-grafiikkaa # Lähde: University of Illinois # zoom:=(f,x0,y0,w)->plot3d(f(x,y),x=x0-w/2..x0+w/2,y=y0-w/2..y0+w/2,axes=BOX): # Esim: f:=(x,y)->x*y; zoom(f,0,0,1);zoom(f,0,0,1/2);zoom(f,0,0,1/4); # ## Tangenttitaso: Israel s. 156, annetaan myös f argumentiksi. # tangtaso:=(f,a,b,hmax)->plot3d(f(a,b)+D[1](f)(a,b)*(x-a)+D[2](f)(a,b)*(y-b),x=a-hmax..a+hmax,y=b-hmax..b+hmax,style=PATCH,color=red,grid=[7,7]): ###################################################### # # Helppokäyttöiset funktiot nuolen piirtämiseen O:sta kohdepisteeseen. # snuoli:=u->arrow([0,0],u,0.01,0.05,0.02,color=blue): pnuoli:=u->arrow([0,0],u,0.01,0.05,0.02,color=red): # # Esim: with(plottools):with(plots): # n1:=snuoli([0,1]):n2:=snuoli([1,0]):n3:=pnuoli([1,1]): # display([n1,n2,n3]); # # Yleisemmin: nuoli:=(alku,loppu,vari)->arrow(alku,loppu,0.01,0.05,0.02,color=vari): # # Esim: display([nuoli([0,0],[0,1],red),nuoli([0,1],[1,1],black), # nuoli([1,1],[0,0],blue)]); # # Gradientti, suunnattu derivaatta, ym. Pieni paketti, joka ei käytä linalg- # pakkausen valmista gradient-funktiota. # Otetaan käyttöön Israel:n evl ja ennen kaikkea &. # evl:=V->convert(evalm(V),list): # Vektori -> lista #`&.`:=(U,V)->linalg[dotprod](U,V,orthogonal): `&.`:=(U,V)->linalg[dotprod](U,V): `&x`:=(U,V)->convert(linalg[crossprod](U,V),list): Gr2:=(f,a,b)->: # 2. muutt. fkt. gradientti LinearAlgebra-tyyli Gr2L:=(f,a,b)->convert(,list): gr2:=f->vector([D[1](f),D[2](f)]): # linalg-tyyli #n:=v->sqrt(v &. v): # Normi (eukl.) Sd2:=(f,a,b,u)->Gr2(f,a,b).LinearAlgebra[Normalize](u,2): # u sarakevektori # LinearAlgebra-tyyli sd2:=(f,a,b,u)->gr2(f)(a,b)&.linalg[normalize](u): # Suunnattu derivaatta # linalg-tyyli jana2dpp:=(p1,p2)->plot([p1,p2],thickness=2,color=green): jana2dpv:=(p,v)->plot([p,evl(vector(p)+vector(v))],thickness=2,color=blue): jana3dpp:=(p1,p2)->spacecurve([p1,p2],thickness=2,color=green): jana3dpv:=(p,v)->spacecurve([p,evl(vector(p)+vector(v))],thickness=2,color=blue): ################ Optimointia ja epälin sys. ################## # # Apufunktio: evm:=proc(mat,symblist,valulist) local Mat,i; Mat:=mat; if type(mat,list) then Mat:=convert(mat,vector) fi; for i to nops(symblist) do Mat:=subs(symblist[i]=valulist[i],op(Mat)) od; map(eval,Mat); end: # Esim. # J:=linalg[jacobian](flist,xlist);Jx0:=evm(J,xlist,x0); ## Elegantimpi tapa, wh.werner@t-online.de (Wilhelm Werner): # with(linalg): # f:=(x,y)->[x^2+sin(y),x*y+tan(x+y)]: # J:=map(unapply,jacobian(f(x,y),[x,y]),x,y): # J(3,1); # # # Optimointi: # Jyrkimmän laskeutuman menetelmä - Steepest descent # SDstep:=proc(f,p0,Lmax) local vars,_x,g,s,fs,lambda,i,t,p,dp; vars:=[seq(_x[i],i=1..nops(p0))]; g:=linalg[grad](f(vars),vars); s:=-normalize(evm(g,vars,p0)); fs:=unapply(f(evalm(p0+lambda*s0)),lambda); lambda:=Lmax; while(fs(lambda)>fs(0)) do lambda:=lambda/2. od; p:=interp([0,lambda/2,lambda],map(fs,[0,lambda/2,lambda]),t); dp:=diff(p,t); lambda:=solve(dp=0,t); [v2l(evalm(p0+lambda*s0)),fs(lambda)]; end: ## SD2 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:=

; 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:=; 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]; # Newton2(F,<-1,0.5>,4); local F,xc,X,i,Jc,Fc,h,J,x,y; F:=Fun(x,y); xc:=evalf(x0); X:=<>; 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:=; od: convert(Transpose(X),listlist); end: # # plotonlypath -- Reitinpiirto. Muuten sama kuin plotpath, mutta ei piirrä # korkeuskäyriä. # Soveltuu erit. optimointiin. # Reitti on 2-rivinen matriisi, esim. Newton2:n, SDstep2:n ym. palauttama # tulos. # 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); # plotpath(F,X,xa[1]-1..xa[1]+1,xa[2]-1..xa[2]+1); # display(%%,view=[1.9..2.1,0.2..0.7]); # # plotpath -- Reitinpiirto. # Reitti on 2-rivinen matriisi, esim. Newton2:n palauttama tulos. # plotpath:=proc(F,X,xrange,yrange) local reitti,reitinkuva,reittinumerot,rnumkuva,kuvaajat,f1,f2,x,y; f1:=F(x,y)[1];f2:=F(x,y)[2]; kuvaajat:=implicitplot({f1=0,f2=0},x=op(1,xrange)..op(2,xrange),y=op(1,yrange)..op(2,yrange)): reitti:=convert(Transpose(X),listlist); 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([kuvaajat,reitinkuva,rnumkuva]); display([kuvaajat,reitinkuva,rnumkuva],view=[xrange,yrange]); end: # Esim: # xa:=<1,1>: # X:=Newton2(F,xa,5); # plotpath(F,X,xa[1]-1..xa[1]+1,xa[2]-1..xa[2]+1);alkupiste=xa; # display(%%,view=[1.9..2.1,0.2..0.7]); # Newton3:=proc(Fun,x0,niter) # Fun on 3. muuttujan funktio # x0 alkupiste, niter iteraatioiden lkm. # Esim: # F:=(x,y,z)->[y*exp(x)-2,x^2+y-4,..]; # Newton3(F,<-1,0.5>,4); local F,xc,X,i,Jc,Fc,h,J,x,y,z; F:=Fun(x,y,z); xc:=evalf(x0); X:=<>; J:=Matrix(linalg[jacobian](F,[x,y,z])); for i from 0 to niter do Jc:=subs(x=xc[1],y=xc[2],z=xc[3],J); Fc:=Vector(subs(x=xc[1],y=xc[2],z=xc[3],F)); h:=LinearAlgebra[LinearSolve](Jc,-Fc); xc:=xc+h; X:=; od: end: ## Pienimmän neliösumman sovitus ### PNS:=(C,ydata)->convert(LinearSolve(Transpose(C).C,Transpose(C).Vector(yd)),list): # Esim: # xd:=evalf([1940, 1950, 1960, 1970, 1980, 1990]); # yd:=[132.165, 151.326, 179.323, 203.302, 226.542, 249.633]; # C2:=Van(xd,nops(xd),3); kertoimet:=PNS(C2,yd); # poly2:=sum(kertoimet[i]*x^(i-1),i=1..3); # ######################### # Grafiikkafunktioita # ######################### ## Koordinaattiristikko 19.4.01 vaakajana:=(x1,x2,y)->plot([[x1,y],[x2,y]],thickness=1,color=grey): pystyjana:=(y1,y2,x)->plot([[x,y1],[x,y2]],thickness=1,color=grey): hila:=(a,b,m,c,d,n)->display([seq(vaakajana(a,b,c+_k*(d-c)/n),_k=0..n),seq(pystyjana(c,d,a+_kk*(b-a)/m),_kk=0..m)]): #Käytämme seq-muuttujia _k, _kk,koska -> tyylissä ei voi lokalisoida # (tietääkseni) # Esim: # hila(-1,1,20,3,10,4); # (siistiä!) ########### Vektorikentät 19.4.01 ################ # # Käyräintegraali cint3:=proc(F,phi,a,b) # Lasketaan R^3:n käyräintegraali, tulos: [symbolinen, numeerinen] # Modifioitu versio fcint3 alempana. # F -- Vektorikenttä: funktiolista (esim. alla) # phi -- Parametrisointifkt. phi:[a,b] -> R^3 # a,b -- Integrointirajat # Esim: F:=map(unapply,[-y,-x*y,0],x,y,z); phi:=t->[cos(t),sin(t),0]; # cint3(F,phi,0,Pi/2); # Tasokenttä: 3. komponentti = 0 # local intfun,k,t,dphi,tarkka,numappr; dphi:=map(diff,phi(t),t); intfun:=add(F[k](op(phi(t)))*dphi[k],k=1..3); tarkka:=int(intfun,t=a..b); numappr:=evalf(Int(intfun,t=a..b)); [tarkka,numappr]; end: fcint3:=proc(F,phi,a,b) # Lasketaan R^3:n käyräintegraali pelkästään numeerisesti. # Modifioitu versiosta cint3 . # F -- Vektorikenttä: funktiolista (esim. alla) # phi -- Parametrisointifkt. phi:[a,b] -> R^3 # a,b -- Integrointirajat # Esim: F:=map(unapply,[-y,-x*y,0],x,y,z); phi:=t->[cos(t),sin(t),0]; # fcint3(F,phi,0,Pi/2); # Tasokenttä: 3. komponentti = 0 # local intfun,k,t,dphi; dphi:=map(diff,phi(t),t); intfun:=add(F[k](op(phi(t)))*dphi[k],k=1..3); evalf(Int(intfun,t=a..b)); end: scint3:=proc(F,phi,a,b) # Lasketaan R^3:n käyräintegraali pelkästään symbolisesti. # Modifioitu versiosta cint3 . # F -- Vektorikenttä: funktiolista (esim. alla) # phi -- Parametrisointifkt. phi:[a,b] -> R^3 # a,b -- Integrointirajat # Esim: F:=map(unapply,[-y,-x*y,0],x,y,z); phi:=t->[cos(t),sin(t),0]; # fcint3(F,phi,0,Pi/2); # Tasokenttä: 3. komponentti = 0 # local intfun,k,t,dphi; dphi:=map(diff,phi(t),t); intfun:=add(F[k](op(phi(t)))*dphi[k],k=1..3); int(intfun,t=a..b); end: # # Potentiaalifunktio 23.4.01 # potentiaali:=proc(F) local x,y,z; # Tarkistettava pyörteettömyys ja 1-yhtenäisyys. int(F[1](t,0,0),t=0..x)+int(F[2](x,t,0),t=0..y)+ int(F[3](x,y,t),t=0..z); unapply(%,x,y,z); end: # Esim: # F:=map(unapply,[2*x*y*z,x^2*z,x^2*y],x,y,z); # potentiaali(F); # grad(%(x,y,z),[x,y,z]); potvers2:=proc(F,p0) local x,y,z; # Tarkistettava pyörteettömyys ja 1-yht. int(F[1](t,p0[2],p0[3]),t=p0[1]..x)+int(F[2](x,t,p0[3]),t=p0[2]..y)+int(F[3](x,y,t),t=p0[3]..z ); unapply(%,x,y,z); end: