情報科学概論 9 回目
2015/12/3
今回の目的
• 数値計算法 その 1 –誤差，方程式の解法–
数値の表現と限界
1
1.1
整数型
• 整数は 2 進数で厳密に表現可能
⇒ 演算結果が整数のものについては，誤差は気にしなくて良い．
• 但し，オーバーフローに注意．表現できる数字の範囲には制限がある．
整数型に入りきれない大きな（小さな）数字を扱う時
−215 ≤ short int ≤ 215 − 1
1.2
浮動小数点型
• 実数は 2 進数で厳密に表現しきれない．⇒ IEEE 形式の浮動小数点（表現方法の一つ）
符号部と指数部と仮数部に分け，2 進数で表わす．
• float 型の場合
4 bytes = 32 bits を 以下のように分ける．
符号部 (S)1 bit, 指数部 (E) 8 bits, 仮数部 (M) 23 bits ⇒ (−1)S × 2E−127 × 1.M
(例): S = 1, E = 10000001, M = 01000000000000000000000 ⇒ (−1)1 × 2129−127 × (1.25) = −5.0
範囲: 1.175494 × 10−38 < |float| < 3.402823 × 10+38
• double 型の場合
8 bytes = 64 bits を符号部 (S)1 bit, 指数部 (E) 11 bits, 仮数部 (M) 52 bits に分ける．
範囲: 2.225074 × 10−308 < |double| < 1.797693 × 10+308
数値計算時の誤差
2
2.1
丸め誤差
• 厳密に表現できない数字は「丸め」られる．
「丸め」: 近似値に変換すること
• 10 進数では有限小数でも 2 進数では無限小数になる場合がある．
(例)10 進数の 0.1 は 2 進数では
0.0001100110011001100110011001... = +1.10011... × 2−4
だから float で表すと S = 0, E = 123 となるので，
0 01111011 10011001100110011001100 [11001100 の循環]
つまり，仮数部の 24 ビット目以下が切り捨てられるので，
0.00000000000000000000000110011001100 × 2−4 ≃ 6 × 10−9 だけ誤差がでたことになる．
• float 型では仮数部の 24 bits 目以降は切り捨て．⇒ 最初の桁の 2−23 ∼ 10−7 倍程度までが信用できる．
• double 型では仮数部の 53 bits 目以降は切り捨て．⇒ 最初の桁の 2−52 ∼ 2 × 10−16 倍程度までが信用で
きる．
• 例
1
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
int i, n = 100;
float x, sum, sum2;
x = 1.0 / n; /* 0.01 を計算しているつもり */
sum = 0;
for (i = 0; i < n; i++){
sum += x;
}
printf("sum = %.10f\n", sum); /* 1 / n * n だから 1 に戻るはず... */
if (sum == 1.0) {
printf("Sum == 1\n");
} else {
printf("Sum != 1\n");
}
return 0;
}
2.2
桁落ち
• 絶対値がほぼ等しい数値同士の減算を行った際，有効数字が減少する．
(例)1.2345678 × 109 (小数点以下 7 桁) −1.2345677 × 109 (7 桁) = 1.0 × 102
左辺の相対誤差 ∼ 10−8 , 右辺の相対誤差 ∼ 10−1 ．
• (大きい数) - (大きい数) = (小さい数) の引き算には注意を要する．
できれば (式を変形するなどして) 避ける．
• 例：2 次方程式の解の公式
2 つの解を x1 , x2 とすると
ax2 + bx + c = 0
√
b2 − 4ac
x1 =
2a
√
−b − b2 − 4ac
x2 =
2a
−b +
√
|b| ≫ |ac| の場合，|b| ≃ b2 − 4ac となるので，x1 , x2 のどちらかの分子の絶対値が極端に小さくなる． ⇒
桁落ちの発生
サンプルプログラム (桁落ちのケアをしていないもの)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(void)
{
char buf[64];
float a, b, c;
float x1, x2, p;
printf("ax^2+bx+c=0: (a, b, c) = ");
fgets(buf, sizeof(buf), stdin);
if (sscanf(buf, "%f,%f,%f", &a, &b, &c) != 3) {
fputs("Input error\n", stderr);
exit(-1);
}
2
if(a==0){ /* 2 次方程式でない場合は終了*/
fputs("a=0",stderr);
exit(-2);
}
p = sqrt(b * b - 4 * a * c);
x1 = (-b + p) / (2 * a); /* b と p が近いと桁落ちが発生 */
x2 = (-b - p) / (2 * a);
printf("x1 = %e, x2 = %e\n", x1, x2);
return 0;
}
桁落ちを防ぐ方法
1. x1 , x2 のどちらかは，桁落ちを起こさない．b > 0 の時は x2 , b < 0 の時は x1 が桁落ちを生じない．
2. 解と係数の関係:
x1 x2 = c/a を用いることで，もう一方の解を桁落ちをさせず求めることができる．
⇒ 実習で
√
• その他の例: 1 + x − 1 を求める場合．x ≪ 1 では桁落ちが発生．(コンパイラーが適切に処理してくれる
場合もある．)
√
x
桁落ちを防ぐために分子の有理化: 1 + x − 1 = √1+x+1
2.3
情報落ち
絶対値の大きく異なる 2 つの数を加えると，小さいほうの数値の持つ情報が失われる．
有効桁数が 6 桁の場合：123456+12.3456 = 123468
S=
1000
∑
k=1
1
k
1 1 1
1
1
+ + + ··· +
+
(1)
1 2 3
999 1000
1
1
1 1 1
S=
(2)
+
+ ··· + + +
1000 999
3 2 1
式 (1) は大きい値から加え，(2) は小さい値から加えている．式 (1) の場合，最後の方は大きな値+小さい値とな
るので，最後の方の小さい数字の誤差が大きくなる．(2) の場合，最初の方は小さい値同士を足して最後に大きな
値を足すので，誤差が少なくなる．
S=
様々な大きさの足し算を行う場合，できる限り小さい値から加える．
3
方程式の解法 –ニュートン・ラフソン (Newton Raphson) 法–
f (x) = 0 の解を求めたい．
解の近似値 x = x0 が分かっている場合
(
)
(
)
( )
1 d2 f
1 dn f
df
2
(x − x0 ) +
(x − x0 ) + · · · +
(x − x0 )n + · · · = 0
f (x) = f (x0 ) +
dx x=x0
2 dx2 x=x0
n! dxn x=x0
1 次の項まで取って変形すると
x = x0 −
f (x0 )
(f ′ (x0 ) ̸= 0)
f ′ (x0 )
この式に基づいて逐次的に真の解に近付けていく．
xn = xn−1 −
3
f (xn−1 )
f ′ (xn−1 )
f(x0)
f(x1)
f(x2)
x2
x1 x0
Figure 1: ニュートン・ラフソン法による解の求め方
• 出発点として、解に近い初期値を与えることができれば素早く収束する．
• 関数 f (x) = 0 が変曲点を持つ場合，その周辺に初期値を設定すると収束しない時がある．
⇒ 繰り返しの上限回数を設定する必要がある．
• 重根の場合，収束が遅くなる: 解に近づくと微分値が 0 に近づくため
• 収束の条件
ある小さな正の数 ϵ に対して
|xn+1 − xn | < ϵ, もしくは，
|xn+1 −xn |
< ϵ となった時計算を止める．後者は，解が 0 に近い時収束しにくくなる．場合に応じて使い分
|xn−1 |
ける．
• 例: x2 − a = 0 の解を求める．
#include <stdio.h>
#include <stdlib.h>
#define ACC 1.0e-12
int main(void){
double a,x,xt;
char buf[30];
printf("input a positive number : ");
if (fgets(buf,sizeof(buf),stdin) == NULL ||
sscanf(buf,"%lf",&a) != 1 || a < 0 ) { /* one positive double */
fputs("input error\n", stderr);
exit(-1);
}
x=a;
do {
xt=x;
x=0.5 * (x + a/x); /* x2 = x1 - f(x1)/f’(x1) */
printf("x=%25.16f xt-x=%25.16f\n", x, xt-x);
} while(x-xt >= ACC || xt-x >= ACC); /* |xt-x|<=ACC */
printf("sqrt(%f)=%25.16f\n",a, x);
return 0;
}
4
© Copyright 2024 ExpyDoc
