Luennot viikolla 7

12-14.2.2000

Kts. myös Lvhtml7.html

Ohjelmatiedosto /p/edu/mat-1.414/maple/ohjelmat.mpl

> read("../v2/ohjelmat.mpl");

Uusia funktioita ohjelmat.mpl -tiedostossa:

> op(polyval);op(fliplr);op(append);

Tietokonearitmetiikkaa, laskennan tehokkuudesta

Polynomien laskennasta

> p:=-8*x^5+7*x^4-6*x^3-5*x^2-4*x+3;nops(p);

> q:=convert(p,horner);

> with(codegen[cost]):

Warning, new definition for fortran

> cost(p);cost(q);

> coeffs(p,x); # Hyödyllinen funktio, mutta jostain syystä järjestys on mitä sattuu => hyödytön.

Määritellään oma:

> op(kertoimet);

> kertoimet(p,x);

> polyval(kertoimet(p,x),x);expand(%);

Aritmetiikkaa

Esimerkki dec2bin-muunnoksesta, kun 0 < x < 1.

> b:='b':

> x:=0.578125:q:=x:

> for i from 1 to 10 do b[i]:=floor(2*q);q:=frac(2*q);od: c:=[seq(b[i],i=1..10)];

> puoli:=1/2:1*puoli*1+0*puoli^2+0*puoli^3+1*puoli^4+0*puoli^5+1*puoli^6;evalf(%);

> cc:=append([0],c);

> polyval0(cc,0.5);

Hyvinhän tuo toimii. Lasketaan vielä luentoesimerkki:

> b:='b':x:=0.1:q:=x:N:=20:
for i from 1 to N do b[i]:=floor(2*q);q:=frac(2*q);od: c:=[seq(b[i],i=1..N)];
polyval0(append([0],c),0.5);

>

> b:='b':x:=0.1:q:=x:N:=40:
for i from 1 to N do b[i]:=floor(2*q);q:=frac(2*q);od: c:=[seq(b[i],i=1..N)];
polyval0(append([0],c),0.5);

Siis 40-bittinen esitys antaa 10 desimaalinumeron tarkkuudella oikein.

Ylivuoto, alivuoto, skaalaus

Esim: , kun a ja b eri kertaluokkaa.

> Digits:=7:

> a:=10.^15;b:=1;sqrt('a'^2+'b'^2)=sqrt(a^2+b^2);

>

> Digits:=32:a:=10.^15;b:=1;c:=sqrt(a^2+b^2);

>

> Digits:=5:

> seq([j,2.^(-j),1.+2^(-j)],j=10..20);

Nähdään, että j:n arvolla 15 jäädään1. kertaa samaan 1:een, joten kone-epsilon on n. . Jos halutaan katsoa ihan tarkkaan, niin: 0.49e-4

> .49e-4+1.0;.5e-4+1.0;

>

Merkitsevien numeroiden kumoutuminen

> f:=x->x*(sqrt(x+1)-sqrt(x));

> Digits:=10:

> matrix([seq([10^k,evalf(f(10^k),6),evalf(f(10^k),20)],k=1..6)]);

> Digits:=6:

> x:=100.:sqrt(x+1),sqrt(x);

> sqrt(x+1)-sqrt(x);

> Digits:=9:

> x:=100.:sqrt(x+1),sqrt(x);

>

> sqrt(x+1)-sqrt(x);

>

6:lla numerolla lasketussa tuloksessa oli vain 3 oikeaa numeroa, 3 alkunumeroa kumoutuivat vähennyslaskussa ja siten lasku tehtiin vain 3:n numeron tarkkuudella.

Tai ehkä paremmin ilmaistuna: Laskut pitää suorittaa 3:n lisänumeron tarkkuudella, jotta tulos on 6:lla numerolla oikein.

Kirjoitetaan kaava muotoon, jossa lavennetaan "liittoluvulla":

> g:=x->x/(sqrt(x+1)+sqrt(x));

> x:='x':simplify(f(x)-g(x));

