> restart:
Błędy obliczeń numerycznych
Zagadki
Zagadka 1
> a1:=evalf(Pi)*10^5;
a1 := 314159.2654
> a2:=evalf(Pi)*10^(-5);
a2 := 0.00003141592654
> a3:=a1+a2:
> is(a1=a3);
true
Zagadka 2
> x:=987654321;
x := 987654321
> y:=evalf(Pi);
y := 3.141592654
> z:=-x;
z := -987654321
> x+y+z;
3.1
Zagadka 3
> restart:
> y:=1/1*(10^x*(1+x*10^(-x))-10^x);
( -ðx)
y := 10x ( 1 ð+ð ðx 10 ) ð-ð ð10x
> eval(y,x=20.);
0.
> simplify(y);
x
> plot(y,x=0..20);
>
Błąd przybliżenia liczby Ą
> restart:
> xd:=Pi; # wartość doładna
xd := pð
> xp:=3.1415; # 3.1416 # wartość przybliżona
xp := 3.1415
def
%ð
Dðx =ð x -ð x
> Delta[x]:=xd-xp;
Dðx := pð ð-ð ð3.1415
def
Dðx
eðx =ð , x Ä…ð 0
x
> epsilon[x]:=Delta[x]/xd;
pð ð-ð ð3.1415
eðx :=
pð
1
eðx <ð 10-ðd
2
> solve(abs(epsilon[x])<1/2*10^(-d),{d});
{ d ð<ð ð4.229257627 }
> xp;evalf(Pi);
3.1415
3.141592654
>
BÅ‚Ä…d odejmowania bliskich sobie liczb
> Digits:=12; # 22
Digits := 12
Pierwsza liczba
> x1d:=sqrt(2);
x1p:=evalf[10](x1d);
Delta[x1]:=evalf(x1d-x1p);
epsilon[x1]:=evalf((x1d-x1p)/x1d);
x1d := 2
x1p := 1.414213562
Dðx1 := 0.37 10-9
eðx1 := 0.261629509038 10-9
Druga liczba
> x2d:=sqrt(2000001/1000000);
x2p:=evalf[10](x2d);
Delta[x2]:=evalf(x2d-x2p);
epsilon[x2]:=evalf((x2d-x2p)/x2d);
2000001
x2d :=
1000
x2p := 1.414213916
Dðx2 := -0.7 10-10
eðx2 := -0.494974623088 10-10
Błąd względny różnicy
Dðr Dðx1 -ð Dðx2
eðr =ð =ð
r x1 -ð x2
> epsilon[r]:=evalf((Delta[x1]-Delta[x2])/(x1d-x2d));
eðr := -0.00124448467021
Natomiast błąd względny sumy
Dðs Dðx1 +ð Dðx2
eðs =ð =ð
s x1 +ð x2
> epsilon[s]:=evalf((Delta[x1]+Delta[x2])/(x1d+x2d));
eðs := 0.106066003920 10-9
Jak zaradzić?
1. Zwiększyć Digits
>
> restart:
Propagacja błędów
> f:=(x1,x2)->4/3*x1*x2^3;
4
f := ( x1, x2 ) ð®ð ð x1 x23
3
2
ìð üð
Å›ð f
Dðf ð
íð żðDðx
åðîðÅ›ð xi %ð
x=ðx i
i=ð1
þð
> Delta[f]:=add((D[i](f)(x1p,x2p))*Delta[x||i],i=1..2);
4
Dðf := x2p3 Dðx1 ð+ð ð4 x1p x2p2 Dðx2
3
2
1 Å›ð f
eð ð Dðxi
f åðÅ›ð %ð
%ð
f (x) xi x=ðx
i=ð1
> epsilon[f]:=expand(Delta[f]/f(x1p,x2p));
Dðx1 3 Dðx2
eðf := ð+ð ð
x1p x2p
> x1p:=3.14;x2p:=5.00;
x1p := 3.14
x2p := 5.00
> f(x1p,x2p);
523.3333333
> Delta[x1]:=0.0016; Delta[x2]:=0.001;
Dðx1 := 0.0016
Dðx2 := 0.001
> 'Delta[f]'=Delta[f];
Dðf ð=ð ð0.5806666667
> 'epsilon[f]'=epsilon[f];
eðf ð=ð ð0.001109554140
>
Błędy obcięcia
> restart:
> f:=exp(-x);
(-ðx )
f := e
Rozwinięcie funkcji wokół x = 0
> sz:=taylor(f,x); # taylor(f,x,6);
1 1 1 1
sz := 1 ð-ð ðx ð+ð ð x2 ð-ð ð x3 ð+ð ð x4 ð-ð ð x5 ð+ð ðO( x6 )
2 6 24 120
> w:=convert(sz,polynom);
1 1 1 1
w := 1 ð-ð ðx ð+ð ð x2 ð-ð ð x3 ð+ð ð x4 ð-ð ð x5
2 6 24 120
> plot([f,w],x=0..2,0..1);
Wyznaczenie wartości funkcji dla x = 0.1
> x0:=0.1;
x0 := 0.1
> fp:=eval(w,x=x0);
fp := 0.9048374167
> fd:=eval(f,x=x0);
fd := 0.9048374180
> Delta['f']:=fd-fp;
Dðf := 0.13 10-8
> epsilon['f']:=Delta['f']/fd;
eðf := 0.1436722194 10-8
Wyznaczenie wartości funkcji w punkcie x = 2.1
> x0:=2.1;
x0 := 2.1
> fp:=eval(w,x=x0);
fp := 0.0314957500
> fd:=eval(f,x=x0);
fd := 0.1224564283
> Delta['f']:=fd-fp;
Dðf := 0.0909606783
> epsilon['f']:=Delta['f']/fd;
eðf := 0.7428003541
Błąd rzędu 74% !!!
>
Jak zaradzić?
1. Zwiększenie liczby wyrazów rozwinięcia (10)
2. Zmienić punkt rozwinięcia
> sz:=taylor(f,x,10);
1 1 1 1 1 1 1 1
sz := 1 ð-ð ðx ð+ð ð x2 ð-ð ð x3 ð+ð ð x4 ð-ð ð x5 ð+ð ð x6 ð-ð ð x7 ð+ð ð x8 ð-ð ð x9 ð+ð ðO( x10 )
2 6 24 120 720 5040 40320 362880
> w:=convert(sz,polynom);
1 1 1 1 1 1 1 1
w := 1 ð-ð ðx ð+ð ð x2 ð-ð ð x3 ð+ð ð x4 ð-ð ð x5 ð+ð ð x6 ð-ð ð x7 ð+ð ð x8 ð-ð ð x9
2 6 24 120 720 5040 40320 362880
> plot([f,w],x=0..4,0..1);
> fp:=eval(w,x=x0);
fp := 0.1220713254
> fd:=eval(f,x=x0);
fd := 0.1224564283
> Delta['f']:=fd-fp;
Dðf := 0.0003851029
> epsilon['f']:=Delta['f']/fd;
eðf := 0.003144815714
Błąd rzędu 0.3%
>
Zmiana punktu rozwinięcia (tylko 6 wyrazów)
> sz:=taylor(f,x=2);
1 1 1 -2 )
1
( -2 ) ( -2) (-2 ) (-2 ) ( (-2 )
sz := e ð-ð ðe ( x ð-ð ð2 ) ð+ð ð e ( x ð-ð ð2 )2 ð-ð ð e ( x ð-ð ð2 )3 ð+ð ð e ( x ð-ð ð2 )4 ð-ð ð e ( x ð-ð ð2 )5 ð+ð ðO( ( x ð-ð ð2 )6 )
2 6 24 120
> w:=convert(sz,polynom);
1 -2)
1 -2 )
1 -2 )
1 -2)
( -2 ) ( -2 ) ( ( ( (
w := e ð-ð ðe ( x ð-ð ð2 ) ð+ð ð e ( x ð-ð ð2 )2 ð-ð ð e ( x ð-ð ð2 )3 ð+ð ð e ( x ð-ð ð2 )4 ð-ð ð e ( x ð-ð ð2 )5
2 6 24 120
> plot([f,w],x=0..4,0..1);
> fp:=evalf(eval(w,x=x0));
fp := 0.1224564279
> fd:=eval(f,x=x0);
fd := 0.1224564283
> Delta['f']:=fd-fp;
Dðf := 0.4 10-9
> epsilon['f']:=Delta['f']/fd;
eðf := 0.3266467964 10-8
Wniosek: przyrost w szeregu Taylora powinien być mały
Błędy zaokrągleń
> restart:
Zagadka 1
> a1:=evalf(Pi)*10^5;
a1 := 314159.2654
> a2:=evalf(Pi)*10^(-5);
a2 := 0.00003141592654
> a3:=a1+a2;
a3 := 314159.2654
> is(a1=a3);
true
Jak zaradzić?
1. Zwiększyć Digits
> Digits:=20;
Digits := 20
> a3:=a1+a2;
a3 := 314159.26543141592654
> is(a1=a3);
false
> restart:
>
Zagadka 2
> restart:
> x:=987654321;
x := 987654321
> y:=evalf(Pi); # y:=Pi;
y := 3.141592654
> z:=-x;
z := -987654321
> #Digits:=20;
> x+y+z; # catastrophic cancellation!
3.1
> x+y;
0.9876543241 109
Jak zaradzić?
1. Zmienić kolejność w dodawaniu
2. Pozostawić (nie ewaluować) symbol Ą
3. Zmienić Digits
> Digits:=10;
Digits := 10
>
>
Zagadka 3
> restart:
> y:=1/1*(10^x*(1+(x)*10^(-x))-10^x);
( -ðx)
y := 10x ( 1 ð+ð ðx 10 ) ð-ð ð10x
> eval(x*10^(-x),x=20.);
0.2000000000 10-18
> 1+%;
1.
> eval(y,x=20.); # catastrophic cancellation!
0.
Jak zaradzić?
1. Uprościć wyrażenie
2. Podać argument w postaci liczby cakowitej
3. Zwiększyć Digits
> simplify(y);
x
> y;
( -ðx)
10x ( 1 ð+ð ðx 10 ) ð-ð ð10x
> eval(y,x=20); # argument całkowity
20
> evalf[25](eval(y,x=20.)); # wyższa dokładność
20.0000
>
Uwarunkowanie numeryczne wyznaczania zer wielomianów
> restart:
> n:=21; # n=20
n := 21
> w:=mul((x-i),i=1..n);
w := ( x ð-ð ð1 ) ( x ð-ð ð2 ) ( x ð-ð ð3 ) ( x ð-ð ð4 ) ( x ð-ð ð5 ) ( x ð-ð ð6 ) ( x ð-ð ð7 ) ( x ð-ð ð8 ) ( x ð-ð ð9 ) ( x ð-ð ð10 ) ( x ð-ð ð11 ) ( x ð-ð ð12 ) ( x ð-ð ð13 ) ( x ð-ð ð14 )
( x ð-ð ð15 ) ( x ð-ð ð16 ) ( x ð-ð ð17 ) ( x ð-ð ð18 ) ( x ð-ð ð19 ) ( x ð-ð ð20 ) ( x ð-ð ð21 )
> w:=expand(w);
w := -ð298631902863216384000 x2 ð+ð ð186244810780170240000 x ð-ð ð181664979520697076096 x4
ð+ð ð284093315901811468800 x3 ð+ð ð83637381699544802976 x5 ð-ð ð28939583397335447760 x6 ð+ð ð7744654310169576800 x7
ð-ð ð1634980697246583456 x8 ð+ð ð276019109275035346 x9 ð-ð ð37600535086859745 x10 ð+ð ð4154823851430525 x11
ð-ð ð373100999802531 x12 ð+ð ð27188611869881 x13 ð-ð ð1599718388730 x14 ð+ð ð75289668850 x15 ð-ð ð2792167686 x16 ð+ð ð79721796 x17
ð-ð ð1689765 x18 ð+ð ð25025 x19 ð-ð ð231 x20 ð+ð ðx21 ð-ð ð51090942171709440000
> wz:=w+10^(-9)*x^n;
wz := -ð298631902863216384000 x2 ð+ð ð186244810780170240000 x ð-ð ð181664979520697076096 x4
ð+ð ð284093315901811468800 x3 ð+ð ð83637381699544802976 x5 ð-ð ð28939583397335447760 x6 ð+ð ð7744654310169576800 x7
ð-ð ð1634980697246583456 x8 ð+ð ð276019109275035346 x9 ð-ð ð37600535086859745 x10 ð+ð ð4154823851430525 x11
ð-ð ð373100999802531 x12 ð+ð ð27188611869881 x13 ð-ð ð1599718388730 x14 ð+ð ð75289668850 x15 ð-ð ð2792167686 x16 ð+ð ð79721796 x17
1000000001
ð-ð ð1689765 x18 ð+ð ð25025 x19 ð-ð ð231 x20 ð+ð ð x21 ð-ð ð51090942171709440000
1000000000
> solve(w,x);
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21
> fsolve(wz,x);
1.000000000, 2.000000000, 3.000000000, 4.000000000, 4.999999999, 6.000000140, 6.999991102, 8.000294164,
8.994419597, 10.08469533, 10.65187913
> fsolve(wz,x,complex);
>
Uwarunkowanie numeryczne rozwiązywania liniowych ukladów równań
> restart:
> r1:=a1*x+a2*y+a3*z=a4;
r1 := a1 x ð+ð ða2 y ð+ð ða3 z ð=ð ða4
> r2:=b1*x+b2*y+b3*z=b4;
r2 := b1 x ð+ð ðb2 y ð+ð ðb3 z ð=ð ðb4
> r3:=c1*x+c2*y+c3*z=c4;
r3 := c1 x ð+ð ðc2 y ð+ð ðc3 z ð=ð ðc4
> M:=<,,>; # macierz współczynników
éða1 a2 a3Å‚ð
Ä™ð Å›ð
Ä™ð Å›ð
M := Ä™ðb1 b2 b3Å›ð
Ä™ð Å›ð
ëðc1 c2 c3ûð
> roz:=solve({r1,r2,r3},{x,y,z});
a2 b3 c4 ð-ð ða2 b4 c3 ð+ð ða3 c2 b4 ð-ð ða3 b2 c4 ð+ð ða4 b2 c3 ð-ð ða4 c2 b3
roz := { x ð=ð ð ,
a1 b2 c3 ð-ð ða1 c2 b3 ð-ð ðb1 a2 c3 ð+ð ðc2 b1 a3 ð-ð ðb2 c1 a3 ð+ð ðc1 a2 b3
a1 b3 c4 ð-ð ða1 b4 c3 ð-ð ðb1 a3 c4 ð-ð ðb3 c1 a4 ð+ð ðb4 c1 a3 ð+ð ðb1 a4 c3
y ð=ð ð-ð ,
a1 b2 c3 ð-ð ða1 c2 b3 ð-ð ðb1 a2 c3 ð+ð ðc2 b1 a3 ð-ð ðb2 c1 a3 ð+ð ðc1 a2 b3
c2 b1 a4 ð-ð ða1 c2 b4 ð+ð ða1 b2 c4 ð-ð ðb2 c1 a4 ð-ð ðb1 a2 c4 ð+ð ðc1 a2 b4
z ð=ð ð }
a1 b2 c3 ð-ð ða1 c2 b3 ð-ð ðb1 a2 c3 ð+ð ðc2 b1 a3 ð-ð ðb2 c1 a3 ð+ð ðc1 a2 b3
> assign(roz);
> x;
a2 b3 c4 ð-ð ða2 b4 c3 ð+ð ða3 c2 b4 ð-ð ða3 b2 c4 ð+ð ða4 b2 c3 ð-ð ða4 c2 b3
a1 b2 c3 ð-ð ða1 c2 b3 ð-ð ðb1 a2 c3 ð+ð ðc2 b1 a3 ð-ð ðb2 c1 a3 ð+ð ðc1 a2 b3
I zestaw danych
> a1:=1;a2:=2;a3:=3;a4:=1;b1:=4;b2:=5;b3:=6;b4:=2;c1:=7;c2:=8;c3:=8.99;c4:=1; #c3:=9.01;
a1 := 1
a2 := 2
a3 := 3
a4 := 1
b1 := 4
b2 := 5
b3 := 6
b4 := 2
c1 := 7
c2 := 8
c3 := 8.99
c4 := 1
> c3;
8.99
> x;
199.6666667
> LinearAlgebra[Determinant](M);
0.03
Zmienić c3
>
Dokładność softwarowa i sprzętowa
> Digits:=10; # wartość domyślna
Digits := 10
> evalf(Pi); # dokładność softwarowa
3.141592654
> evalhf(Pi); # dokładność sprzętowa
3.14159265358979312
> Digits:=50;
Digits := 50
> evalf(Pi); # dokładność softwarowa
3.1415926535897932384626433832795028841971693993751
> evalhf(Pi); # dokładność sprzętowa (18 cyfr znaczących)
3.14159265358979312
>
Zalecenia
Dodawać najpierw liczby mniejsze
> restart:
> Digits:=5;
Digits := 5
1000
S =ð i
åð
i=ð1
> S:=0:
> for i from 1000 by -1 to 1 do
S:=S+evalf(sqrt(i)):
end do:
> S;
21087.
> S:=0:
> for i to 1000 do
S:=S+evalf(sqrt(i)):
end do:
> S;
21094.
Jak zaradzić?
1. Zwiększyć Digits
2. Nie ewaluować pierwiastków
> Digits:=10;
Digits := 10
> S:=0:
> for i to 1000 do
S:=S+evalf(sqrt(i)):
end do:
> S; # wartość dokładniejsza
21097.45585
> Digits:=20;
Digits := 20
> S:=0:
> for i to 1000 do
S:=S+evalf(sqrt(i)):
end do:
> S; # wartość bardzo dokładna
21097.455887480735351
>
Podczas iteracji prowadzić obliczenia na liczbach zmiennoprzecinkowych
> restart:
> a:=10^6; n:=15;
a := 1000000
n := 15
> x:=1;
x := 1
> for i to 5 do
x:=1/2*(x+a/x);
end do;
1000001
x :=
2
1000006000001
x :=
4000004
1000028000070000028000001
x :=
8000056000056000008
1000120001820008008012870008008001820000120000001
x :=
16000560004368011440011440004368000560000016
1000496035960906202518364512465793311436201080861435825792904512250518300906192035960000496000001
x :=
32004960201379365884048929024827374165723285723067373729024508048803365856201376004960000032
>
> x:=1.;
x := 1.
> for i to n do
x:=1/2*(x+a/x);
end do;
x := 500000.5000
x := 250001.2500
x := 125002.6250
x := 62505.31242
x := 31260.65553
x := 15646.32231
x := 7855.117546
x := 3991.211544
x := 2120.881016
x := 1296.191593
x := 1033.841239
x := 1000.553871
x := 1000.000153
x := 1000.000000
x := 1000.000000
Próbować rozwiązywać problem w sposób dokładny
> restart:
e-ð x -ð x =ð 0
> fsolve(exp(-x)-x=0,x);
0.5671432904
> solve(exp(-x)-x=0,x);
LambertW( 1 )
> evalf[100](%);
0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946
> plot(exp(-x)-x,x=-infinity..infinity);
>