C2.mw

数値計算:第2章 誤差

関西学院大学・理工学部・情報科学科 西谷滋人

Significant digits(有効桁数)

> restart;

1ワードの整数の最大値とその2進数表示.

> 2^(4*8-1)-1;
convert(2^(4*8-1)-1,binary);

2147483647

1111111111111111111111111111111

1ワードの整数の最大桁.

> length(2^(4*8-1)-1);

10

64bitの場合の整数の最大桁.

> length(2^(8*8-1)-1);

19

Mapleでは多倍長計算するので,通常のプログラミング言語で起こるintの最大数あたりでの

奇妙な振る舞いは示さない.

> 2147483647+100;

2147483747

> ?length;

1ワードの実数では,仮数部2進数23bit,二倍長実数で54bit.有効桁数は以下の通り.

> length(2^(23));
length(2^(52));

7

16

Quadratic equation(二次方程式)

> restart;

ケタ落ちを起こしそうな,二次方程式の係数を指定.

> a:=1;b:=10000000;c:=1;

a := 1

b := 10000000

c := 1

解をsolveで解析的に求め,その後,浮動小数点に直す.

デフォルトでは有効桁数が10桁なので思った答えを返さない.

大きな有効桁数を指定すると精確な値が返ってくる.

> s1:=solve(a*x^2+b*x+c=0,x);
evalf(s1);

evalf(s1,40);

s1 := -5000000+24999999999999^(1/2), -5000000-24999999999999^(1/2)

0., -10000000.00

-0.100000000000001000000000000e-6, -9999999.999999899999999999999000000000000

数値的に解を求めるfsolveを使うと有効数字10桁で正しい答えを返す.

> fsolve(a*x^2+b*x+c=0,x);

-10000000.00, -0.1000000000e-6

fsolveのoptionにあるfulldigitsを指定しても解は改善されない.

> fsolve(a*x^2+b*x+c=0,x,fulldigits);

-10000000.00, -0.1000000000e-6

大域パラメータのDigits(Defaultは10)を大きくすると,fsolveでも精確な解を求めてくれる.

> Digits:=40;
(x1,x2):=fsolve(a*x^2+b*x+c=0,x);

Digits := 40

x1, x2 := -9999999.999999899999999999999000000000000, -0.1000000000000010000000000000200000000000e-6

printfで見やすく表示させる.

> printf("%45.40e \n%25.20e \n%25.20e\n",x1,x1,x2);

-9.9999999999998999999999999990000000000000e+06
-9.99999999999990000000e+06

-1.00000000000001000000e-07

その他補足

> a:=100001;
b:=100000;

suma:=0;

sumb:=0;

sum1:=0;

for i from 1 to 100000 do

 suma:=suma+a;

 sumb:=sumb+b;

 sum1:=sum1+a-b;

end do:

sum1;

suma-sumb;

a := 100001

b := 100000

suma := 0

sumb := 0

sum1 := 0

100000

100000

> normal(1/(1-x)-1/(1+x));

-2*x/((-1+x)*(1+x))

> combine(sin(x)^2,trig);

1/2-1/2*cos(2*x)