ratk4.html

Numsym Harj. 4

18.3.04 ratkaisuja

> 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ä.

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

>

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

>