H/harj12ohje.mws

25.4.02 HA

Nämä viimeiset LV harjoitukset ovat kenties pitkät, mutta myös hyvin (=liian?) pitkälle neuvotut tällä työarkilla.

Ratkaisussa tulee olla selvillä siitä, mitä ENTER-painallukset tarkoittavat. Toki on jätetty sinne tänne täydennettäviä aukkoja.

Itselläni oli aika paljon työtä, mutta asiat ovat niin valmiiksi pureskeltuja, että rastien pitäisi olla lopulta aika kevyitä.

Alustukset

> restart:

> #read("c:\\opetus\\maple\\v202.mpl"): # HA:n kotikone (Win)

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

> #read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl"): # HA:n Linux

> # Myös voit ottaa: www../maple/v202.mpl - tieodoston ja sijoittaa sen sopivan polun varteen. Siinä tapauksessa

> # kannattaa päivittää, koska v202.mpl elää.

> with(plots):

Huom! Ennenkuin suoritat komentoja, poista read-lauseista soveltuvasti kommentit ja tarpeen mukaan muuta polkua.

1.

Kun käytettävissä on harj12av.mws, niin opettavaisinta lienee laskea ensin aivan samalla tavalla ottamalla vain

isommat gausstaulukot käyttöön. Samalla saamme hyödyn AV teht. 5:ssä laskemastamme muunnoksesta.

> restart:

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

> gaussmn:=Sum(Sum(wx[i]*wy[j]*intf(sx[i],sy[j]),j=1..n),i=1..m);

Tässä on yleinen tulokaava. Lasketaan nyt AV teht. 5:n tilanne

Katso funktioiden gaussolmut ja gausspainot sisältö:

> op(gaussolmut); op(gausspainot);

> m:=5: sx:=gaussolmut(m);wx:=gausspainot(m);

> n:=4: sy:=gaussolmut(n):wy:=gausspainot(n):

> f:=(x,y)->ln(x+2*y);

> x:=... y:= ... # se lineaarinen muunnos + siirto

> dxdy_per_dudv:= ... # pinta-alkioiden muuntosuhde

> F:=f(x,y); # Integroitava funktio u:n ja v:n avulla.

Muunsimme integroinnin yksikköneliön yli (kuten AV 5:ssä).

> intf:=unapply(F*dxdy_per_dudv,u,v); # Kerrotaan vielä muuttujanvaihdoksen muuntosuhteella.

> value(gaussmn);

> int(int(f(xx,yy),xx=1.4..2),yy=1..1.5);

> f_arvoja=m*n;

Testaile !

Kokeillaan samalla myös valmiita funktioitamme

> gtaulux:=gausstaulukko(5): gtauluy:=gausstaulukko(4):

> Gaussnelio(gtaulux,gtauluy,intf);

> Gaussnelio(gausstaulukko(3),gausstaulukko(3),intf); # Näin ei tehdä käytönnössä, mutta "teoriassa" kätevää.

> with(LinearAlgebra):

> Matrix(3,3,(m,n)->Gaussnelio(gausstaulukko(m),gausstaulukko(n),intf));

> Matrix(5,5,(m,n)->Gaussnelio(gausstaulukko(m),gausstaulukko(n),intf));

Tämä on tosi eleganttia, kylläkin tuhlailevaa laskentaa. Lähes sama eleganssi voitaisiin saavuttaa laskemalla ensin valmiit taulukot ja sitten hakemalla niistä. Jos tehtäisiin isompi taulukko, niin sitten kannattaisi jossain vaiheessa.

Gaussin kaavoja käytetään elävässä elämässä tietysti niin, että muodostetaan ensin sopivan kokoinen taulukko ja sitten

haetaan aina sieltä. Siksi funktiomme onkin rakennettu niin, että taulukko annetaan argumenttina. (Taulukko voitaisiin vaikka

lukea tiedostosta.)

Gaussxy :n kokeilu:

Tämä on sikäli helppokäyttöisempi, että ei tarvita muuttujan vaihtoa:

> map(nops,gtaulux); map(nops,gtauluy);

> op(f);

> Gaussxy(gtaulux,gtauluy,f,1.4,2,1,1.5);

> Gaussnelio(gtaulux,gtauluy,intf);

Kylläpä toimivat kauniisti, vai mitä! (Toivottavasti jäi jotain tehtävää.)

