計算機言語 II 第 12 回 誤差

計算機言語

II

12

http://www.math.u-ryukyu.ac.jp/~suga/gengo/2018-2/12.pdf

1

今回は実数double 型のデータに現れる計算誤差の話をします. まず始めに, 計算機イプシロンを求めます. 計算機イプシロンとは

1.0 +ε >1.0

となる最小の正の数 εのことです. もちろん, 実数の場合は, (有理数でも)上の式をみたす εはいくらでも0 に近くとれるので意味がありません. (数学では最小値は存在しないというのでした.)

計算機の場合, double 型でも変数はトビトビの値を取りますから,上記のεが定まります. 次のプログラム

を12-1.cという名前で入力し,コンパイルして実行してみて下さい. 処理内容については,各自考えて下さい.

/* File name 12-1.c */

/* Machine epsilon */

#include <stdio.h>

main() {

int i; /* 繰り返しの数を数える */

double e; /* 計算機イプシロン */

e = 1.0; /* 実数型では, 代入する数も小数にします */

/* 2で割っていき, 1.0 + e = 1.0 となったらループを抜けます */

for (i = 0; 1.0 + e > 1.0; i++) { e /= 2.0;

}

/* ループを抜ける一つ手前の値が計算機イプシロンです */

printf("machine epsilon = 2^(%d)\n", -(i-1));

}

2

実際上の問題で,コンピュータなどを用いて数値計算を行う場合には, 近似値を用います. そのため,必然的 に誤差の問題が生じます.

2.1

真の値Aの近似値aに対して,近似値と真の値の差の絶対値|a−A|を近似値aの誤差といいます. また, ある小さい正の数εに対して

|a−A|≤ε

が成り立つとき,εを近似値aの誤差の限界といいます. この時,真の値の範囲は, 次の様になります. a−ε≤A≤a+ε

例えば,√

2 = 1.41421356· · ·の小数第6位を四捨五入した値a= 1.41421を√

2の近似値にとると a−0.000005≤√

2< a+ 0.000005 よって,aは不等式|a−√

2|≤0.000005を満たし,この近似値の誤差の限界は0.000005,すなわち5×10⁻⁶ となります.

2つの地点P, Q間の距離の測定値は 100m,別の2地点R, S 間の距離の測定値は 42195mで,これらの誤 差の限界はともに10cmであったとします. このとき,誤差の限界は同じでも,P, Q 間の距離と,R, S 間の距 離の大きさが違いすぎるため,測定値の近似の程度はかなり異なると考えます.

このようなとき, 真の値Aの近似値aに対して,近似の程度をみるには,次の相対誤差を考えます. a−A

A

2.2

コンピュータや電卓では, 入力された10進数は内部で2進数に変換されて, 計算などの処理が行われてま す. 2進数でも, 10進数の場合と同様に小数が考えられ, 2進数の101.011とは,

1×2²+ 0×2¹+ 1×2⁰+ 0×2⁻¹+ 1×2⁻²+ 1×2⁻³

と表される数のことであり, 10進数で表すと 5.375です. 10進数の有限小数は, 2進数に変換すると無限小数 になる場合があります. 例えば, 10進数の小数0.1を2進数に変換すると

0.00011001100110011001100· · · ·

と,無限小数になります. コンピュータや電卓には, 0と1を記憶する場所がたくさんありますが, 無限ではな

コンピュータは,非常に多い回数の計算を高速で行いますが,無限回の計算をするわけではありません. 弧度 法で表された角xの三角関数sinxは

sinx=x−x³ 3! +x⁵

5! −x⁷

7! +· · · ·

と,項が無限に続く式になります. この式を用いてコンピュータでsinxの値を求める場合,無限回の計算は行 えないので,右辺を有限個の項で打ち切って,その近似値を求めることになります. このような事情のために生 じる誤差を,打ち切り誤差といいます.

コンピュータを利用するときに生じる上のような誤差のほかに, 有効数字の極端な桁落ちによる誤差もあり ます. 例えば,有効数字が5桁の非常に近い2つの数

a= 2.1438, b= 2.1425

の差はa−b= 1.3×10⁻³ となり,有効数字が2桁に落ちます. そのため,この差にbを掛けた (a−b)×b= 1.3×10⁻³×2.1425 = 2.78525×10⁻³

では,有効数字は最初の2つだけになり,結果は2.8×10⁻³とすべきです. 即ち,相対誤差が,拡大された事に なります.

以上のような誤差は, 1回の計算では微小であっても, 例えば, コンピュータが得意とする繰り返し計算など では,計算回数が多くなるに従って誤差が累積し,無視できなくなる可能性が多くあります.

3

2次方程式x²−50000.0001x+ 5 = 0を考えます. 因数分解をすると, (x−50000)(x−0.0001) = 0なの で,この方程式の根は, 50000と0.0001です. これを根の公式を素朴に利用して計算しますと, 0.0001の方が 正確に計算できません. より正確な計算をさせるには,絶対値の大きい根はそのまま根の公式で計算し,小さい 方は,根と係数を利用した割り算で計算させます.

数学的には,簡単な計算もコンピュータでやらせるには注意が必要であることの,最も分かりやすい例です.

/* File name 12-2.c */

#include <stdio.h>

#include <math.h>

int main() {

double alpha, beta;

/* 素朴な根の公式による計算 */

alpha = (50000.0001 + sqrt(50000.0001*50000.0001 - 4*0.0001))/2.0;

beta = (50000.0001 - sqrt(50000.0001*50000.0001 - 4*0.0001))/2.0;

printf("alpha=%15.10f, beta=%15.10f\n", alpha, beta);

beta = 5.0/alpha; /* 根と係数の関係で計算 */

printf("alpha=%15.10f, beta=%15.10f\n", alpha, beta);

}

4

∑N n=1

1

n^{を考えます}. N → ∞^{でこれは発散しますが},有限和は当然値を持ちます. これを, 1からNまで和 をとった時とN から1 まで和をとった時では, Cを用いて計算すると, 上の桁落ちの効果で違った値になり ます. より正確な値は,大きい方から和を計算したsum2です.

/* FIle name 12-3.c */

#include <stdio.h>

#define MAX 10000000

int main() {

double sum1=0.0;

double sum2=0.0;

int i;

for (i=1; i <= MAX; i++){

sum1 += 1.0/i;

}

for (i=MAX; i >0; i--){

sum2 += 1.0/i;

}

printf("sum1=%17.16f\nsum2=%17.16f\n", sum1, sum2);

}

レポート問題

:

1

10

(

)

1. 2次方程式x²+ax+bの数値解を求めるプログラムを書け. 件名: Quadratic equation

• a,bを入力させて,数値解を出力する.

• ^{虚数解を持つ時には}, そのことを表示する.

• 計算は上で述べたように, 桁落ちしないようなプログラムにする. (絶対値の大きい解を先に計算 し,小さい解は割り算で計算する.)

2.

1000000∑

n=1

(−1)ⁿ⁻¹

n = 1−1 2+1

3 −1

4 +· · · の値を大きい方から和をとった時と, 小さい方から和をとっ た時で比較せよ. どちらがlog 2により近いか? log 2の値は, /usr/include/math.hにM_LN2で定義さ れている. 件名: log 2

ちなみに, log 2を上の式で計算したりはしない. log 2 =−log1

2 ^{を用いれば},効率よく計算できる. テキストの置き場所: ftp://ftp.math.u-ryukyu.ac.jp/pub/gengo/2018-2/