Numeerista integrointia, osa 2, L/numint2.mws

V2/2002, 17.2.02 HA

Alustukset

> restart:

Warning, the name changecoords has been redefined

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

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

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

Muutettiin vähän tyyliä. Otetaan sittenkin gausstaulukko parametrilistaan, eikä käytetä globaalia. Vaihdetaan nimi isolla G:llä alkavaksi (vanhakin on tallella.)

Tässä koko arsenaali:

> op(gausstaulukko);op(gaussolmut); op(gausspainot);op(Gaussint);

proc (n) options operator, arrow; [[gaussolmut(n)],...

proc (k) options operator, arrow; fsolve(orthopoly[...

proc (n) options operator, arrow; seq(int(L1(i,[gau...

proc (taulukko, f) local s, w, i; s := taulukko[1];...

> Gaussint(gausstaulukko(3),cos); #"Kertakäyttötapa" pienellä n:llä.

1.683003547

> int(cos(x),x=-1..1): evalf(%);

1.682941970

Integrointi mielivaltaisen välin yli, testausta

Otetaan mukaan muuttujanvaihto, jolla saadaan integrointi mieleivaltaisen välin [a,b] yli:

> 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));

Gaussab := proc (taulu, f, a, b) options operator, ...

Testataan:

> Gaussab(gausstaulukko(3),cos,-1,1);

1.683003547

> Gaussab(gausstaulukko(3),cos,-2,2);

1.825780686

> int(cos(x),x=-2..2): evalf(%);

1.818594854

Esim:

> int(4/(1+x^2),x=0..1);

Pi

> f:=x->4/(1+x^2);

f := proc (x) options operator, arrow; 4*1/(x^2+1) ...

> n:=2: gtaulu2:=gausstaulukko(n):n:=4: gtaulu4:=gausstaulukko(n):n:=8: gtaulu8:=gausstaulukko(n):

> Gaussab(gtaulu2,f,0,1);

3.147540983

> Gaussab(gtaulu4,f,0,1);

3.141611905

> Gaussab(gtaulu8,f,0,1);

3.141592655

Entä, jos jaetaan väli [0,1] kahteen osaan, joilla kummallakin sovelletaan 4:n pisteen Gaussia:

> Gaussab(gtaulu4,f,0,1/2)+Gaussab(gtaulu4,f,1/2,1);

3.141592699

> evalf(Pi);

3.141592654

9:nnessä numerossa on eroa (4:n verran), "globaali" on tässä esimerkissä senverran parempi. Sietää tutkia lisää.

Tasointegraali ja avaruusintegraali

Tästä aiheesta ei ole KRE-kirjassa (muista "oppikirjoistamme" puhumattakaan) mitään. Tärkeimmät lähteemme:

[KMN] Kahaner-Moler-Nash: Numerical Methods and Software

[Rec] Press-Flannery-Teukolski-Vetterling: Numerical Recipes, The Art of Scientific Computing

[BF] Burden-Faires: Numerical Analysis

Tehtävä on merkittävästi vaikeampi kuin 1 muutt, fkt.

Monte-Carlo on hyvä ehdokas.

Jos reuna on sileä ja yksinkertainen ja funktio on sileä, niin peräkköiset 1-ulotteiset integraalit toimivat.

Lukuteoreettiset menetelmät: Jokunen sana tuonnempana.

Funktion ja alueen luonteesta riippuen on paikallaan toisinaan jakaa alue eri osiin ja laskea integraalit yhteen. (Adaptiivinen algoritmi, mutta toisinaan on mahdotonta

automatisoida, tarvitaann käsityötä.

Kirjoitetaan ensin integrointi yli yksikköneliön. Muutetaan gaussintnelio iso-G-alkuiseksi (ja lyhennetään vähän nimeä). Siinä taas siirretään aiemmat globaalit parametriilistaan.

> op(Gaussnelio);

proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...
proc (gtaulux, gtauluy, f) local sx, wx, sy, wy, i,...

> gtaulux:=gausstaulukko(3);gtauluy:=gausstaulukko(4):

gtaulux := [[-.7745966692, 0., .7745966692], [.5555...

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^5*y^7);

0.

> m:=3: n:=4: d1:=2*m-1; d2:=2*n-1;

d1 := 5

d2 := 7

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^d1*y^d2);

> int(int(x^d1*y^d2,x=-1..1),y=-1..1): evalf(%);

0.

0.

Tämä ei merkitse, että mikä tahansa astetta 12 oleva kahden muuttujan polynomi integroituisi tarkasti. Jos vaikka d1:een lisätään 1 ja d2:sta vähennetään 1, niin mikään

ei enää takaa, että integroituminen olisi tarkkaa (koska x-puolella ylitettiin tarkan integroituvuuden asteluku).

Tiedetään vain, että kaikki astelukua d=min{d1,d2} olevat 2:n muuttujan polynomit integroituvat tarkkaan ja lisäksi eräät korkeampiasteiset, kuten x^d1 y^d2

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^(d1+1)*y^(d2-1)),int(int(x^(d1+1)*y^(d2-1),x=-1..1),y=-1..1): evalf(%);

.6857142854e-1, .8163265306e-1

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^4*y^4),int(int(x^4*y^4,x=-1..1),y=-1..1): evalf(%);

.1599999999, .1600000000

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^5*y^6),int(int(x^5*y^5,x=-1..1),y=-1..1): evalf(%);

0., 0.

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^6*y^6),int(int(x^6*y^6,x=-1..1),y=-1..1): evalf(%);

.6857142854e-1, .8163265306e-1

> Gaussnelio(gtaulux,gtauluy,(x,y)->x^6*y^6),int(int(x^6*y^6,x=-1..1),y=-1..1): evalf(%);

.6857142854e-1, .8163265306e-1

Nyt voitaisiin tehdä muunnos integroitavalta alueelta yksikköneliölle ja palauttaa integraali sillä tavoin tähän. Ehkäpä otetaan harjoituksen aiheeksi.

Seuraavassa suoritetaan peräkkäiset integroinnit ensin x-projisioituvalle tasoaluleelle ja sitten xy-projisioituvalle avaruuskappaleelle.

Siinä emme käytä Gaussnelio-funktiota, vaan rakennamme Gaussab :n varaan.

Tasointegraali peräkkäisinä integraaleina

> a2:=x->-sqrt(1-x^2);b2:=-a2;

a2 := proc (x) options operator, arrow; -sqrt(1-x^2...

b2 := -a2

> gtaulux;f:=(x,y)->1;

[[-.7745966692, 0., .7745966692], [.5555555555, .88...

f := 1

> F:=x->Gaussab(gtauluy,f,a2(x),b2(x));

F := proc (x) options operator, arrow; Gaussab(gtau...

> F(0.5);

1.732050808

> Gaussab(gtaulux,F,-1,1);

3.183234516

Mahtava juttu, nyt Maple tottelee ihan normaalia matematiikkaa!

> Gaussxy:=(gtaulux,gtauluy,f,a1,b1,a2,b2)->Gaussab(gtaulux,x->Gaussab(gtauluy,y->f(x,y),a2(x),b2(x)),a1,b1);

Gaussxy := proc (gtaulux, gtauluy, f, a1, b1, a2, b...

> gtaulux:=gausstaulukko(3):

> Gaussxy(gtaulux,gtaulux,(x,y)->x^2*y^2,-1,1,a2,b2);

.1124365390

> n:=5: gtaulux:=gausstaulukko(n): Gaussxy(gtaulux,gtaulux,1,-1,1,x->-sqrt(1-x^2),x->sqrt(1-x^2));

3.151812671

> for n from 3 to 15 do gtaulux:=gausstaulukko(n): ala[n]:=Gaussxy(gtaulux,gtaulux,1,-1,1,x->-sqrt(1-x^2),x->sqrt(1-x^2)): od:

>

> seq(ala[n],n=3..15);

3.183234516, 3.160555054, 3.151812671, 3.147727965,...
3.183234516, 3.160555054, 3.151812671, 3.147727965,...

Avaruusintegraali peräkkäisinä integraaleina

Nyt homma osataan. Kyse on vain siitä, tehdäänkö pitkä "one-liner" vai pari lyhyempää riviä.

Laskettavana on

> restart: Int(Int(Int(f(x,y,z),z=a3(x,y)..b3(x,y)),y=a2(x)..b2(x)),x=a1..b1);

Warning, the name changecoords has been redefined

Int(Int(Int(f(x,y,z),z = a3(x,y) .. b3(x,y)),y = a2...

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

>

> # G:=(x,y)->Gaussab(gtaulux,z->f(x,y,z),a3(x,y),a2(x,y)); # Tämä on kehittelyä.

> 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);

Gaussxyz := proc (gtaulux, gtauluy, gtauluz, f, a1,...
Gaussxyz := proc (gtaulux, gtauluy, gtauluz, f, a1,...

Houkutus one-lineriin oli liian suuri!

> a3:=(x,y)->x-y: b3:=(x,y)->x+y:

> Int(Int(Int(z,z=a3(x,y)..b3(x,y)),y=x^2..x),x=0..1): %=evalf(value(%));

Int(Int(Int(z,z = x-y .. x+y),y = x^2 .. x),x = 0 ....

> gtaulux:=gausstaulukko(5):

>

> Gaussxyz(gtaulux,gtaulux,gtaulux,(x,y,z)->z,0,1,x->x^2,x->x,a3,b3);

.8333333343e-1

>

Eipä huono ensimmäiseksi testiksi!