数値解析 第1回 (4) 誤差の発生メカニズム
このページでは、コンピュータが誤差を産み出す原因をいくつか紹介します。
1. 丸め誤差
C言語や Java などのプログラミング言語で実数を表すデータ型には float 型と double 型があり、
多くの処理系で float 型は 4 バイトデータ、double 型は 8 バイトデータです。
IEEE規格 ( IEEE754 ) における double 型 8 バイト （ = 64 ビット ）の内訳は
- 符号 : 1 ビット
- 指数部 : 11 ビット
- 仮数部 : 52 ビット
です。
問 64ビットのデータは何通りありますか？
答え
$2^{64} = 18446744073709551616$ 通り。
問 64ビットのデータは何通りありますか？
答え
$2^{64} = 18446744073709551616$ 通り。
本当は無限個存在する実数を、
コンピュータ上ではこの有限通りのどれかに当てはめて表現しますので必然的に誤差を生じます。
これを「丸め誤差」と呼びます。
「丸める」の英語は round、10 進数なら四捨五入のことで、Excel の関数にもありましたね。
Ex.1
$e^{\pi\sqrt{163}}$ は、本当は無理数
$262537412640768743.999999999999250072597198185688879\cdots$
ですが、充分な精度をもって計算しないと
$262537412640768744$
という整数のように勘違いしています。
# -*- coding: utf-8 -*-
# exppisqrt163.py
from decimal import *
for p in range(10, 61, 5):
getcontext().prec = p
pi = Decimal('3.141592653589793238462643383279502884197169399(以下略)')
x = Decimal('163')
y = (pi * x.sqrt()).exp()
print(u'精度 =', p, ": ", y)
精度 = 10 : 2.625374097E+17
精度 = 15 : 2.62537412640764E+17
精度 = 20 : 262537412640768744.17
精度 = 25 : 262537412640768743.9999990
精度 = 30 : 262537412640768744.000000000024
精度 = 35 : 262537412640768743.99999999999924980
精度 = 40 : 262537412640768743.9999999999992500725944
精度 = 45 : 262537412640768743.999999999999250072597198209
精度 = 50 : 262537412640768743.99999999999925007259719818568865
精度 = 55 : 262537412640768743.9999999999992500725971981856888793537
精度 = 60 : 262537412640768743.999999999999250072597198185688879353856337
2. オーバーフロー・アンダーフロー
- 表現できる最大の数を超えてしまうエラーをオーバーフローと言います。
- 表現できる最小の数を下回って 0 に置き換えられてしまう現象をアンダーフローと言います。
float 型、double 型のオーバーフローはプログラムがエラー表示を出してくれますのでまだましです。
アンダーフローや、int 型のオーバーフローは気付かないところで起こっても
そのまま実行し続けるのでタチが悪いです。
Ex.2
$10$ のべき乗を double 型で計算してゆくと、やがてオーバーフローが起こり、
inf ( 無限大 = infinity を表す記号) が表示されます。
/* error_overflow.c */
#include<stdio.h>
int main()
{
int i;
double x = 1;
for(i = 1; i <= 311; i++){
x *= 10;
printf("10^%d = %e\n", i, x);
}
return 0;
}
(中略)
10^306 = 1.000000e+306
10^307 = 1.000000e+307
10^308 = 1.000000e+308
10^309 = inf
10^310 = inf
10^311 = inf
Ex.3
$2$ のべき乗を 4 バイトの int 型で計算するとこうなります。
/* error_int_overflow1.c */
#include<stdio.h>
int main()
{
int i, x = 1;
for(i = 1; i <= 33; i++){
x *= 2;
printf("2^%2d = %d\n", i, x);
}
return 0;
}
(中略)
2^29 = 536870912
2^30 = 1073741824
2^31 = -2147483648
2^32 = 0
2^33 = 0
理由
int 型で $2$ を掛けることは、ビット列を左シフトすることです。
従って $2$ の $31$ 乗はビット列として
10000000000000000000000000000000
であり、「$2$ の補数」としてこれは $-2^{31}$ のことになります。
さらに $2$ の $32$ 乗はそれを左シフトして
00000000000000000000000000000000
になります。
Ex.4
今度は $3$ のべき乗を int 型で計算してみましょう。
/* error_int_overflow2.c */
#include<stdio.h>
int main()
{
int i, x = 1;
for(i = 1; i <= 24; i++){
x *= 3;
printf("3^%2d = %d\n", i, x);
}
return 0;
}
(中略)
3^18 = 387420489
3^19 = 1162261467
3^20 = -808182895
3^21 = 1870418611
3^22 = 1316288537
3^23 = -346101685
3^24 = -1038305055
$20$ 乗以上はオーバーフローが起こってケッタイなことになっています。
Ex.5
$10^{-n}$ に $10^n$ を掛けると理論上は $1$ に戻るはずですが、
$n$ が大きくなると $10^{-n}$ がアンダーフローを起こしてやはりケッタイなことになります。
/* error_underflow.c */
#include<stdio.h>
int main()
{
int i, n;
float x;
for(n = 38; n <= 47; n++){
x = 1.0;
for(i = 1; i <= n; i++) x /= 10;
for(i = 1; i <= n; i++) x *= 10;
printf("n = %d のとき %f\n", n, x);
}
return 0;
}
n = 38 のとき 1.000000
n = 39 のとき 1.000000
n = 40 のとき 0.999995
n = 41 のとき 0.999967
n = 42 のとき 1.000527
n = 43 のとき 0.994922
n = 44 のとき 0.980909
n = 45 のとき 1.401299
n = 46 のとき 0.000000
n = 47 のとき 0.000000
3. 累積誤差
いわゆる「塵も積もれば山となる」というやつです。
Ex.6
$1/n$ を $n$ 個足すと $1$ のはずですが ...
/* error_accumulation3.c */
#include<stdio.h>
int main()
{
int i, j, n = 1;
double x, sum;
for(i = 0; i < 13; i++){
n *= 5;
x = 1.0 / n;
sum = 0.0;
for(j = 0; j < n; j++) sum += x;
printf("n = %10d のとき 和 = %13.10lf\n", n, sum);
}
return 0;
}
(中略)
n = 1953125 のとき 和 = 1.0000000000
n = 9765625 のとき 和 = 1.0000000001
n = 48828125 のとき 和 = 1.0000000004
n = 244140625 のとき 和 = 0.9999999973
n = 1220703125 のとき 和 = 1.0000000227
4. 積み残し
絶対値が大きく異なる2つの数を足したり引いたりすると、
小さい方の数の下位の桁が無視される現象です。
例えば、ともに有効数字が5桁のふたつの数 $x=1.2345$ と $y=0.0056789$ を足すとき、
$x$ は最大で $5 \times 10^{-5}$ の絶対誤差を持っていますので、
$y$ の下位3桁は書いても意味がない、ということです。
$\begin{array}{cl}
&1.2345 \\
+&0.0056789 \\
\hline
&1.2402\\
\end{array}$
Ex.7
$\displaystyle{\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+\cdots}$ の
理論値は $\displaystyle{\frac{\pi^2}{6} = 1.6449340668 \cdots}$ ですが ...
/* error_zeta2.c */
#include<stdio.h>
#include<math.h>
#define pi 3.14159265358979324
int main()
{
int i, j, n;
float x;
printf("理論値は %f だが\n\n", pi * pi / 6);
printf("昇順に加えると\n");
for(j = 2; j <= 7; j++){
n = (int)pow(10, j);
x = 0.0;
for(i = 1; i <= n; i++) x += 1.0 / i / i;
printf("%8d項目までの和 = %f\n", n, x);
}
printf("\n");
printf("降順に加えるとましで\n");
for(j = 2; j <= 7; j++){
n = (int)pow(10, j);
x = 0.0;
for(i = n; i >= 1; i--) x += 1. / i / i;
printf("%8d項目までの和 = %f\n", n, x);
}
return 0;
}
理論値は 1.644934 だが
昇順に加えると
100項目までの和 = 1.634984
1000項目までの和 = 1.643935
10000項目までの和 = 1.644725
100000項目までの和 = 1.644725
1000000項目までの和 = 1.644725
10000000項目までの和 = 1.644725
降順に加えるとましで
100項目までの和 = 1.634984
1000項目までの和 = 1.643934
10000項目までの和 = 1.644834
100000項目までの和 = 1.644924
1000000項目までの和 = 1.644933
10000000項目までの和 = 1.644934
昇順に足すと $1$ 以上ある部分和に対して小さな $1/(n^2)$ を足すことになるので積み残しが多いのに対し、
降順に足すと、まだ比較的小さい部分和に対して $1/(n^2)$ を足してゆくので積み残しが少なくて済みます。
5. 桁落ち
値の近いふたつの数を引くと、有効数字の桁数が減ってしまう現象です。例えば
$\begin{array}{cr}
&\sqrt{10}=3.1622 \\
-&\pi=3.1416 \\
\hline
&0.0202\\
\end{array}$
Ex.8
$\displaystyle{\left(\frac{1}{n} - \frac{1}{n+1}\right)\times n \times (n+1) = 1}$ のはずですが ...
/* error_cancellation1.c */
#include<stdio.h>
int main()
{
int n, i;
float x, y, z, w;
n = 10;
for(i = 0; i < 8; i++){
x = 1.0 / n;
y = 1.0 / (n + 1);
z = x - y;
w = z * (n * (n + 1.0));
printf("n = %9d のとき [1/n-1/(n+1)]*(n*(n+1)) = %f\n", n, w);
n *= 10;
}
return 0;
}
n = 10 のとき [1/n-1/(n+1)]*(n*(n+1)) = 1.000000
n = 100 のとき [1/n-1/(n+1)]*(n*(n+1)) = 0.999999
n = 1000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 1.000075
n = 10000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 0.999817
n = 100000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 1.000454
n = 1000000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 1.023183
n = 10000000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 1.421086
n = 100000000 のとき [1/n-1/(n+1)]*(n*(n+1)) = 0.000000
$1/n$ と $1/(n+1)$ が近い値なので桁落ちが起こっています。