Kotiin

Luentoja viikolla 7 (14.2 - 16.2. 2000)

Ti 15.2. esitettiin

Ke 16.2.

Sisällys

Kirjallisuutta

Tietokonearitmetiikkaa

Katsotaan vielä polynomien evaluointiasiaa. Vrt. myös Lv7.mws Kirjoitimme Matlab-funktion horner.m (Otettin suoraan [BB]-kirjasta.) Tämä on matlab-hakemistossamme. Ellei, niin cut/paste ja UNIX:ssa: 
    maalaus hiiren vasemmalla
    cat > horner.m
    hiiren keskimmäinen kopioi tähän
    CTR-D   % Tiedoston loppu

function p=horner(c,x)
l=length(c);p=c(1);
for k=2:l
   p=p.*x+c(k);
end;

Istunto, jossa Verrataan flopseja oman hornerin ja Matlabin polyvalin kesken 


%Ensin skalaarilaskenta:
%%%%%%%%%%%%%%%%%%%%%%%%
x=5;
count1=zeros(1,10); count2=zeros(1,10);
value1=zeros(1,10);value2=zeros(1,10);
for n=1:10
   c=rand(1,n+1);
   flops(0);value1(n)=polyval(c,x);count1(n)=flops;
   flops(0);value2(n)=horner(c,x);count2(n)=flops;