> matrix([seq([10^k,evalf(g(10^k),6),evalf(g(10^k),20)],k=1..6)]);

Nytpä tuli hienosti !

Katsotaanpa Taylorin polynomin laskentaa:

> tayexp:=seq(convert(taylor(exp(x),x=0,j),polynom),j=0..10);

> evalf(subs(x=-6.0,[tayexp]),5);exp(-6.);

> ykspertayexp:=seq(1/(convert(taylor(exp(x),x=0,j),polynom)),j=1..10);

> evalf(subs(x=6.0,[ykspertayexp]),5);exp(-6.);

Katsotaan alkuperäistä laskuatapaa pitemmälle. Laskemme niin pitkälle, että osasummat eivät enää muutu.

> tayexp:=seq(convert(taylor(exp(x),x=0,j),polynom),j=25..30):evalf(subs(x=-6.0,[tayexp]),5);exp(-6.);

> (.247875218e-2-.68984e-3)/.247875218e-2;

Suhteellisn virheen lopullinen arvo on yli 70%.

Miksi tarkuus katoaa?

Alkupäässä lasketaan varsin suuria erimerkkisiä arvoja, lopputulos on hyvin pieni. Kun lasketaan kiinteällä suhteellisella tarkkuudella, niin absoluuttinen virhe tulee suurten termien kohdalla myös suureksi. Kun lopputulos on pieni, niin sen suhteellinen virhe tulee suureksi. (Pahimmassa tapauksessa isojen termien absoluuttinen virhe voi olla samaa suuruusluokkaa kuin koko lopppusumma, no kaikkihan tietää, ettei siitä hyvää seuraa.)

Yleensä, jos lasketaan summaa, jossa on pieniä ja suuria positiivisia, niin paras tulos saadaan laskemalla pienimmästä suurimpaan. (Päinvastoin laskettaessa pienen termin lisääminen suureen summaan voi hukuttaa pienen kokonaan (vrt. kone-eps).)

>

Kts. myös Lvhtml7.html

Epälineaariset yhtälöt, Newtonin menetelmä, iteraatiot

Kts. myös NCnewton.mws

[RA] s. 97

[HAM] s. 67-69

[Israel] s. 55 ->

> f:='f':N:=x->evalf(x-f(x)/D(f)(x));

> f:=x->x*exp(-x)-abs(sin(10*x));

> plot(f(x),x=-1..3);

> x1:=(N@@5)(0.25);

> x1:=(N@@5)(x1);

> plot(f(x),x=0..0.5);

Newtonin menetelmä merkitsee yllä olevan funktion N iterointia. Iterointi Maplessa saadaan aikaan @@:lla

>

Tarkastellaan vähän siistimmin käyttäytyvää funktiota

> f:=x->1+sin(2*x)-x;

Iteraatiojono voitaisiin muodostaa elegantisi

> seq((N@@k)(1.2),k=0..6);

Newton suppenee yleensä nopeasti, mutta on tämä hieman tuhlailevan näköistä, kun aina iteroidaan alusta. Vähemmän elegantti, mutta vähemmän tuhlaava on yksinkertainen for-rakenne:

> x[0]:=2.0;

> for i from 0 to 30 do x[i+1]:=N(x[i]) od:

> with(plots):

> display([plot([[x[0],0],[x[0],f(x[0])],[N(x[0]),0]],color=blue),plot(f,0..2.5,color=red)]);

Graafista havainnollistusta varten tehdään Robert Israelin elegantin ideologian mukaiset grafiikka-arvoiset funktiot.

> Delta:=x0->plot([[x0,0],[x0,f(x0)],[N(x0),0]],color=blue);

> fkuva:=display([plot(f,0..3)]):

> display([seq(display([fkuva,Delta(x[k])]),k=0..30)],insequence=true);

> display([seq(display([fkuva,Delta(x[k])]),k=0..30)]);

>

Normaalitapauksen luonteeseen kuuluu, että animaatio häipyy nopeasti näkyvistä.