>
restart:
Przykład 1.
>
a1:=evalf(Pi)*10^5;
:=
a1
314159.2654
>
a2:=evalf(Pi)*10^(-5);
:=
a2
0.00003141592654
>
a3:=a1+a2:
>
is(a1=a3);
true
Przykład 2.
>
y:=1/1*(10^x*(1+(x)*10^(-x))-10^x);
:=
y
−
10
x
(
)
+
1 x 10
(
)
−x
10
x
>
eval(y,x=10.5);
0.
>
simplify(eval(y,x=21/2));
21
2
>
>
restart:
Błąd przybliżenia liczby
π
>
xs:=Pi;
:=
xs
π
>
xp:=3.1415; # 3.1416
:=
xp
3.1415
>
Delta[x]:=xs-xp;
:=
∆
x
−
π 3.1415
>
epsilon[x]:=Delta[x]/xs;
:=
ε
x
−
π 3.1415
π
>
solve(abs(epsilon[x])<1/2*10^(-d),d);
(
)
RealRange
,
−∞
(
)
Open 4.229257627
>
evalf(Pi);
3.141592654
>
Błąd odejmowania bliskich sobie liczb
>
Digits:=12;
:=
Digits
12
Pierwsza liczba
>
x1s:=sqrt(2);
:=
x1s
2
>
x1p:=evalf[10](x1s);
:=
x1p
1.414213562
>
Delta[x1]:=evalf(x1s-x1p);
:=
∆
x1
0.37 10
-9
>
epsilon[x1]:=evalf((x1s-x1p)/x1s);
:=
ε
x1
0.261629509038 10
-9
Druga liczba
>
x2s:=sqrt(201/100);
x2p:=evalf[10](x2s);
Delta[x2]:=evalf(x2s-x2p);
epsilon[x2]:=evalf((x2s-x2p)/x2s);
:=
x2s
201
10
:=
x2p
1.417744688
:=
∆
x2
-0.12 10
-9
:=
ε
x2
-0.846414739035 10
-10
Błąd względny różnicy
2
1
2
1
x
x
x
x
r
r
r
−
∆
−
∆
=
∆
=
ε
>
epsilon[r]:=evalf((Delta[x1]-Delta[x2])/(x1s-x2s));
:=
ε
r
-0.138765953975 10
-6
Natomiast błąd względny sumy
2
1
2
1
x
x
x
x
s
s
s
+
∆
+
∆
=
∆
=
ε
>
epsilon[s]:=evalf((Delta[x1]+Delta[x2])/(x1s+x2s));
:=
ε
s
0.882781375672 10
-10
>
restart:
Propagacja błędów
>
f:=(x1,x2)->4/3*x1*x2^3;
:=
f
→
(
)
,
x1 x2
4
3
x1 x2
3
2
1
i
i
i
f
f
x
x
=
=
∂
∆ ≈
∆
∂
∑
x x
>
Delta[f]:=add((D[i](f)(x1p,x2p))*Delta[x||i],i=1..2);
:=
∆
f
+
4
3
x2p
3
∆
x1
4 x1p x2p
2
∆
x2
2
1
1
( )
f
i
i
i
f
x
f
x
=
=
∂
ε ≈
∆
∂
∑
x x
x
>
epsilon[f]:=expand(Delta[f]/f(x1p,x2p));
:=
ε
f
+
∆
x1
x1p
3
∆
x2
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;
:=
∆
x1
0.0016
:=
∆
x2
0.001
>
'Delta[f]'=Delta[f];
=
∆
f
0.5806666667
>
'epsilon[f]'=epsilon[f];
=
ε
f
0.001109554140
>
Błędy obcięcia
>
restart:
>
f:=exp(-x);
:=
f
e
(
)
−x
Rozwinięcie funkcji wokół x = 0
>
sz:=taylor(f,x); # taylor(f,x,6);
:=
sz
− +
−
+
−
+
1 x
1
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
( )
O x
6
>
w:=convert(sz,polynom);
:=
w
− +
−
+
−
1 x
1
2
x
2
1
6
x
3
1
24
x
4
1
120
x
5
>
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;
:=
∆
f
0.13 10
-8
>
epsilon['f']:=Delta['f']/fd;
:=
ε
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;
:=
∆
f
0.0909606783
>
epsilon['f']:=Delta['f']/fd;
:=
ε
f
0.7428003541
Błąd rzędu 74% !!!
>
Zwiększenie liczby wyrazów rozwinięcia (10)
>
sz:=taylor(f,x,10);
:=
sz
− +
−
+
−
+
−
+
−
+
1 x
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
>
w:=convert(sz,polynom);
:=
w
− +
−
+
−
+
−
+
−
1 x
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
>
plot([f,w],x=0..2,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);
sz
e
( )
-2
e
( )
-2
(
)
−
x 2
1
2
e
( )
-2
(
)
−
x 2
2
1
6
e
( )
-2
(
)
−
x 2
3
1
24
e
( )
-2
(
)
−
x 2
4
1
120
e
( )
-2
−
+
−
+
−
:=
(
)
−
x 2
5
(
)
O (
)
−
x 2
6
+
>
w:=convert(sz,polynom);
w :=
−
+
−
+
−
e
( )
-2
e
( )
-2
(
)
−
x 2
1
2
e
( )
-2
(
)
−
x 2
2
1
6
e
( )
-2
(
)
−
x 2
3
1
24
e
( )
-2
(
)
−
x 2
4
1
120
e
( )
-2
(
)
−
x 2
5
>
plot([f,w],x=0..2,0..1);
>
fp:=evalf(eval(w,x=x0));
:=
fp
0.1224564279
>
fd:=eval(f,x=x0);
:=
fd
0.1224564283
>
Delta['f']:=fd-fp;
:=
∆
f
0.4 10
-9
>
epsilon['f']:=Delta['f']/fd;
:=
ε
f
0.3266467964 10
-8
Wniosek: przyrost w szeregu Taylora powinien być mały
Błędy zaokrągleń
>
restart:
>
x1d:=1.12345678987654321234567898765432123456789;
:=
x1d
1.12345678987654321234567898765432123456789
>
x1p:=1.*x1d;
:=
x1p
1.123456790
>
Delta['x1']:=evalf[12](x1d-x1p);
:=
∆
x1
-0.12 10
-9
>
x2d:=1/6;
:=
x2d
1
6
>
x2p:=evalf(1/6);
:=
x2p
0.1666666667
>
Delta['x2']:=evalf[11](x2d-x2p);
:=
∆
x2
-0.3 10
-10
>
x3d:=sqrt(2);
:=
x3d
2
>
x3p:=evalf(sqrt(2));
:=
x3p
1.414213562
>
Delta['x3']:=evalf[11](x3d-x3p);
:=
∆
x3
0.4 10
-9
Przykłady katastroficznych redukcji
Przykład 1
>
x:=987654321;
:=
x
987654321
>
y:=evalf(Pi); # pozostawic symbol !!!
:=
y
3.141592654
>
z:=-x;
:=
z
-987654321
>
#Digits:=20;
>
x+y+z; # catastrophic cancellation
3.1
>
x+y; # zmienic Digits dwa wiersze wyzej
0.9876543241 10
9
>
Digits:=10;
:=
Digits
10
Przykład 2
>
x:='x':
>
y:=1/1*(10^x*(1+(x)*10^(-x))-10^x);
:=
y
−
10
x
(
)
+
1 x 10
(
)
−x
10
x
>
eval(y,x=10.5);
0.
>
eval(x*10^(-x),x=10.5);
0.3320391543 10
-9
>
1+%;
1.000000000
>
evalf[25](eval(y,x=10.5));
10.50000000000001
>
for i from 10 to 15 do evalf[i](eval(y,x=10.5)); end do;
0.
9.
10.4
10.50
10.499
10.5000
>
y;
−
10
x
(
)
+
1 x 10
(
)
−x
10
x
Jak zaradzic?
>
eval(y,x=21/2);
−
10000000000 10
+
1
21 10
200000000000
10000000000 10
>
simplify(%);
21
2
>
evalf(%);
10.50000000
>
Uwarunkowanie numeryczne wyznaczania zer wielomianów
>
restart:
>
n:=10; # n=20
:=
n
10
>
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
>
w:=expand(w);
w
x
10
55 x
9
1320 x
8
18150 x
7
157773 x
6
902055 x
5
3416930 x
4
8409500 x
3
−
+
−
+
−
+
−
:=
12753576 x
2
10628640 x 3628800
+
−
+
>
wz:=w+10^(-9)*x^n;
wz
1000000001
1000000000
x
10
55 x
9
1320 x
8
18150 x
7
157773 x
6
902055 x
5
3416930 x
4
−
+
−
+
−
+
:=
8409500 x
3
12753576 x
2
10628640 x 3628800
−
+
−
+
>
solve(w,x);
, , , , , , , , ,
1 2 3 4 5 6 7 8 9 10
>
fsolve(wz,x);
1.000000000 2.000000000 3.000000006 3.999999757 5.000003391 5.999979005
,
,
,
,
,
,
7.000065391 7.999893480 9.000086473 9.999972441
,
,
,
>
fsolve(wz,x,complex);
1.000000000 2.000000000 3.000000006 3.999999757 5.000003391 5.999979005
,
,
,
,
,
,
7.000065391 7.999893480 9.000086473 9.999972441
,
,
,
>
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:=<<a1|a2|a3>,<b1|b2|b3>,<c1|c2|c3>>; # macierz współczynników
:=
M
a1 a2 a3
b1 b2 b3
c1 c2 c3
>
roz:=solve({r1,r2,r3},{x,y,z});
roz
=
z
−
+
−
−
+
c2 b1 a4 a1 c2 b4 a1 b2 c4 b2 c1 a4 b1 a2 c4 c1 a2 b4
−
−
+
−
+
a1 b2 c3 a1 c2 b3 b1 a2 c3 c2 b1 a3 b2 c1 a3 c1 a2 b3
,
{
:=
=
y
−
−
−
−
+
+
a1 b3 c4 a1 b4 c3 b1 a3 c4 b3 c1 a4 b4 c1 a3 b1 a4 c3
−
−
+
−
+
a1 b2 c3 a1 c2 b3 b1 a2 c3 c2 b1 a3 b2 c1 a3 c1 a2 b3
,
=
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
}
>
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
>
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
9.01
>
c3;
9.01
>
x;
-200.3333333
>
LinearAlgebra[Determinant](M);
-0.03
Zmienić c3
>
a1:=1;a2:=3;a3:=5;a4:=1;b1:=2;b2:=1;b3:=6;b4:=2;c1:=3;c2:=-5;c3:
=8.99;c4:=1; c3:=9.01;
:=
a1
1
:=
a2
3
:=
a3
5
:=
a4
1
:=
b1
2
:=
b2
1
:=
b3
6
:=
b4
2
:=
c1
3
:=
c2
-5
:=
c3
8.99
:=
c4
1
:=
c3
9.01
>
x;
1.998080614
>
LinearAlgebra[Determinant](M);
-26.05
Zminenić c3
Dokladność softwarowa i sprzętowa
>
Digits:=10; # wartosc domyślna
:=
Digits
10
>
evalf(Pi);
3.141592654
>
evalhf(Pi);
3.14159265358979312
>
Digits:=50;
:=
Digits
50
>
evalf(Pi);
3.1415926535897932384626433832795028841971693993751
>
evalhf(Pi); # 16 cyfr znaczących
3.14159265358979312
>
Zalecenia
Dodawać najpierw liczby najmniejsze
>
restart:
>
Digits:=5;
:=
Digits
5
>
S:=0:
for i from 10^3 by -1 to 1 do
S:=S+evalf(sqrt(i)):
end do:
S;
21087.
>
S:=0:
for i to 10^3 do
S:=S+evalf(sqrt(i)):
end do:
S;
21094.
>
evalf[10](add(sqrt(i),i=1..1000)); # wartość dokładniejsza
21097.45590
Podczas iteracji prowadzić obliczenia na liczbach zmiennoprzecinkowych
>
restart:
>
Digits:=10:
>
a:=10^6; n:=15;
:=
a
1000000
:=
n
15
>
x:=1;
:=
x
1
>
for i to n do
x:=1/2*(x+a/x);
end do:
>
evalf(x);
1000.000000
>
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
>