> 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);
( -x )
y := 10x (1 + x 10 ) - 10x
> 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;
Ä„ - 3.1415
µx :=
Ä„
> 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);
201
x2s :=
10
x2p := 1.417744688
"x2 := -0.12 10-9
µx2 := -0.846414739035 10-10
Błąd względny różnicy
"r "x1 - "x2
µr = =
r x1 - x2
> epsilon[r]:=evalf((Delta[x1]-Delta[x2])/(x1s-x2s));
µr := -0.138765953975 10-6
Natomiast błąd względny sumy
"s "x1 + "x2
µs = =
s x1 + x2
> 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;
4
f := (x1, x2 ) x1 x23
3
2
Å„Å‚üÅ‚
" f
"f H"
òłżł"x
"ół" xi
x=x i
i=1
þÅ‚
> Delta[f]:=add((D[i](f)(x1p,x2p))*Delta[x||i],i=1..2);
4
"f := x2p3 "x1 + 4 x1p x2p2 "x2
3
2
1 " f
µ H" "xi
f "" xi
x=x

f (x)i=1
> epsilon[f]:=expand(Delta[f]/f(x1p,x2p));
"x1 3 "x2
µ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;
"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);
( -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;
"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);
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..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);
1 -2 )
1 -2 )
1 -2 )
1 -2 )
( -2 ) ( -2 ) ( ( ( (
sz := e - e (x - 2 ) + e (x - 2 )2 - e ( x - 2)3 + e (x - 2 )4 - e
2 6 24 120
( x - 2 )5 + O( ( x - 2 )6 )
> w:=convert(sz,polynom);
w :=
1 1 -2 )
1 1 -2 )
( -2 ) (-2 ) (-2 ) ( (-2 ) (
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..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;
1
x2d :=
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 109
> Digits:=10;
Digits := 10
Przykład 2
> x:='x':
> y:=1/1*(10^x*(1+(x)*10^(-x))-10^x);
( -x )
y := 10x (1 + x 10 ) - 10x
> 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;
( -x )
10x (1 + x 10 ) - 10x
Jak zaradzic?
> eval(y,x=21/2);
ëÅ‚ öÅ‚
21 10
ìÅ‚ ÷Å‚
10000000000 10 1 +
ìÅ‚ ÷Å‚ - 10000000000 10
íÅ‚ 200000000000 Å‚Å‚
> 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 := x10 - 55 x9 + 1320 x8 - 18150 x7 + 157773 x6 - 902055 x5 + 3416930 x4 - 8409500 x3
+ 12753576 x2 - 10628640 x + 3628800
> wz:=w+10^(-9)*x^n;
1000000001
wz := x10 - 55 x9 + 1320 x8 - 18150 x7 + 157773 x6 - 902055 x5 + 3416930 x4
1000000000
- 8409500 x3 + 12753576 x2 - 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:=<,,>; # macierz współczynników
îÅ‚a1 a2 a3Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
M := ïÅ‚b1 b2 b3śł
ïÅ‚ śł
ðÅ‚c1 c2 c3ûÅ‚
> roz:=solve({r1,r2,r3},{x,y,z});
c2 b1 a4 - a1 c2 b4 + a1 b2 c4 - b2 c1 a4 - b1 a2 c4 + c1 a2 b4
roz := { z = ,
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
a2 b3 c4 - a2 b4 c3 + a3 c2 b4 - a3 b2 c4 + a4 b2 c3 - a4 c2 b3
x = }
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
>