end;
n=1:10;
[n' count1' count2']

plot(n,count1,n,count2)
title('Skalaaritapaus')
xlabel('Polynomin asteluku'),ylabel('Liukulukuoperaatioiden lukumäärä')
text(n(5),count1(5),'Matlabin polyval')
text(n(5),count2(5),'Oma Horner')
value1-value2

% Sitten vektori:
%%%%%%%%%%%%%%%%%%%%

x=linspace(0,1,200);
N=30
count1=zeros(1,N); count2=zeros(1,N);
%value1=zeros(1,10);value2=zeros(1,10);
for n=1:N
   c=rand(1,n+1);
   flops(0);polyval(c,x);count1(n)=flops;
   flops(0);horner(c,x);count2(n)=flops;
end;
n=1:N;
[n' count1' count2']


plot(n,count1,n,count2)
title('Vektorilaskenta')
xlabel('Polynomin asteluku'),ylabel('Liukulukuoperaatioiden lukumäärä')
text(n(floor(N/3)),count1(floor(N/2)),'Matlabin polyval')
text(n(floor(N/3)),count2(floor(N/3)),'Oma Horner')
grid

type polyval

IEEE-standardi

32-bittinen

merkkieksp (8 bittiä) mantissa (23 bittiä)

Esim:

01 0 0 1 1 0 1 11 0 1 0 1 1 1 1 1 1 1 1 0 0 1 1 0 1 0 1 1 1 1

Normalisointi: m=1.b1b2 ... b23 , 1 <= m < 2

-127 < e <= 127

Lukualue 


           [L,U] = [2-126, 2127] = [10-38, 1038]
Pienin: (1.000 ... 0)2 2-126 = 1.17549 10-38
Suurin (1.111 ... 1)2 2127 = 3.4 1038

64-bittinen (Matlab)

merkkieksp (11 bittiä) mantissa (52 bittiä)

Lukualue: [10-307, 10307], tarkkuus n. 15 numeroa.

IEEE-standardiin kuuluvat symbolit NaN ("Not A Number") ja Inf ("Infinity"). 

1/0 -> Inf
0/0, 0*Inf, arcsin 2, ... antavat NaN.
Viimeksi mainittu on esimerkki funktion soveltamisesta lukuun, joka ei kuulu funktion määrittelyjoukkoon.

Pyöristys

Yli/alivuoto johtuu eksponenttialueen rajoista (OFL, UFL), kun taas pyöristysvirhe aiheutuu mantissan rajallisesta pituudesta.

Pyöristystä on kahta lajia, katkaisu ("truncate") ja oikea pyöristys ("round"). Jälkimmäinen on huomattavasti parempi, ei vain siksi, että se on hieman tarkempi, vaan ennen kaikkea siksi, että siinä tapahtuu tilstollisluonteista pyöristysvirheiden kumoutumista, jolloin laskentatulokset ovat tavallisesti ylärajaennusteita selvästi parempia.

Pyöristys tehdään koulussa opitun säännön mukaan. Yksi esimerkki riittäköön: Lasketaan 5-numeroisessa 10-järjestelmän aritmetiikassa: Olkoon x=2.64575... ja esitetään se 5:llä numerolla normalisoidussa muodossa ja oikein pyöristettynä. , joten float(x)=2.6458 . Normalisoinnin takia |x| > 1 ja siten 

         |x-float(x)|/|x| < 5*10-5  = 1/2*10-4
Yleistyy välittömästi, jos 5:n sijasta on yleinen mantissan pituus t ja 10:n sijalla yleinen kantaluku beta: 
         |x-float(x)|/|x| < 1/2*beta-t+1
Erityisesti binäärijärjestelmässä: 
         |x-float(x)|/|x| < 2-t

Kone-epsilon

Merk eps=1/2 beta1-t

Katsotaan edelleen 10-järj. 5-pituista aritmetiikkaa. (beta=10, t=5). 

                1.0000|
                0.0000|49999
               --------------
                1.0000|4
Yhteenlaskun tulos = 1.0 Kone-epsilon voidaan myös määritellä suurimpana liukulukuna eps , joka toteuttaa ehdon 1+eps=1 (tai pienimpänä positiivisena, jolle 1+eps > 1 .

Suuruusluokka on tärkeä, 2*eps on ihan yhtä hyvä kuin eps.

Ensin näin (editoidaan outputista turhia pois): 

>> format long; n=50;1+2.^-(1:n)

ans =

  Columns 1 through 4 

   1.50000000000000   1.25000000000000   1.12500000000000   1.06250000000000

  ...............

  Columns 45 through 48 

   1.00000000000003   1.00000000000001   1.00000000000001   1.00000000000000

  Columns 49 through 50 

   1.00000000000000   1.00000000000000
Nyt olemme viisaampia: 
format long;N=40:50; [N' (1+2.^-N)']

>> format long;N=40:50; [N' (1+2.^-N)']

ans =

  40.00000000000000   1.00000000000091
  41.00000000000000   1.00000000000045
  42.00000000000000   1.00000000000023
  43.00000000000000   1.00000000000011
  44.00000000000000   1.00000000000006
  45.00000000000000   1.00000000000003
  46.00000000000000   1.00000000000001
  47.00000000000000   1.00000000000001
  48.00000000000000   1.00000000000000
  49.00000000000000   1.00000000000000
  50.00000000000000   1.00000000000000

>> 2^-48
ans =

     3.552713678800501e-15

>> eps

ans =

     2.220446049250313e-16
Matlabin sisäänrakennettu kone-epsilon näyttää olevan hieman liian pieni, mutta lähes 10 ei ole näissä hommissa iso kerroin.

Merkitsevien numeroiden katoaminen vähennyslaskussa

Jos vähennämme luvut 
    x=3.7214
    y=3.7202
--------------
  x-y= 0.00013 = 1.3000*10-4
Kolme viimeistä 0:aa eivät sisällä mitään informaatiota. Suhteellista tarkkuutta katoaa erotuksesta likimain senverran, kuin alkupäässä on samoja numeroita. Asia on helppo ymmärtää myös siten, että vähennyslaskussa absoluuttinen virhe on samaa luokkaa (kork. 2 kertaa) kuin termien abs. virhe, mutta kun luvut ovat lähes samat, on erotuksen suhteellinen virhe suuri.

Esim: Lasketaan exp(x) Taylorin kehitelmän avulla.

exp(x)=1+x+x2/2!+ x3/3! + ...

Lasketaan e-5.5 

    |\^/|     Maple V Release 5.1 (Helsinki University of Technology)
._|\|   |/|_. Copyright (c) 1981-1998 by Waterloo Maple Inc. All rights
 \  MAPLE  /  reserved. Maple and Maple V are registered trademarks of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.

Digits:=5:
x:=-6.5:
seq(x^k/k!,k=0..30);
seq(add(x^k/k!,k=0..n),n=0..30);
exp(x);

> Digits:=5:
> x:=-6.5:
> seq(x^k/k!,k=0..30);
1, -6.5, 21.125, -45.772, 74.379, -96.692, 104.75, -97.266, 79.030, -57.077,

    37.100, -21.923, 11.875, -5.9373, 2.7566, -1.1946, .48526, -.18555,

    .067003, -.022922, .0074500, -.0023059, .00068128, -.00019254, .000052145,

                          -5            -6           -6            -7
    -.000013558, .33894 10  , -.81598 10  , .18942 10  , -.42457 10  ,

             -8
    .91992 10

> seq(add(x^k/k!,k=0..n),n=0..30);
1, -5.5, 15.625, -30.147, 44.232, -52.460, 52.290, -44.976, 34.054, -23.023,

    14.077, -7.846, 4.029, -1.9083, .8483, -.3463, .13896, -.04659, .020413,

    -.002509, .0049410, .0026351, .0033164, .0031239, .0031760, .0031624,

    .0031658, .0031650, .0031652, .0031652, .0031652

> exp(x);
                                   .0015034

>  (.0015034-.0031652)/.0015034;
                                    -1.1054
Tämän pitemmälle ei kannata laskea, koska summa ei enää muutu. Suhteellinen virhe on yli 110% .

Entä, jos laskemme kaavalla exp(-6.5)=1/exp(6.5) . 

Digits:=5:
x:=6.5:
seq(1/add(x^k/k!,k=0..n),n=0..30);
exp(-x);
> Digits:=5:
> x:=6.5:
> seq(1/add(x^k/k!,k=0..n),n=0..30);
1, .13333, .034934, .013441, .0067213, .0040738, .0028553, .0022347, .0018993,

    .0017135, .0016111, .0015561, .0015279, .0015141, .0015078, .0015051,

    .0015040, .0015036, .0015034, .0015034, .0015034, .0015034, .0015034,

    .0015034, .0015034, .0015034, .0015034, .0015034, .0015034, .0015034,

    .0015034


> exp(-x);
                                   .0015034

Kiintopisteiteraatio

Kts. myös Maple: Lv7iter.mws

Alla oleva Maple-teksti on saatu FILE-valikon EXPORT -> txt.

Cobweb Maplella 


#     Määritellään kiintopisteiteraatiota havainnollistavia
# apufunktioita.
#  Määrittelemme ensin funktion porras, joka palauttaa
# grafiikkaobjektin, joka puolestaan edustaa "porrasta"
# [x,x],[x,f(x)],[f(x),f(x)]. Annamme f:n olla yleinen, toistaiseksi
# määrittelemätön globaali symboli.
# 
# Tämän tyylin idea on peräisin Robert Israelilta.
#                                                        
> porras:=x->plot([[x,x],[x,f(x)],[f(x),f(x)]]):
# Käyttöesimerkki, leikkaa/liimaa sovellukseesi.
# Ensin funktiomäärittely ja alkupiste.
> with(plots):
> f:=x->3.0*x-x^2: x0:=1.4:
> 
# Seuraavaksi muodostamme iteraatiojonon ja piirrämme koko prosessin.
# Erityisesti seuraava solu kannattaa leikata/liimata
> X:=[seq((f@@k)(x0),k=0..30)]; # Tämän voit synnyttää myös
> for-silmukalla
> fjaxkuva:=plot([x,f(x)],x=1.3 .. 2.6,color=[blue,black]):
> display({fjaxkuva,seq(porras(X[i]),i=1..nops(X))});
> 
> display([seq(display([fjaxkuva,porras(X[i])]),i=1..nops(X))],insequenc
> e=true);

Cobweb Matlabilla

Luennolla 1.3. esitettyyn tyyliin ajateltuna

Newtonin menetelmä ja fraktaalit

Suoritetaan Newtonin menetelmä vaikkapa polynomille lähtien kompleksitason pistehilasta. Matlabissa standarditapa xy-tason pistehilan muodostamiseen on meshgrid . Sitä käytellään ahkerasti kurssin kakkososassa, aina kun piirretään pintakuvia Matlabilla. 

clear
x=linspace(-1,1,6);y=x;
[X,Y]=meshgrid(x,y)
cgrid=X+i*Y

% Varustetaan jokainen ruutu omalla värillään, ensin numeroidaan 
% kokonaisluvuilla 1:36:

Z=reshape(1:36,6,6);Z=Z'

clf
plot(X(:),Y(:),'.r');hold on % Piirretään väri kerrallaan
Oikeassa tilanteessa valitaan alkup. hilasta se alkupistejoukko, joka näyttää suppenevan kohti "punaista 0-kohtaa". Jne.

Toinen tapa on käyttää pintapiiroksen colormappiä: 

clf
surf(X,Y,Z);view(2)
%shading flat
%colormap(flipud(jet))
colorbar

Otetaan enemmän pisteitä (ja muistetaan nyt puolipisteet) 

clear
x=linspace(-1,1,50);y=x;
[X,Y]=meshgrid(x,y);
cgrid=X+i*Y;

Z=reshape(rem(1:50*50,8),50,50);Z=Z';
clf
surf(X,Y,Z);view(2)
shading flat
colormap(flipud(jet))
%colorbar
Huomaa, että kun valitaan jotain ehtoja toteuttavia alkioita, niin usein on kätevää luiskauttaa pitkiksi vektoreiksi tyylin X(:) . Lopuksi voidaan koota takaisin matriisiksi reshape:llä

Kirjassa Eva PE Anders Sjöberg s. 296-298 on ohjelma mandelbrotprog.m

Binäärihaku, eli bisektio (deka- etc sektio)

Alkuperäinen väli olkoon [a,b] ja f vaihtakoon merkkiä.

  1. alpha:=a: beta:=b: Laske fal:=f(alpha): fbe:=f(beta):
  2. gamma:=1/2(alpha+beta): Laske fgam:=f(gamma):
  3. Jos sign(fgam)=sign(fal) niin alpha:=gamma: muuten beta:=gamma;
  4. Jos beta - alpha > tol, palaa kohtaan 2.
Toteuta Matlabilla (ja/tai Maplella). Selvitä, monellako askeleella saavutetaan haluttu tarkkuus.

Skalaariorientoitunut Matlab-toteutus on tässä

Miksei jaeta vaikka sataan osaan?

Matlabilla tunnetusti vektorioperaatiot viuhuvat lujaa. Katsotaan malliksi vaikka sin-funktion 0-kohdan hakua välillä [-1,1]. 
a=-1; b=1
x=linspace(a,b);y=sin(x);I=find(sign(y(1)) ~= sign(y)); I1=I(1);a=x(I1-1),b=x(I1)
Voit keksiä elegantimmankin, mutta periaate on tänänkaltainen. Kenties olisi mahdollista paikantaa kaikki merkinvaihtokohdat samanaikaisesti. No mutta hyvä näinkin.

Tästä funktio! Voit tietysi tehdä n_osaa, jossa n on parametrina. (Ei ole synonyymi ilmaisulle "en osaa" (edellisestä puuttuu aspiraatio).) Formuloi ja todista "sataosaalause", "n_osaa-lause" (ilman aspiraatiota). (Ai pahus, nykyisin kai harrastetaan ännää eikä ennää, hyvä sanaleikki meni osittain pilalle.)

Mikä on virhearvio? Ehkä olis kiintoisaa verrata aikoja bisektion ja sataosaa:n välillä. (Matlab:ssa tic, toc).

Nämä menetelmät rajaavat (sulkevat) 0-kohdan ("bracket a zero"). Alkuun

Hybridimenetelmät

Alkuun

Luku4

Alkuun

KotiinKotiin

sivuja ylläpitää Heikki Apiola -- etusivu http://www.math.hut.fi/v2/ -- uutisryhmä opinnot.mat.v2