Első labor: A lebegőpontos számítások alapjai

Contents

Túlcsordulás

A Matlab dupla pontosságú lebegőpontos ábrázolásában 52 bit jut a mantisszára (értékes jegyek száma), 1 előjelbit, 11 jegy a karakterisztikára (de ebből egy előjelbit). Ennek következménye:

(2^53+1)-2^53
ans =

     0

$\epsilon_0, \epsilon_1$

A legkisebb ábrázolható szám: $\epsilon_0$

realmin
ans =

    2.225073858507201e-308

Viszont ez nem azt jelenti, hogy ilyen finomsággal tudjuk ábrázolni a lebegőpontos számokat. Ezt az 1 utáni legkisebb ábrázolható szám adja meg: $\epsilon_1$

format long
1+2^-52
ans =

   1.000000000000000

1+2^-53
ans =

     1

53/10*3
ans =

  15.899999999999999

Ez nagyjából 15 értékes jegynek felel meg, ennyi a maximális szóbajövő számítási pontosság. Sajnos ez sem mindig elérhető.

Összeadás

s=1;
p=0;
for i=1:10
    s=s+10^(-16);
    p=p+10^(-16);
end,
s
s =

     1

p=p+1
p =

   1.000000000000001

Azaz érdemes először a kicsi tagokat összeadni, nem mindegy a sorrend.

Kivonás

Vegyük a két elvben egyforma mennyiséget: sqrt(a)-sqrt(a-1), illetve 1/(sqrt(a)+sqrt(a-1)).

c=sqrt(10^6)-sqrt(999999)
c =

     5.000001250436981e-04

d=1/(sqrt(10^6)+sqrt(999999))
d =

     5.000001250000625e-04

Kiegyszerűsödés: két közeli szám kivonasákor fellépő jegyveszteség Ha a két számban az első k jegy azonos, akkor a különbségükben az első k jegy értektelen Az elsőben az utolsó 6 jegy nem stimmel, a második a jó. Tanulság: numerikusan megbízhatatlan egymáshoz közeli számok kivonása.

Vegyük az x^2+b*x+1 egyenlet gyökeit. Ekkor a diszkrimináns: sqrt(b^2-4*a*c), azaz

b=10^9; %esetén
D=sqrt(b^2-4) % a diszkrimináns
D =

     1.000000000000000e+09

Ekkor az egyik gyök: (-b+D)/2, itt:

x1=(-b+D)/2 % ez rossz, viszont helyette az előbbi módszerrel
x1 =

     0

x11=-2/(b+D) % ami sokkal közelebb van a pontos gyökhöz (bár nem annyi)
x11 =

    -1.000000000000000e-09

Szorzás

Sajnos itt sem számíthatunk pontos eredményre minden szituációban:

(1-10^-9)*(1+10^-9)
ans =

     1

A nagy számmal való szorzás is veszélyes művelet, ugyanis felnagyíthatja a hibákat.

Osztás

A 0 közeli számokkal veszélyes osztani, például:

1/(sqrt(10^10)-sqrt(10^10-1))
ans =

     2.000002233314028e+05

helyette:

sqrt(10^10)+sqrt(10^10-1)
ans =

     1.999999999950000e+05

Saját függvények

Mi magunk is készíthetünk saját függvényeket, a már beépített függvények felhasználásával. Ennek módja function kimenet(ek)=függvénynév(bemenet(ek)) ez egy külön .m fájlba kell menteni, hogy mûködjön, és a fájl nevén kell meghívni.

A függvényüket a parancssorból tudjuk használni, ha abban a könyvtárban vagyunk, ahová mentettük. Ezt a pwd paraccsal (vagy a felsõ sorban) ellenõrizhetjük.

Vigyázat, minden a függvény belsejében szereplõ változó lokális, azaz a függvény futása után nem tudjuk elérni, például v-t sem látjuk!

Több bemenet és kimenet is lehetséges, ez a függvény például visszaadja a bementeként kapott számok összegét és szorzatát

function [osszeg, szorzat]=szamol(a,b)

osszeg=a+b;
szorzat=a*b;

Órai feladat

Készítsünk egy olyan függvényt, amely minél pontosabban oldja meg a másodfokú egyenletet (vizsgálja, hogy lehet-e kiegyszerűsödés). Az első sor legyen function [gyok1, gyok2]=masodfoku(a,b,c).

function [gyok1, gyok2]=masodfoku(a,b,c)
% Thf nem hiányos, azaz a~=0, b~=0, c~=0
D=sqrt(b^2-4*a*c);

if b>0
    gyok1=-2*c/(b+D);
    gyok2=(-b-D)/(2*a);
elseif b<0
    gyok1=(-b+D)/(2*a);
    gyok2=(-2*c)/(b-D);
end