Numsym Harj. 4
18.3.04 ratkaisuja
1.
> | restart: |
Warning, the name changecoords has been redefined
> | read("/p/edu/mat-1.192/ns04.mpl"); |
> | #read("/home/apiola/opetus/numsym/04/maple/ns04.mpl"); |
> | #read("c:\\usr\\heikki\\numsym04\\maple\\ns04.mpl"); |
> | op(ode23step); |
> |
y'=1
> | f:=(t,y)->1; |
> | dsolve({diff(y(t),t)=f(t,y(t)),y(0)=0}); |
> | ytarkka:=unapply(subs(%,y(t)),t); |
> | ode23step(f,t[n],ytarkka(t[n]),h); |
> | Y[n+1]:=%[2]; |
> | ytarkka(t[n]+h); |
Samat ovat, joten lokaali virhe = 0
y'=t
> | restart: |
Warning, the name changecoords has been redefined
> | #read("/home/apiola/opetus/numsym/04/maple/ns04.mpl"); #read("c:\\usr\\heikki\\numsym04\\maple\\ns04.mpl"); read("/p/edu/mat-1.192/ns04.mpl"); |
> | f:=(t,y)->t; |
> | dsolve({diff(y(t),t)=f(t,y(t)),y(0)=0}); |
> | ytarkka:=unapply(subs(%,y(t)),t); |
> | ode23step(f,t[n],ytarkka(t[n]),h); |
> | Y[n+1]:=%[2]; |
> | subs(h=t[n+1]-t[n],%); |
> | simplify(%); #factor(%); |
> | ytarkka(tn+h); |
Yhtyvät, joten lokaali virhe on taas = 0.
y'=t^2
> | f:=(t,y)->t^2; |
> | dsolve({diff(y(t),t)=f(t,y(t)),y(0)=0}); |
> | ytarkka:=unapply(subs(%,y(t)),t); |
> | ode23step(f,t[n],ytarkka(t[n]),h); |
> | Y[n+1]:=%[2]; |
> | subs(h=t[n+1]-t[n],%); |
> | simplify(%); |
> | ytarkka(tn+h); |
Uskomatonta, mutta totta!
> | f:=(t,y)->t^3; |
> | dsolve({diff(y(t),t)=f(t,y(t)),y(0)=0}); |
> | ytarkka:=unapply(subs(%,y(t)),t); |
> | ode23step(f,t[n],ytarkka(t[n]),h); |
> | Y[n+1]:=%[2]; |
> | virhearvio:=%%[3]; |
> |
> | ytarkka(t[n]+h); |
> | expand(%); |
> | LKV[n+1]:=%-Y[n+1]; |
> | LKV[n+1]:=simplify(%); |
> | simplify(virhearvio); |
> | expand(%); |
Hiukan oudolta näyttää, että virhearvio antaa O(h^3).
Selitys on siinä, että Virhe-estimaattorina on 3:nnen ja toisen kertaluvun menetelmien tulosten erotus.
Se on siten tarkka funktiolle f(t,y)=t^2, mutta ei funktiolle t^3. Virhe-estimaattori arvioi alemman kertaluvun menetelmän virhettä,
mutta menetelmä laskee ylemmällä (3;lla). Menetelmä on siten tarkka funktiolle f(t,y)=t^3, mutta se ei "tiedä" olevansa niin hyvä.
2.
> | dy:=diff(y(x),x)=2/sqrt(Pi)*exp(-x^2); |
> | dsolve({dy,y(0)=0},y(x)); |
> | numratk:=dsolve({dy,y(0)=0},y(x),numeric); |
> | [seq([op(numratk(k*0.2)),"erf(x)=",erf(k*0.2)],k=0..10)]; |
> |
> | numy:=u->subs(numratk(u),y(x)); |
Numeerisen ratkaisun jälkikäsittely, kts. HAM s. 175 (myös ...maple/ns04.mpl)
> | numy(1)-erf(1.); |
> | plot([numy,erf+0.01],0..3); |
> | xk:=k*h: h:=0.2: |
> | array([seq([xk,numy(xk),evalf(erf(xk))],k=0..20)]); |
> |
> |
3.
> |
4.
> | restart: |
Warning, the name changecoords has been redefined
Taylorin menetelmä yhtälölle y'=f(t,y) tarkoittaa, että muodostetaan ratkaisufunktion y(t) derivaattoja pisteessä t ja käytetään niiden avulla laskettavaa, pisteessä t kehitettyä n:nnen asteen Taylorin polynomia ratkaisufunkton approksimoimiseen pisteessä t+h. Jos n=1, on kyseessä Eulerin menetelmä.
Koska y'(t)=f(t,y(t)), saadaan derivaatat muodostetuksi ketjusäännön avulla. Jos Maplelle annetaan määrittelemättömiä funktio(lausekke)ita
y(t) ym. sisältävä lauseke ja pyydetään derivoimaan, se soveltaa normaaleja derovoimissääntöjä (kuten ketjusääntöä).
Suoritetaan ensi ideointia.
> | y1:=f(t,y(t)); |
> | y2:=diff(y1,t); |
> | y2:=subs(diff(y(t),t)=y1,%); |
> | y3:=diff(y2,t); |
> | y3:=subs(diff(y(t),t)=y1,%); |
Näin voidaan jatkaa. Yleisiä kaavoja ei kannata rakentaa pitkälle, vaan otetaan konkreettinen funktio.
> | f:=(t,y)->1/y^2-y*t; |
> | y1; |
> | y2; |
> | y3; |
Näiden kokeilujen pohjalta johdutaan tehtäväpaperissa olevaan koodiin:
> | #read("c:\\usr\\heikki\\numsym04\\maple\\ns04.mpl"); |
> | #read("/home/apiola/opetus/numsym/04/maple/ns04.mpl"); |
> | read("/p/edu/mat-1.192/ns04.mpl"); |
> | op(taystep2); |
Testataan, kuten tehtäväpaperissa:
> | f:=(t,y)->y*t; |
> | dsolve({diff(y(t),t)=t*y(t),y(0)=1},y(t)); |
> | ytarkka:=subs(%,y(t)); |
> | Y[0]:=1: h:=0.1: T:=seq(k*h,k=0..10); for k to 11 do Y[k]:=taystep2(f,T[k],Y[k-1],h): end do: |
> | pisteet:=seq([T[k],Y[k-1]],k=1..11); |
> | plot([[pisteet],ytarkka+0.01],t=0..1); # Nostetaan tarkkaa hiukan, jotta erottuvat |
Nyt on helppoa jatkaa Taylorin 4:nnen asteen menetelmään:
> | restart: |
Warning, the name changecoords has been redefined
> | #read("c:\\usr\\heikki\\numsym04\\maple\\ns04.mpl"); |
> | #read("/home/apiola/opetus/numsym/04/maple/ns04.mpl"); read("/p/edu/mat-1.192/ns04.mpl"); |
> | op(tay4step); |
Kaksi viimeistä riviä edustavat "laiskan Maple-miehen" tapaa kerätä kunkin h:n potenssin kertoimet. (Varmaan löytyy suorempiakin.). Tässä täytyy vielä konvertoida tulos polynomiksi, sillä taylor palauttaa "sarja-tietorakenteen", jossa on jäännöstermi mukana.
(Kts. ?taylor)
> | f:=(t,y)->y*t; |
Kiintoisaa on verrata tay4step:iä ja rk4step:iä
> | op(rk4step); |
> | h:='h': |
> | rk:=rk4step(f,t0,y0,h): |
> | t4:=tay4step(f,t0,y0,'h'): |
> | rk; |
> | t4; |
Tässä tapauksessa taylorin menetelmä antaa jopa yksinkertaisemman lausekkeen.
> | tay4step(f,t0,y0,0.1); |
> | f:=(t,y)->1/y^2-y*t; |
> | t4:=tay4step(f,t0,y0,h); |
> | t4step_f:=unapply(t4,t0,y0,h); |
Meillä on nyt funktio t4step_f , joka suorittaa Taylorin 4:nnen asteen menetelmän askeleen juuri tähän konkreettiseen funktioon f liittyen. Huomataan, että unapplyn käyttö on olennaista. Tässä vaiheessa voidaan
tehdä haluttuja sievennysoperaatioita.
Kirjoitetaan nyt funktio, joka laskee n askelta käyttäen tätä askelfunktiota. Nyt siis funktio f ei enää ole argumenttilistalla, vaan sen vaikutus tulee mukaan kutsuttaessa t4step_f -funktiota.
> | op(tay4_f); |
Kutsuparametreina ovat vali , alkuarvo välin alkupisteessä ja askelten lukumäärä. (Huomaa, että vali j a y0 ovat aivan samoin kuin Matlabin ODE-funktioissa. Askeleen suorittava t4step_f on "kolvattu" kiinteän nimisenä koodiin.
> | T4pisteet:=tay4_f([1,2],1,10); |
Siistiä!
Verrataan RK4-pisteisiin:
> | # op(RK4); |
> | RK4(f,[1,2],1,10); |
> |
> | T4pisteet[-1][-1]; # Viimeisen parin viimeinen alkio, eli y(2):n approksimaatio |
> | RK4(f,[1,2],1,10)[-1][-1]; |
> | tay4_f([1,2],1,10)[-1][-1]; # Yllä oleva yhdellä rivillä. |
5:n numeron tarkkuudella samat (6:nnessa yhden numeron ero). Kumpi on tarkempi?
Nyt päästään tutkimaan tehtävän tilannetta
Katsotaan tehtävässä kysyttyja askelpituuksia.
> | Digits:=10: # Oletusarvo |
> | tayarvot:=seq(tay4_f([1,2],1,2^k)[-1][-1],k=2..10); |
> | RKarvot:=seq(RK4(f,[1,2],1,2^k)[-1][-1],k=2..10); |
Mikä on tarkin arvo? Turha toivoa, että kahta viimeistä numeroa voisi saada oikein. Askelpituuden puolittaminen ei loputtomasti lisää tarkkuutta, vaan jostain alkaen tulokset alkavat huonontua.
Otetaan viimeistä edellisten keskiarvo:
> | tarkinkai:=.8235678976; |
> | tayvirheet:=seq((tarkinkai-tayarvot[k]),k=1..nops([tayarvot])-1); |
> | RKvirheet:=seq((tarkinkai-RKarvot[k]),k=1..nops([tayarvot])-1); |
> | seq(tayvirheet[i+1]/tayvirheet[i],i=1..nops([tayarvot])-2); |
> | seq(RKvirheet[i+1]/RKvirheet[i],i=1..nops([tayarvot])-2); |
> | taulukko:=<<h,seq(2.^(-k),k=2..9)>|<Taylor4,tayarvot[1..-2]>|<RK4,RKarvot[1..-2]>|<Tay4Virhe,tayvirheet>|<RK4Virhe,RKvirheet>>; |
> | RKarvot[1..-2]; |
> |
> |
> | tayarvot20:=seq(tay4_f([1,2],1,2^k)[-1][-1],k=2..10); |
> | RKarvot20:=seq(RK4(f,[1,2],1,2^k)[-1][-1],k=2..10); |
> | melkeintarkka:=tayarvot20[-1]; |
> | Digits:=10: |
> | tayvirheet:=seq((melkeintarkka-tayarvot[k]),k=1..nops([tayarvot20])); |
> | RKvirheet:=seq((melkeintarkka-RKarvot[k]),k=1..nops([tayarvot])); |
> | seq(tayvirheet[i+1]/tayvirheet[i],i=1..nops([tayarvot20])-1); |
> | seq(RKvirheet[i+1]/RKvirheet[i],i=1..nops([tayarvot])-1); |
> |
> | taulukko:=<<h,seq(2.^(-k),k=2..10)>|<Taylor4,tayarvot>|<RK4,RKarvot>|<Tay4Virhe,tayvirheet>|<RK4Virhe,RKvirheet>>; |
Kysymys: Kuinka monta oikeaa numeroa?
Kun taulukoa katsoo, niin h=2^(-5) on pienin askel, jolla virhe pienenee normaalin näköisesti. Siitä alaspäin tulokset jopa huononevat.
Lisätään ensin laskentatarkkuutta.
> | tarkinkai; # Ainakin 2 viimeistä numeroa ovat "satunnaislukuja". |
Sama tarkemmin.
> | Digits:=20: |
> | tayarvot20:=seq(tay4_f([1,2],1,2^k)[-1][-1],k=2..10); |
> | RKarvot20:=seq(RK4(f,[1,2],1,2^k)[-1][-1],k=2..10); |
> | tarkinkai20:=tayarvot20[-1]; |
> | tayvirheet20:=seq((tarkinkai20-tayarvot[k]),k=1..nops([tayarvot20])); |
> | RKvirheet20:=seq((tarkinkai20-RKarvot[k]),k=1..nops([tayarvot20])); |
> | 2.^(-4); |
Globaalin virheen tulisi pienetyä suunnilleen tällä kertoimella, kun askel puolittuu.
> | seq(tayvirheet20[i+1]/tayvirheet20[i],i=1..nops([tayarvot20])-1); |
> | seq(RKvirheet20[i+1]/RKvirheet20[i],i=1..nops([tayarvot20])-1); |
Virhekäytös on järkevännäköistä arvoon i=5 saakka, eli h=1/2^6 = 1/64.
> | tayarvot20[5],RKarvot20[5]; |
Kertalukua lisää, yleinen Taystep
> | restart: |
Warning, the name changecoords has been redefined
> | #read("c:\\usr\\heikki\\numsym04\\maple\\ns04.mpl"); |
> | #read("/home/apiola/opetus/numsym/04/maple/ns04.mpl"); read("/p/edu/mat-1.192/ns04.mpl"); |
> |
> | op(Taystep); |
> | f:=(t,y)->1/y^2-y*t; |
> | MainTaylor(4,f,[1,2],1,10); |
> | MainTaylor(5,f,[1,2],1,10); |
> | MainTaylor(6,f,[1,2],1,10); |
> | MainTaylor(6,f,[1,2],1,20); |
> | melkeintarkka; |
> | MainTaylor(7,f,[1,2],1,10); |
> |
Kertalukua 7 olevalla Taylorilla päästään jo 10:llä askeleella tarkkaan (10 numeroa) tulokseen.
> |
Analyyttinen ratkaisu on tolkuton, "numeric" voi aiheuttaa työarkin korruptoitumista.
> | #nratk:=dsolve({diff(y(t),t)=f(t,y(t)),y(1)=1},y(t),numeric); |
> |