#Tässä on listaus ohjelmakoodeista, jotka voi hakea alla olevasta #URL:stä. Koneella olevaan tiedostoon tulee lisäyksiä ja parannuksia #nopemmin kuin kirjan sivulle. ## Tämä elää ja on jo kehittynyt kirjan versioista. ########################################################## ## Viimeinen päivitys 17.3. 1999 ## ########################################################## #Jos haluat saada nämä funktiot käyttöösi, niin hae ne www-selaimesi #{\tt File}-valikon {\tt Save as} valinnan kautta. #Kaikki proseduurit saat käyttöön \maple-istunnossa komentamalla # #> read(`ohjelmat.mpl`); #Alla tiedoston kommenteissa kerrotaan, miten polku annetaan \win-ympäristössä, #\unix:ssa se tehdään, kuten aina. #Tämänhetkinen ohjelmalistaus on puutteellinen ja keskeneräinen. #Tiedostoa ylläpidetään ja sen voi ladata www-sivun kautta. # # http://www.math.hut.fi/~apiola/maple/opas/ohjelmat.mpl #========================================================= # # Download this file into your working directory. # In your Maple session do # > read("path/ohjelmat.mpl"); # In pc the path syntax is of the form # > read("c:\\windows\\mymaple\\ohjelmat.mpl"); # # Yleiskäyttöisiä apufunktioita # # Vector-list-conversions etc.: ############################### v2l:=vek->convert(vek,list): # kirjassa 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))]): m2l:=M->map(op,convert(M,listlist)): # rowwise "ravel" the # matrix into a long list # Compare to Matlab's M(:) which does it columnwise # #### Diskretointi vrt. Matlab:n 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 s. 117. Vaihdoimme oletusarvon 10:een, Matlab:n 100;n # sijasta. # #### ITERAATIOT ################ # porras := x -> plot([[x, x], [x, f(x)], [f(x), f(x)]]): # Kts. 2.5 Funktioiden yhdistäminen ja iterointi, ss. 64 - 67 (Kuva 2.1 s. 67) # # Newtonin menetelmä # ================== # 1-dim: # ====== N:=x->evalf(x-f(x)/D(f)(x)): # Käyttöesimerkki: # f:=x->x*exp(-x)-abs(sin(10*x)); # # x1:=(N@@5)(0.25); # # Operaattorityyli: MakeNewtonstep:=proc(flauseke::algebraic,x::name) local askel; askel:=simplify(x-flauseke/diff(flauseke,x)); unapply(askel,x); end: # Esim. # f:=x-2*sqrt(x); # N:=MakeNewtonstep(f,x); # Iteroidaan, kuten edellä. # #Auxiliary function, useful in implementing multidim Newton: # Multidim Newton can be found on www-page 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); ########### # Ngen -- Generate Matlab code for 1-dim Newton # Things to add: # - include file name in parameter list # - Make "terminal" default (or N.m) # Look at example worksheets in # /p/edu/mat-1.192/k98/luennot-ja-harjoitukset/ # especially harj2sol3.mws # ########### Ngen:=proc(f,x) local df; df:=diff(f,x); interface(echo=0); system(`rm N.m`); writeto(`N.m`); # Comment for testing lprint(`function x1=N(x)`); lprint(`f=`,f,`;`); lprint(`df=`,df,`;`); lprint(`x1 = x -f ./df;`); writeto(terminal); end: ############# 5.4. Numeerinen integrointi ################# # 5.4 s. 117 -> # Huom: Kirjan kaavassa on virhe jäännöstermissä, tässä korjattu (2. der.) trapetsi:=proc(f,a,b,n) local h,S,R; h:=(b-a)/n; S:=evalf((h/2)*(f(a)+2*sum(f(a+j*h),j=1..n-1)+f(b))); R:=-(h^2*(b-a)/12)*(D@D)(f)(xi); [S,R] # Palutetaan lista: [Intappr, virhe] end: # Esim: # f:=x->exp(x)*cos(x); # trapetsi(f,0,Pi,10); # # Ekonomisempaa saattaa olla tarkastella ensin virhearviota ja vasta # virhatarkastelun jälkeen laskea integraalin approksimaatio: trapetsivirhe:=proc(f,a,b,n) local h; h:=(b-a)/n; -(h^2*(b-a)/12)*(D@D)(f)(xi) end: # # Vastaavalla tavalla voidaan toteuttaa Simpsonin kaava. # ### 5.4. Interpolaatio ######### # #lag:=(xd,x)->seq(product((x-xd[k])/(xd[j]-xd[k]),k=1..j-1)* # product((x-xd[k])/(xd[j]-xd[k]),k=j+1..nops(xd)),j=1..nops(xd)); # "Tuotantoversio kirjoitetaan proc-syntaksilla, jotta voidaan lokalisoida # indeksimuuttujat. lag:=proc(xd,x) local k,j; seq(product((x-xd[k])/(xd[j]-xd[k]),k=1..j-1)* product((x-xd[k])/(xd[j]-xd[k]),k=j+1..nops(xd)),j=1..nops(xd)) end: ######## # # Esim 1: Pisteisiin [a,b,c] liittyvät Lagrangen kertojapolynomit. # Huom: Ne riippuvat vain "x-datasta". # lagjono:=lag([a,b,c],x); # Niiden avulla voidaan kirjoittaa interpolaatiopolynomi, joka datapisteissä # saa halutut arvot, tässä y[0],y[1], y[2] . # # p:=sum(y[i-1]*lagjono[i],i=1..3); # # Esim 2: # lag(linspace(-1,1),x); # Väli [-1,1] jaetaan 10 osaan ja muod. # # jakopisteisiin liittyvät lag-polynomit. # # HUOM! Maplessa on myös valmis interp, jota kutsutaan: # interp(xdata,ydata,x); # # Yllä oleva lag on ennen kaikkea opettelutarkoitukseen (ja lagrangen # polynomeja voidaan tarvita muissakin yhteyksissä. # Toinen seikka on, että interp-funktioon on rakennettu sisään # expand, # mikä ei läheskään aina ole järkevää. # lagjono jättää polynomin "Lagrangen muotoon" # # Vaihtoehtoinen tapa ########## # # Voidaan tehdä myös tähän tapaan (kenties hieman elegantimmin). # Kirjoitetaan lauseke (x-x0)*(x-x1)...(x-xn) ja jaetaan se # (x-xj):llä. # Tämä antaa osoittajan. Nimittäjä saadaan tästä korvaamalla x arvolla # xj. Näin johdutaan koodiin: 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: # Tässä siis kirjoitimme yhden Lagrangen polynomin kerrallaan ja # aloitimme indeksoinnin 0:sta. # Lagrangen kertojajono muodostetaan nyt näin: # # seq(L(j,xd,x),j=0..nops(xd)-1); # ############# Kompleksianalyysi ################### # Kompleksitason käyräintegraali # cintf:=proc(f,p,a,b) # f on integroitava kompleksifunktio, määriteltävä funktioksi # p on polun määrittelevä kompleksifunktio, kuten yllä # a ja b ovat parametrivälin alku- ja loppupisteet. # local intfun,t,dp; dp:=diff(p(t),t); intfun:=evalc(f(p(t))*dp); evalf(Int(intfun,t=a..b)); end: ############ Tavalliset differntiaaliyhtälöt ##################### # #################################################################### # Eulerin menetelmä (skalaaritapaus) y'=f(t,y). y_{n+1}=y_n+h*f(t_n,y_n) #################################################################### Euler:=proc(f,a,b,ya,m) local n,h,t,y; h:=evalf((b-a)/m); t[0]:=a;y[0]:=ya; for n from 0 to m do y[n+1]:=y[n]+h*f(t[n],y[n]); t[n+1]:=t[n]+h; od; [seq([t[n],y[n]],n=0..m)]; end: #Esim. ########### #`Example`; #`f:=(t,y)->t-y^2;`; #`e3:=Euler(f,0,5,1,3);`; #`plot(e3);`; #`For more, see nslib98.mpl-file`; #Useita piirroksia: ################## #e6:=Euler(f,0,5,1,6):e12:=Euler(f,0,5,1,12): #ple3:=plot(e3):ple6:=plot(e6):ple12:=plot(e12): #with(plots): #display({ple3,ple6,ple12}); # ############ Vektoriversio ################ # EulerV:=proc(F,a,b,Ya,m) local n,h,t,Y; h := evalf((b - a)/m); t[0] := a; Y[0] := Ya; for n from 0 to m do Y[n + 1] := evalm(Y[n] + h*F(t[n],Y[n])); t[n + 1] := t[n] + h od; [seq([t[n], op(convert(Y[n],list))], n = 0 .. m)] end: # # ############################################################## # Parannettu Euler eli Heun (Improved Euler) ############################################################## Heun:=proc(f,a,b,ya,m) local n,h,t,y,ytahti; h:=evalf((b-a)/m); t[0]:=a;y[0]:=ya; for n from 0 to m do t[n+1]:=t[n]+h; ytahti:=y[n]+h*f(t[n],y[n]); y[n+1]:=y[n]+(h/2)*(f(t[n],y[n])+f(t[n+1],ytahti)); od; [seq([t[n],y[n]],n=0..m)]; end: ############ Vektoriversio ################ # HeunV:=proc(F,a,b,Ya,m) local n, h, t, Y, Ytahti; h := evalf((b-a)/m); t[0] := a; Y[0] := Ya; for n from 0 to m do t[n+1] := t[n]+h; Ytahti := evalm(Y[n]+h*F(t[n],Y[n])); Y[n+1] := evalm(Y[n]+1/2*h*(F(t[n],Y[n])+F(t[n+1],Ytahti))) od; [seq([t[n],op(convert(Y[n],list))],n = 0 .. m)] end: ########################################################### # Runge-Kutta 4 ########################################################### RK4:=proc(f,a,b,ya,m) local n,h,t,y,k1,k2,k3,k4; h:=evalf((b-a)/m); t[0]:=a;y[0]:=ya; for n from 0 to m do k1:=h*f(t[n],y[n]); k2:=h*f(t[n]+h/2,y[n]+k1/2); k3:=h*f(t[n]+h/2,y[n]+k2/2); k4:=h*f(t[n]+h,y[n]+k3); y[n+1]:=y[n]+(1/6)*(k1+2*k2+2*k3+k4); t[n+1]:=t[n]+h; od: [seq([t[n],y[n]],n=0..m)]; end: ########################################################## # Backward Euler (Implicit Euler) ########################################################## # BE with Newton iteration, works fine ! ############################################## BEN:=proc(f,a,b,Y0,N,niter) local y0,t0,y1,polyg,G,F,h,i,t1; h:=evalf((b-a)/N); y0:=Y0;t0:=a; polyg:=[a,y0]; for i from 0 to N do t1:=t0+h; F:=unapply(y-y0-h*eval(f(t1,y)),y); G:=unapply(simplify(y-F(y)/D(F)(y)),y); y1:=y0+h*eval(f(t0,y0)); # Initial value for # iteration (foreuler) y1:=(G@@niter)(y1); # Newton iteration polyg:=polyg,[t1,y1]; t0:=t1;y0:=y1; od: [polyg]; end: ############################################### # Example: # f:=(t,y)->-10*(t-1)*y; # sol:=t->exp(-5*(t-1)^2); # map(sol,[seq(1+i*.25,i=0..4)]); # [[BEN(f,1,2,1,4,5)],[Euler(f,1,2,1,4)]];map(matrix,"); # with(plots):gfor:=plot([Euler(f,1,4,1,4)],color=GREEN): # gb:=plot([BEN(f,1,4,1,4,5)],color=BLACK): # gt:=plot(sol,1..5): # display({gb,gt}); # display(gfor); #################################################### ### Matrix for solving Laplace or Poisson equ using 5-point # difference mesh lapmatr:=proc(p) local B,C,v; B:=band([1,-4,1],p); C:=diag(seq(B,i=1..p)); v:=vector(2*p+1,0); v[1]:=1;v[2*p+1]:=1; evalm(C+band(v,p^2)) end: ## Grid lines for a 2d plot ## Andersson p. 140 ## grid2D:=proc(rx,ry) local x,y; CURVES(seq([[x,op(1,ry)],[x,op(nops(ry),ry)]],x=rx), seq([[op(1,rx),y],[op(nops(rx),rx),y]],y=ry), COLOUR(RGB,.5,.5,.5) ) end: ## Uniform grid lines, example: ## grl:=ugrid(0..5,10,-3..7,15) ## ugrid:=proc(rx,nx,ry,ny) local a,b,c,d,hx,hy,i,j; a:=op(1,rx);b:=op(2,rx); c:=op(1,ry);d:=op(2,ry); hx:=(b-a)/nx;hy:=(d-c)/ny; grid2D([seq(a+i*hx,i=0..nx)],[seq(c+j*hy,j=0..ny)]) end: ###################### #Fourier-sarjat ########################## # # Jaksollinen laajennus: MapleTech Vol 3 No. 3 1996 # PeriodicExtender:=proc(f,d::range) subs({'F'=f,'L'=lhs(d),'D'=rhs(d)-lhs(d)}, proc(x::algebraic) local y; y:=floor((x-L)/D); F(x-y*D); end) end: # Esim: #sw:=PeriodicExtender(signum,-1..1); #plot(sw,-4..4); # parillinenjatko:=f->unapply(piecewise(x>0,f(x),x<0,f(-x)),x): paritonjatko:=f->unapply(piecewise(x>0,f(x),x<0,-f(-x)),x): # Esim. Ei vielä testattu, mutta vrt. työarkkiin fourier.mws # fe:=parillinenjatko(x->x);fo:=paritonjatko(x->x); # plot([PeriodicExtender(fe,-1..1),PeriodicExtender(fo,-1..1)],-5..5); # parillinen_osa:=f->unapply((1/2)*(f(x)+f(-x)),x): pariton_osa:=f->unapply((1/2)*(f(x)-f(-x)),x): # Parillinen osa on se, josta syntyy cos-sarja, pariton se, josta syntyy # sin-sarja. # #### 7.2. Vektorilaskentaa ######### # 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]); # # Yleiemmin: 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)]); # ############ Luku 10 Laplace-muunnokset ########### # # Konvoluutio # konvoluutio:=(f,g)->unapply(int(f(tau)*g(t-tau),tau=0..t),t): # # Esim: Vakiofunktion f(t)=1 ja funktion g(t)=t konvoluutio # # h:=konv(t->1,t->t); #