Luennot viikolla 6

8-10.2.2000

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

Kun koodeja kertyy, niitä kannattaa kerätä tekstitiesostoon, jonka voi lukea Mapleen. komennolla:

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

> op(lag);

Diskretointiapufunktio:

Mikään ei meitä estä määrittelemästä Matlabista tuttua linspace-funktiota: (Otamme Maplessa oletukseksi 10.)

> op(linspace);

>

Interpolointivirhe voidaan ohjelmoida funktioksi, jolle annetaan interpoloitavan funktion lauseke ja interpolointipisteet sekä

symboli, jolla muuttujaa merkitään. Otetaan globaali symboli samaan rooliin, jossa se aina on.

> op(interR);

Testataan sitä pienellä esimerkillä: sin-funktion interpolointi pisteissä 0,1,2 .

> interR(sin,[0,1,2],x);

Käyttösimerkki oli työarkilla Lv5.mws . Tässä kirjoitimme interR-funktion niin, että xd ja x sijoitettiin parametrilistaan, se on toki turvallisempaa ja oikeaoppisempaa ohjelmointia kuin aiempi globaalien suruton viljelytyyli,

Interpolaatioesimerkki (kertausta Lv5:stä)

Kerrataan nyt vielä esimerkki (ei tatrvitse katsoa, jos kyllästyttää)

Max-virhe koko välillä

Yläraja max-virheelle koko välillä

Splinit

> readlib(spline):

> f:=x->1/(1+25*x^2):

> n:=4:xd:=linspace(-1,1,n);

> yd:=map(f,xd);

> s:=spline(xd,yd,x);

> plot([f(x),s],x=-1..1);

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

> n:=10:xd:=linspace(-1,1,n);
yd:=map(f,xd);
s:=spline(xd,yd,x):
plot([f(x),s],x=-1..1);
plot(f(x)-s,x=-1..1);

Numeerinen integrointi

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

Määrittelimme funktiot trapetsi ja simpsoni. Lisäksi student :ssa on trapezoid ja simpson , mutta ilman virhetermiä. Muuten ne ovat yhdenmukaiset int -funktion kanssa, mikä on kaunista.

> op(trapetsi);op(simpsoni);print(simpsonvirhe);

> with(student);

Katsotaan vielä esimerkkinä integraalin laskemista trapetsilla + virhearviota.

> f:=x->exp(x)*cos(x);plot(f(x),x=0..Pi);

> SR:=trapetsi(f,0,Pi,10);

> S:=SR[1];R:=SR[2];

> plot(abs(R),xi=0..Pi,thickness=2);

> maxvirhe:=0.4;

> maxsuhtvirhe:=maxvirhe/abs(S);

>

Trapetsin ja Simpsonin vertailua

Lasketaan yksi esimerkki, joka näyttää laskentatöiden suuruusluokkavertailun.

> Isi:=Int(sin(x),x=0..Pi);tol:=2*10^(-5);

> n:='n':tol:='tol':trapetsivirhe(sin,0,Pi,n);Mtr:=eval(subs(xi=Pi/2,%));
simpsonvirhe(sin,0,Pi,n):Ms:=abs(eval(subs(xi=Pi/2,%)));

> ntr:=solve(Mtr=tol,n)[1];ns:=solve(Ms=tol,n)[1];

> evalf(subs(tol=2*10^(-5),ntr));evalf(subs(tol=2*10^(-5),ns));

Siis t rapetsilla tarvitaan 360 osaväliä ja Simpsonilla riittää 18 ! (Muista, oltava parillinen)

Lasketaan kolme integraalia molemmilla menetelmillä ja katsotaan, miten virhe käyttäytyy askelpituuden funktiona. Käytämme virheenä oikeaa virhettä, emme tässä selvittele virheen arvioimista jäännöstermillä, vaan oikean virheen käyttäytymistä.

> I1:=Int(exp(-x^2),x=0..1)=int(exp(-x^2),x=0..1);
I2:=4*Int(1/(1+x^2),x=0..1)=4*int(1/(1+x^2),x=0..1);
I3:=Int(1/(2+cos(x)),x=0..2*Pi)=int(1/(2+cos(x)),x=0..2*Pi);

> f1:=x->exp(-x^2);f2:=x->4/(1+x^2);f3:=x->1/(2+cos(x));

> N:=8:

> trapS1:=seq(trapetsi(f1,0,1,2^k)[1],k=1..N);
tarkka1:=evalf(lhs(I1));trapV1:=seq(tarkka1-trapS1[i],i=1..N);

> tsuhteet1:=seq(trapV1[i]/trapV1[i+1],i=1..N-1);

Tämä on hyvin sopusoinnussa trapetsivirheen O(h^2); - käyttäytymisen kanssa.

> trapS2:=seq(trapetsi(f2,0,1,2^k)[1],k=1..N);
tarkka2:=evalf(lhs(I2));trapV2:=seq(tarkka2-trapS2[i],i=1..N);tsuhteet2:=seq(trapV2[i]/trapV2[i+1],i=1..N-1);

> Digits:=20:trapS3:=seq(trapetsi(f3,0,2*Pi,2^k)[1],k=1..N);
tarkka3:=evalf(lhs(I3));trapV3:=seq(tarkka3-trapS3[i],i=1..N);tsuhteet3:=seq(trapV3[i]/trapV3[i+1],i=1..N-1);

> with(student):trapezoid(f3(x),x=0..2*Pi);evalf(%);

Warning, new definition for D

> evalf(trapezoid(f3(x),x=0..2*Pi,4));trapetsi(f3,0,2*Pi,4);

>

Tämä on hyvin sopusoinnussa trapetsivirheen - käyttäytymisen kanssa.

