# Numeerista integrointia, L/numint.mws # V2/2002, viimeinen viikko HA # Alustukset > #read("c:\\opetus\\maple\\v202.mpl"): > read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl"): > #op(L); > op(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 proc > seq(L1(i,[x1,x2,x3],x),i=1..3); (x - x2) (x - x3) (x - x1) (x - x3) (x - x1) (x - x2) -------------------, -------------------, ------------------- (x1 - x2) (x1 - x3) (x2 - x1) (x2 - x3) (x3 - x1) (x3 - x2) # Kertausta: Lagrangen interpolaatio # Lagrangen interpolaatiomenetelmässä muodostetaan annettuun xdataan # [x0,...,xn] liittyvät polynomit L[0],...,L[n], siten että # L[i](x[j])=delta[i,j] > L:='L': L[i](x[j])=delta[i,j]; > read("c:\\opetus\\maple\\v202.mpl"): > xd:=[x0,x1,x2,x3]; > L(0,xd,x); > L(1,xd,x); # jne. Yleisesti: KRE s. 849 Lagrange interpolation. # Lagrangen (kertoja)polynomit on valmiina koodattu tiedostossamme # v202.mpl, joka luettiin yllä. Siksi nuo kaavatkin tulivat # yllä ihan oikein, kun vaan tarjottiin L-kirjainta! # Lagrangen polynomien avulla voidaan kirjoittaa suoraan yleinen # interpolaatiotehtävän ratkaisu: # Tehtävä: # Annettu xdata=[x0,...,.xn] ja ydata=[y0,..,yn]. Määrättävä # korkeintaan astetta n oleva polynomi näiden datapisteiden # kautta. # Tiedämme, että polynomin kertoimet saataisiin ratkaisemalla # lineaarinen yhtälösysteemi, jonka kerroinmatiirina on Vandermonden # matriiisi. Kuten PNS:n yhteydestä muistamme, Vandermonde on # ei-singulaarinen ja tässä tapauksessa # neliömatriisi, joten 1-käs ratkaisu on. # Tätä tietoa emme tarvitse Lagrangen menetelmässä, emme myöskään # joudu ratkaisemaan yhtälösysteemiä, vaan nerokkuus # piilee Lagrangen polynomien ortogonaalisuustyyppisessä # ominaisuudessa (kts. yllä). # Ratkaisu: > L:='L': > 'p(x)'=sum(y[i]*L[i](x),i=0..n); # Jäännöstermille voidaan johtaa lauseke: > R=product((x-x[i]),i=0..n)*D[n+1](f)(xi)/n!; # Tällöin käytettävissä on funktio (n+1) kertaa derivoituva funktio f, # jonka arvot x-datapisteissä ovat y-data-arvoja: > y[i]=f(x[i]), i=0..n; > # Kertausta : Gaussin integrointi 1-dim tapaus # Erilaisia integrointisääntöjä "quadrarure rules" saadaan integroimalla # interpolaatiopolynomeja. # Ensimmäiseksi mieleen juolahtava ajatus on jakaa väli tasavälisesti ja # muodostaa interpolaatiopolynomi, joka integroidaan. # Näin saadaan ns. Newton-Cotes-kaavat. Kaksi alinta astetta olevaa # Newton-Cotes-kaavaa ovat trapetsisääntö (puolisuunnikas-) ja # Simsonin sääntö. Niissä approksimidaan pinta-alaa puolisuunnikkaalla # ja vastaavasti 2. adteen polynomin reunustamalla alalla. # # Yhdistetyt säännöt: ("Composite rules") # Alhaista kertalukua olevilla kaavoilla laskettaessa jaetaan koko # integroimisalue osiin, joilla kullakin sovelletaan ao. sääntöä. Koko # integraali on osien summa. # # Optimaalinen solmujen valinta, Gaussin säännöt. # # Yleisemmin voidaan valita mikä tahansa (ei välttämättä tasavälinen) # pisteistö x1,...,xn integrointiväliltä [a,b]. # Huom! Vaihdoimme indeksoinnin alkamaan 1:stä, tämä korostaa sitä, että # johtamamisen alaisena olevat kaavat ovat yleensä ns. "avoimia # kaavoja", # joissa solmupisteet x[i] ovat välin sisäpisteitä. # > L:='L':p:=sum(y[i]*L[i](x),i=1..n); n ----- \ p := ) y[i] L[i](x) / ----- i = 1 > intsaanto:=sum(f(x[i])*Int(L[i](x),x=a..b),i=1..n); n b ----- / \ | intsaanto := ) f(x[i]) | L[i](x) dx / | ----- / i = 1 a # Merk: > w[i]=Int(L[i](x),x=a..b); b / | w[i] = | L[i](x) dx | / a > intsaanto:=sum(w[i]*f(x[i]),i=1..n); n ----- \ intsaanto := ) w[i] f(x[i]) / ----- i = 1 # Yleisesti saadaan painotettu summa funktion arvoista. # Gaussin idea: Annetaan solmujen x[i] olla muuttujina # optimointitehtävässä: # Valitse solmut x[i] ja painot w[i] siten, että # integrointisääntö on tarkka mahdollisimman korkea-asteisilla # polynomeilla. # Korkein asteluku, joka voidaan toivoa saavutettavaksi, on 2n-1 # (2n tuntematonta, 2n määrättävää kerrointa). # # Ratkaisutapa 1: Muodostetaan (epälineaarinen) yhtälösysteemi # integroimalla monomeja 1,x,x^2, ... . # Tällä saadaan alhaista kertalukua olevia (ainakin n=2). Korkeampia # voidaan yrittää vaikka Newtonin menetelmällä # ratkaista, mutta se ei johda pitkälle. # # Ratkaisutapa 2: Normeerataan integrointi välille [-1,1]. # Osoittautuu, että solmut x1,...,xn saadaan ns. Legendren # polynomien nollakohtina. Kun solmut tunnetaan, saadaan painot # yksinkertaisesti kaavasta: # > w[i]=Int(L[i](x),x=a..b); b / | w[i] = | L[i](x) dx | / a # Integrointi yleisen välin [a,b] yli hoituu yksinkertaisella # lineaarisella muuttujanvaihdolla. (Nokkelia kun ollaan, # niin muodostetaan 1. asteen interpolaatiopolynomi, joka a:ssa saa # arvon -1 ja b:ssä arvon 1.) > read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl"): > (-1)*L(0,[a,b],x)+1*L(1,[a,b],x);simplify(%); x - b x - a - ----- + ----- a - b b - a -2 x + b + a ------------ a - b > t:='t': > tyht:=t=%%; -2 x + b + a tyht := t = ------------ a - b > xx:=solve(%,x); xx := - 1/2 t a + 1/2 t b + 1/2 b + 1/2 a > dx:=diff(xx,t); dx := - 1/2 a + 1/2 b # Saadaan integrointikaava: > Int(f(x),x=a..b)=Int(f(xx)*dx,x=-1..1); b 1 / / | | | f(x) dx = | | | / / a -1 f(- 1/2 t a + 1/2 t b + 1/2 b + 1/2 a) (- 1/2 a + 1/2 b) dx # Maple-toteutus > with(orthopoly); Warning, the name L has been redefined [G, H, L, P, T, U] > #read("c:\\opetus\\maple\\v202.mpl"): > read("/home/apiola/opetus/peruskurssi/v2-3/maple/v202.mpl"): # > Nimikonflikti on lähellä, emme kuitenkaan tarvitse polynomia L - > Laguerren polynomi, vaan Legenren polynomia, jonka Maplen orthopoly > nimeää P:ksi. # ( L-matemaatikoita riittää: Lagrange, Legendre, Laguerre, Laplace, ... # ) # Gaussin integroinnin yhteydessä on tapana indeksoida Lagrangen # polynomit (siis solmut) 1..n. Siis datapisteitä on # n kpl, jolloin interpoloidaan n-1 asteisella polynomilla. Tavoitteena # siis kaava, joka on tarkka asteeseen 2n-1 saakka. # Jatkon indeksipeliä helpottaa, kun kirjoitamme Lagrangen # kertojapolynomista version L1, jossa indeksointi menee # 1..n. (Koodikin on hiukan yksinkertaisempi.) > op(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 proc > > P(10,x); 46189 10 109395 8 45045 6 15015 4 3465 2 63 ----- x - ------ x + ----- x - ----- x + ---- x - --- 256 256 128 128 256 256 > Digits:=15: > fsolve(P(10,x)=0,x); -.973906528517172, -.865063366688985, -.679409568299024, -.433395394129247, -.148874338981631, .148874338981631, .433395394129247, .679409568299024, .865063366688985, .973906528517172 > Digits:=10: > gausssolmut:=seq([fsolve(P(k,x)=0,x)],k=1..5); gausssolmut := [0.], [-.5773502692, .5773502692], [-.7745966692, 0., .7745966692], [-.8611363116, -.3399810436, .3399810436, .8611363116], [-.9061798459, -.5384693101, 0., .5384693101, .9061798459] > %[2,1]; > gaussolmut:=k->fsolve(P(k,x)=0,x); gaussolmut := k -> fsolve(P(k, x) = 0, x) > L1(3,[a,b,c],x); (x - a) (x - b) --------------- (c - a) (c - b) > gausspainot:=n->seq(int(L1(i,[gaussolmut(n)],x),x=-1..1),i=1..n); gausspainot := 1 / | n -> seq( | L1(i, [gaussolmut(n)], x) dx, i = 1 .. n) | / -1 > gausspainot(3); .5555555555, .8888888889, .5555555555 > gausstaulukko:=n->[[gaussolmut(n)],[gausspainot(n)]]; gausstaulukko := n -> [[gaussolmut(n)], [gausspainot(n)]] > gausstaulukko(20); [[-.9931285992, -.9639719273, -.9122344283, -.8391169718, -.7463319065, -.6360536807, -.5108670020, -.3737060887, -.2277858511, -.07652652113, .07652652113, .2277858511, .3737060887, .5108670020, .6360536807, .7463319065, .8391169718, .9122344283, .9639719273, .9931285992], [ .01761400712, .04060142981, .06267204838, .08327674164, .1019301199, .1181945321, .1316886385, .1420961094, .1491729865, .1527533872, .1527533872, .1491729866, .1420961096, .1316886385, .1181945321, .1019301199, .08327674164, .06267204838, .04060142977, .01761400713]] > sw:=gausstaulukko(4); sw := [[-.8611363116, -.3399810436, .3399810436, .8611363116], [.3478548453, .6521451549, .6521451549, .3478548450]] > 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: > > gtaulu:=gausstaulukko(3); gtaulu := [[-.7745966692, 0., .7745966692], [.5555555555, .8888888889, .5555555555]] > gaussint(3,x->x^4),int(x^4,x=-1..1); evalf(%); .3999999998, 2/5 .3999999998, .4000000000 > n:=3: gtaulu:=gausstaulukko(n): > seq([gaussint(3,x->x^k),evalf(int(x^k,x=-1..1))],k=1..2*n-1); [0., 0.], [.6666666664, .6666666667], [0., 0.], [.3999999998, .4000000000], [0., 0.] > > display(nelio,plot([ss],style=point)); # Integrointi 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); gtaulux := [[-.7745966692, 0., .7745966692], [.5555555555, .8888888889, .5555555555]] > gtauluy:=gtaulux: > gaussintnelio(3,3,(x,y)->x^4*y^4),int(int(x^4*y^4,x=-1..1),y=-1..1); > evalf(%); .1599999998, 4/25 .1599999998, .1600000000 > sx:=gtaulux[1]:sy:=gtauluy[1]: > with(plottools): with(plots): > n:=3: > nelio:=polygon([[1,1],[-1,1],[-1,-1],[1,-1]]): > ss:=seq(seq([sx[i],sy[j]],j=1..n),i=1..n): > display(nelio,plot([ss],style=point),axes=frame); > # #