2.

Tässä on tilaisuus harjoitella muunnoksia alueilta toisille.

> restart:with(plots): with(plottools): 

Warning, the name changecoords has been redefined

> kolmio:=[[0,0],[1,0],[1,1]]:kolmiokuva:=polygon(kolmio,filled=true,color=yellow):

> display(kolmiokuva);"xy-taso";

> nelio01:=[[0,0],[0,1],[1,1],[1,0]]: # [0,1] x [0,1]

> nelio01kuva:=polygon(nelio01,filled=true,color=green):

> display(nelio01kuva);"uv-taso";

> nelio_11:=[[1,1],[-1,1],[-1,-1],[1,-1]]: # [-1,1] x [-1,1]

> nelio_11kuva:=polygon(nelio_11,filled=true,color=cyan):

> display(nelio_11kuva); "zw-taso";

Olkoot:

> g:=(u,v)->[u,u*v]; gv:=uv->g(uv[1],uv[2]); # Kuvaus [0,1] x [0,1] --> Delta (kolmio)

> h:=(w,z)->[..,..]; hv:=uv->h(uv[1],uv[2]); # Kuvaus [-1,1] x [-1,1] --> [0,1] x [0,1] .

Aivan samoin kuin optimoinnin yhteydessä (Muistatko SD:n?) on kätevää määritellä kahden muuttujan funktioista myös lista/vektori-versiot gv ja hv .
Niiden avulla on vaivatonta katsoa, miten eri pisteet kuvautuvat.

> map(hv,nelio_11);map(gv,%);

Näinhän niiden pitikin mennä.

Katsotaan Gaussin pisteitä. Voi olla, että olemme käyttäneet restarttia, no tiedämme, mitä pitää olla 2:n pisteen tapauksessa:

> s[1]:=-1/sqrt(3); s[2]:=-s[1];

> ss:=[seq(seq([s[i],s[j]],j=1..2),i=1..2)];

> gpnelio01:=map(hv,ss);

> display(nelio01kuva,plot(gpnelio01,style=point,symbol=cross)); "uv-taso";

> gpkolmio:=map(gv@hv,ss); # Näin kuvautuvat neliön [-1,1] x [-1,1] Gauss2-pisteet kolmioon Delta.

> display(kolmiokuva,plot(gpkolmio,style=point,symbol=cross));

Gaussin pisteiden sijainti kolmiossa ei näytä optimaaliselta. Oikeassa laidassa pisteet ovat paljon harvemmassa.
Miettimisen paikka: Mitä kannattaisi tehdä. Jaetaan alue sopivasti osiin. Ehkä kannattaisi ottaa suht. iso N vaakasuunnassa ja kutakin
pystyjanaa kohti kasvatettaisiin N:ää lineaarisesti pystysuunnassa. Näillä Gauss-palikoilla olisi helppo toteuttaa.
Parempi vielä: Siirrytään käyttämään ns. pinta-alakoordinaatteja: Kolmion sisäpiste yhdistetään janalla kuhunkin kärkeen ja otetaan syntyneiden kolmioiden
pinta-alat pisteen koordinaateiksi.
FEM-laskennassa suoritetaan usein integrointia kolmion yli, Pisteet valitaan usein reunalta, mediaanien leikkauspisteestä jne., riippuen kantafunktiotyypistä.
Periaatteessa tehdään samanlainen mahd. korkean asteluvun vaatimus.
Tutkimuksen ja opiskelun aihe ... FEM-kirjallisuus ... (FEM = Finite Element Method)
Testataan integraalikaavoja esimerkeillä.

> f:='f':Ixy:=Int(Int(f(x,y),y=0..x),x=0..1);

> Iuv:=Int(Int(f(u,u*v)*u,v=0..1),u=0..1);

