wykład 06 - Błędy Obliczeń, Inżynieria Oprogramowania - Informatyka, Semestr IV, Metody Obliczeniowe, ...
[ Pobierz całość w formacie PDF ]
>
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
10
()
x
10
x
(
10
x
y
:=
1
)
>
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
:=
>
xp:=3.1415; # 3.1416 # wartość przybliżona
:=
xp
3.1415
def
x xx
>
Delta[x]:=xd-xp;
x
:=
3.1415
x
def
,
x
0
x
x
>
epsilon[x]:=Delta[x]/xd;
3.1415
x
:=
1
10
2
d
x
>
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
:=
0.37 10
-9
x1
0.261629509038 10
-9
x1
:=
Druga liczba
>
x2d:=sqrt(2000001/1000000);
x2p:=evalf[10](x2d);
Delta[x2]:=evalf(x2d-x2p);
epsilon[x2]:=evalf((x2d-x2p)/x2d);
2000001
1000
x2d
:=
x2p
:=
1.414213916
:=
-0.7 10
-10
x2
-0.494974623088 10
-10
x2
:=
Błąd względny różnicy
>
epsilon[r]:=evalf((Delta[x1]-Delta[x2])/(x1d-x2d));
:=
r
x
x
1
2
r
r
x
x
1
2
r
-0.00124448467021
Natomiast błąd względny sumy
x
x
s
s xx
>
epsilon[s]:=evalf((Delta[x1]+Delta[x2])/(x1d+x2d));
:=
1
2
s
1
2
0.106066003920 10
-9
s
Jak zaradzić?
1. Zwiększyć Digits
>
>
restart:
Propagacja błędów
>
f:=(x1,x2)->4/3*x1*x2^3;
4
3
x1 x2
3
f
:=
(
x1 x2
,
)
2
f
f
x
x
i
x
i
1
i
>
Delta[f]:=add((D[i](f)(x1p,x2p))*Delta[x||i],i=1..2);
4
3
x2p
3
x1
4
x1p x2p
2
x2
f
:=
1
()
2
f
x
f
xx
i
>
epsilon[f]:=expand(Delta[f]/f(x1p,x2p));
f
x
x
i
1
i
x1
x1p
3
x2
x2p
f
:=
>
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;
x1
:=
0.0016
:=
0.001
x2
>
'Delta[f]'=Delta[f];
0.5806666667
f
>
'epsilon[f]'=epsilon[f];
0.001109554140
f
>
Błędy obcięcia
>
restart:
>
f:=exp(-x);
e
()
x
f
:=
Rozwinięcie funkcji wokół x = 0
>
sz:=taylor(f,x); # taylor(f,x,6);
1
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
O
x
6
sz
:=
1
x
()
>
w:=convert(sz,polynom);
1
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
w
:=
1
x
>
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;
0.13 10
-8
f
:=
>
epsilon['f']:=Delta['f']/fd;
0.1436722194 10
-8
f
:=
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;
f
:=
0.0909606783
>
epsilon['f']:=Delta['f']/fd;
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
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
1
720
x
6
1
5040
x
7
1
40320
x
8
1
362880
x
9
O
x
10
sz
:=
1
x
(
)
>
w:=convert(sz,polynom);
1
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
1
720
x
6
1
5040
x
7
1
40320
x
8
1
362880
x
9
w
:=
1
x
>
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;
f
:=
0.0003851029
>
epsilon['f']:=Delta['f']/fd;
f
:=
0.003144815714
Błąd rzędu 0.3%
>
Zmiana punktu rozwinięcia (tylko 6 wyrazów)
>
sz:=taylor(f,x=2);
1
2
e
()
1
6
e
()
1
24
e
()
1
120
e
()
e
()
-2
e
()
-2
-2
-2
-2
-2
2
3
4
5
6
sz
:=
(
x
2
)
(
x
2
)
(
x
2
)
(
x
2
)
(
x
2
)
O(
(
x
2
)
)
>
w:=convert(sz,polynom);
1
1
1
1
e
()
-2
e
()
-2
2
e
()
-2
6
e
()
-2
24
e
()
-2
120
e
()
-2
2
3
4
5
w
:=
(
x
2
)
(
x
2
)
(
x
2
)
(
x
2
)
(
x
2
)
>
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;
0.4 10
-9
f
:=
>
epsilon['f']:=Delta['f']/fd;
0.3266467964 10
-8
f
:=
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 10
9
Jak zaradzić?
1. Zmienić kolejność w dodawaniu
2. Pozostawić (nie ewaluować) symbol
π
3. Zmienić Digits
>
Digits:=10;
[ Pobierz całość w formacie PDF ]