Entäpä Simpson:

> Digits:=10:N:=8:simpS1:=seq(simpsoni(f1,0,1,2^k)[1],k=1..N);
tarkka1:=evalf(lhs(I1));simpV1:=seq(tarkka1-simpS1[i],i=1..N);
ssuhteet1:=seq(simpV1[i]/simpV1[i+1],i=1..N-1);

Error, division by zero

> Digits:=20:N:=8:simpS1:=seq(simpsoni(f1,0,1,2^k)[1],k=1..N);
tarkka1:=evalf(lhs(I1));simpV1:=seq(tarkka1-simpS1[i],i=1..N);
ssuhteet1:=seq(simpV1[i]/simpV1[i+1],i=1..N-1);

> Digits:=10:N:=8:simpS2:=seq(simpsoni(f2,0,1,2^k)[1],k=1..N);
tarkka2:=evalf(lhs(I2));simpV2:=seq(tarkka2-simpS2[i],i=1..N);
ssuhteet2:=seq(simpV2[i]/simpV2[i+1],i=1..N-1);

Error, division by zero

> Digits:=20:N:=8:simpS2:=seq(simpsoni(f2,0,1,2^k)[1],k=1..N);
tarkka2:=evalf(lhs(I2));simpV2:=seq(tarkka2-simpS2[i],i=1..N);
ssuhteet2:=seq(simpV2[i]/simpV2[i+1],i=1..N-1);

> Digits:=10:N:=8:simpS3:=seq(simpsoni(f3,0,2*Pi,2^k)[1],k=1..N);
tarkka3:=evalf(lhs(I3));simpV3:=seq(tarkka3-simpS3[i],i=1..N);
ssuhteet3:=seq(simpV3[i]/simpV3[i+1],i=1..N-1);

> Digits:=20:N:=8:simpS3:=seq(simpsoni(f3,0,2*Pi,2^k)[1],k=1..N);
tarkka3:=evalf(lhs(I3));simpV3:=seq(tarkka3-simpS3[i],i=1..N);
ssuhteet3:=seq(simpV3[i]/simpV3[i+1],i=1..N-1);

Error, division by zero

> Digits:=50:N:=8:simpS3:=seq(simpsoni(f3,0,2*Pi,2^k)[1],k=1..N);
tarkka3:=evalf(lhs(I3));simpV3:=seq(tarkka3-simpS3[i],i=1..N);
ssuhteet3:=seq(simpV3[i]/simpV3[i+1],i=1..N-1);

>

> Digits:=10: # Palataan normaaliin.

Luentotehtävä: (8.2.)

Gaussin kvadratuuri

> f0:=x->1;f1:=x->x;f2:=x->x^2;f3:=x->x^3;

> yht0:=int(f0(x),x=-1..1)=c1*f0(x0)+c2*

> f0(x1);

> yht1:=int(f1(x),x=-1..1)=c1*f1(x0)+c2*f1(x1);
yht2:=int(f2(x),x=-1..1)=c1*f2(x0)+c2*f2(x1);
yht3:=int(f3(x),x=-1..1)=c1*f3(x0)+c2*f3(x1);

>

> solve({yht0,yht1,yht2,yht3},{x0,x1,c1,c2});

> allvalues(%);

Loistavaa!

Korkeamman kertaluvun menetelmissä tarvitaan ortogonaalipolynomeja

with(orthopoly); Legendren polynomit

> with(orthopoly);

Warning, new definition for L

> P(0,x);P(1,x);P(2,x);

> solve(P(2,x)=0,x);

Hämmästyttävää! Jos olisi niin ihmeellisesti, että solmut saataisiin Legengden polynomien nollakohtina, niin painokertoimet olisi helppo ratkaista. Niiden suhteenhan yhtälösysteemi on lineaarinen. Toisaalta ne epäilemättä saadaan myös solmuihin liittyvien Lagrangen kertojapolynomien integraaleina. Kumpi on parempi tapa laskea?

Tehtävä viikolle 7

Moniulotteinen integrointi

Yhteenveto menetelmistä

Tietokonearitmetiikkaa, laskennan tehokkuudesta

Polynomien laskennasta

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

> q:=convert(p,horner);

> with(codegen[cost]):

Warning, new definition for fortran

Warning, new definition for makeproc

> cost(p);cost(q);

>

Aritmetiikkaa

> x:=0.57812;q:=x:

> b[1]:=floor(2*q);q:=frac(2*q);

> b[2]:=floor(2*q);q:=frac(2*q);

> b[3]:=floor(2*q);q:=frac(2*q);b[4]:=floor(2*q);q:=frac(2*q);

> b[5]:=floor(2*q);q:=frac(2*q);b[6]:=floor(2*q);q:=frac(2*q);

> b[7]:=floor(2*q);q:=frac(2*q);b[8]:=floor(2*q);q:=frac(2*q);

> b[9]:=floor(2*q);q:=frac(2*q);b[10]:=floor(2*q);q:=frac(2*q);

> b[11]:=floor(2*q);q:=frac(2*q);b[12]:=floor(2*q);q:=frac(2*q);

> b:='b':

> b:='b':x:=0.578125:q:=x:

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

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

> polyval0(cc,0.5);

> polyval([u,v,w],x);

> append:=(L1,L2)->[op(L1),op(L2)];

> Digits:=7:

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

>

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

> c;

> Digits:=20:

> Digits:=5:

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

> .49e-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.);

Mitenkähän siinä huonossa ensimmäisessä tavassa lopulta käy?

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

No on ne nyt sinnepäin, mutta on siinä heittoakin melkoisesti. Ei tunnu uskottavalta, että enää voisi korjaantua.

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

>

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

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

>