> Iwz:=(..)*Int(Int(f((..,..)*(),z=-1..1),w=-1..1);

> f:=(x,y)->x*y; fv:=xvek->f(xvek[1],xvek[2]);

> value(Ixy);

> value(Iuv);

> value(Iwz);

Tehdään nyt varsinainen tämän tehtävän esimerkki.

> f:=(x,y)->exp(-x^2*y^2); fv:=xvek->f(xvek[1],xvek[2]);

> uv:=h(w,z); xy:=gv(uv);

> intfunwz:=...;

> Int(Int(intfunwz,z=-1..1),w=-1..1): %=value(%);evalf(rhs(%));

> x:='x': y:='y':Int(Int(f(x,y),y=0..x),x=0..1): %=value(%); evalf(rhs(%));

Saasaan samalla näyttö Maplen taidoista laskea erikoisfunktioilla (erf ja hypergeom).

No nyt sitten Gaussiin.

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

> gtaulux:=gausstaulukko(5);

> gtauluy:=gausstaulukko(4);

> fzw:=unapply(intfunwz,w,z);

> op(Gaussnelio);

> Gaussnelio(gtaulux,gtauluy,fzw);

Otetaan vähemmän pisteitä

> gtaulu2:=gausstaulukko(2):gtaulu3:=gausstaulukko(3):gtaulu4:=gausstaulukko(4):

> Gaussnelio(gtaulu2,gtaulu3,fzw);

> Gaussnelio(gtaulu3,gtaulu2,fzw);

3.

Lieriökoordinaatteihin siirtyminen auttaa merkittävästi laskennassa, kuten numint3.mws :ssä näemme.

Lasketaan kartion z=r ja tason z=2 rajoittaman kappaleen massa ja momentti xy-tason suhteen, kun kartion tiheys rho(x,y,z) = sqrt(x^2+y^2) .

Edellinen saadaan copy/paste:lla:

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

> M:=Int(Int(Int(..,z=..),r=..),Theta=..);

> Mxy:=Int(Int(Int(..,z=..),r=..),Theta=..);

> gtaulu:=gausstaulukko(5):

> Mg:=Gaussxyz(gtaulu,gtaulu,gtaulu,(Theta,r,z)->r^2,0,2*Pi,0,2,(Theta,r)->r,2);

> value(M);evalf(%);

> Mgxy:=Gaussxyz(gtaulu,gtaulu,gtaulu,(Theta,r,z)->z*r^2,0,2*Pi,0,2,(Theta,r)->r,2);

> value(Mxy);evalf(%);

> z0g:=.. ; #Gaussilla laskettu

> z0:=.. ;evalf(%); # Tarkka intergaali ja sen pyöristys.

Tämä saa riittää, x0 ja y0 ovat ilman muuta nollia.

4.

> restart:with(linalg):

Torus:

> x:=(a+w*cos(v))*cos(u);
y:=(a+w*cos(v))*sin(u);
z:=w*sin(v);

> Jac:=jacobian([x,y,z],[u,v,w]);

> Jdet:=simplify(det(%));

> M:=Int(Int(Int(),..),..);

> M:=value(%);

Keskiön suhteen kiinnostava on x0. Symmetriasyistä on oltava y0=0, z0=0.

>

> Mxy:=Int(Int(Int(..)..)..);

> value(Mxy);

> x0:=..

>

5.

> restart:with(plots): with(plottools):

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

> T:=[[-1,1],[-1,1]];

> satu:=rand(0..10^10)/10.^10;

Kun minä ystävällisesti annan sopivat kutsut sellaisenaan, niin selvitä sinä ystävällisesti, miksi kutsu on juuri tällainen. (Pelkkä aivoton ENTER aiheuttaa minulle allergisen reaktion.)

> N:=10: [MonteCarlo2d(1,(x,y)->x^2+y^2,T,N^2),Riemann2d(1,(x,y)->x^2+y^2,T,N)]; N^2*f_arvoa;

> N:=50: [MonteCarlo2d(1,(x,y)->x^2+y^2,T,N^2),Riemann2d(1,(x,y)->x^2+y^2,T,N)]; N^2*f_arvoa; # Talleta ennenkuin kutsut kovin suurella N:llä (älä tuhlaa koko yhteisön CPU-aikaa).

Joskus voi olla järkevämpää laskea suht. pienellä N:llä ja muodostaa keskiarvo saaduista tuloksista. Sekin käy kätevästi:

> N:=20: lkm:=20: keskiarvo:=add(MonteCarlo2d(1,(x,y)->x^2+y^2,T,N^2),i=1..lkm)/lkm;lkm*N^2*f_arvoa;

> N:=30: lkm:=10: keskiarvo:=add(MonteCarlo2d(1,(x,y)->x^2+y^2,T,N^2),i=1..lkm)/lkm;lkm*N^2*f_arvoa;

Monte-Carlo näyttää olevan joka kerta parempi kuin Riemann, vai miltä sinusta näyttää. Ero näkyy erityisen selvästi pienillä laskentamäärillä. "Hyvällä onnella" saadaan

Monte Carlolla varsin kohtuullinen likiarvo hyvinkin pienellä laskentamäärällä. Toisaalta Monte Carlo saattaa isollakin laskennalla antaa yllättävästi huonomman tulloksen kuin

pienemmällä. Kannattaako siis uhkapeli?

6.

> restart:

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

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

Toruspinta:

> x:=(a+b*cos(v))*cos(u):
y:=(a+b*cos(v))*sin(u):
z:=b*sin(v):

Toki on luonnollista laskea toruskoordinaateissa, kuten teht. 4. Tyypillinen Monte Carlo olisi sellainen, jossa reunat olisivat hankalammat, kuten luentoesimerkissä.

Tällöin olisi luontevaa laskea suorak. koordinaateissa. (Samaten tehtävässä, joka olisi suoraan annettu suorak. koord.) Sitäpaitsi särmiön yli integroinnissa ei ole

kovin paljon vitsiä Monte Carlon suhteen.

Simuloidaan todellisuutta ja lasketaan suorakulmaisissa.

Muunnetaan torus siis suorakulmaisiin koordinaatteihin. (Tässä on tilanne, jossa on helpompi johtaa pinnan (ja kappaleen) esitys muunnetussa koordinaatistossa,

toki sama tilanne esiintyy pallolla pallokoordinaatistossa ja lieriöllä lieriökoordinaatistossa.)

Tämä tapahtuu etupäässä käsin. Käytämme Maplea tässä kohden lähinnä matemaattiseen tekstinkäsittelyyn, toki hyödynnämme

hänen sievennystaitojaan.

Mieti seuraavat vaiheet mielellään kynää ja paperia käyttäen, niin ne on tehtykin.

> x2y2z2:=simplify(x^2+y^2)+z^2;

x2y2z2 := a^2+2*a*b*cos(v)+b^2*cos(v)^2+b^2*sin(v)^...

> x:='x': y:='y': z:='z':

> x2y2z2:=a^2+2*a*(sqrt(x^2+y^2)-a)+b^2*cos(v)^2+b^2*sin(v)^2;

>

x2y2z2 := a^2+2*a*(sqrt(x^2+y^2)-a)+b^2*cos(v)^2+b^...

> expand(x2y2z2);

-a^2+2*a*sqrt(x^2+y^2)+b^2*cos(v)^2+b^2*sin(v)^2

> x^2+y^2+z^2+a^2=expand(2*a*(sqrt(x^2+y^2)-a))+b^2;

x^2+y^2+z^2+a^2 = 2*a*sqrt(x^2+y^2)-2*a^2+b^2

> x^2+y^2+z^2+a^2-2*a*(sqrt(x^2+y^2))=b^2;

x^2+y^2+z^2+a^2-2*a*sqrt(x^2+y^2) = b^2

Kun ajatellaan

> r=sqrt(x^2+y^2);

r = sqrt(x^2+y^2)

nähdään vasemmalla binomin neliö. Siten saadaan:

> (sqrt(x^2+y^2)-a)^2+z^2=b^2;

(sqrt(x^2+y^2)-a)^2+z^2 = b^2

Toruksen sisusta on siten:

> (sqrt(x^2+y^2)-a)^2+z^2<=b^2;

(sqrt(x^2+y^2)-a)^2+z^2 <= b^2

> x:=(a+b*cos(v))*cos(u);
y:=(a+b*cos(v))*sin(u);
z:=b*sin(v);

> a:=2: b:=1:

> plot3d([x,y,z],u=-Pi/2..Pi/2,v=0..2*Pi,scaling=constrained,labels=["x","y","z"]);

Kuvassa on suoraan näkyvissä ympäröivä laatikko (etenkin, kun kääntelet sitä).

> T:=[[..]],[..],[..]];

> satu:=rand(0..10^9)/10.^9;

> g:=(x,y,z)->(sqrt(x^2+y^2)-a)^2+z^2;

> MonteCarlo3d(1,g,T,10),evalf(a*b^2*Pi^2); # Mitä kaikkea tässä tapahtuu.

> MonteCarlo3d(1,g,T,20),evalf(a*b^2*Pi^2);

Nyt voit vielä omatoimisesti laskea keskiön, jos jaksat.