第1章　数値計算

1.1 数値計算と誤差

1.1.1 誤差・絶対誤差・相対誤差

数値誤差を生じることのない計算を理論計算と呼ぶことにしよう。理論計算とは、いわば無限桁計算である。数値計算においては、浮動小数点10進6桁（単精度計算）、あるいは10進15桁（倍精度計算）程度の有限桁数の数値を用いて計算するので、誤差は避けがたい。誤差が小さいうちは、実用上は困難が生じること少ないが、大きくなると数値計算の結果は信頼できない。数値計算の学習を始めるにあたって、まず誤差についての考察から始める。

ある数の真の値を $a$ 、その近似値を $x$ とするとき

\[e = x - a \tag{1.1}\]

を誤差（error）という。誤差の絶対値

\[\left\vert e \right\vert = \left\vert x - a \right\vert \tag{1.2}\]

を絶対誤差（absolute error）という。絶対誤差がある（小さな）正の数 $\varepsilon$ より小さいとき、すなわち、

\[\left\vert x - a \right\vert < \varepsilon \tag{1.3}\]

ならば、真の値 $a$ は近似値 $x$ の近傍の（小さな）範囲

\[x - \varepsilon < a < x + \varepsilon \tag{1.4}\]

の中にあることがわかるので、この $\varepsilon$ を誤差の限界という。誤差 $e$ が真の値 $a$ と何桁あっているかを言い表すために、誤差と真の値との比

\[e_R = \frac{e}{a} = \frac{x - a}{a} \tag{1.5}\]

を相対誤差（relative error）という。 $e_R$ が $10^{-5}$ なら、 $x$ は 10進5桁の精度（accuracy）をもつ。真の値が0あるいは0に近いときは、誤差の程度を相対誤差で言い表すことはできないから絶対誤差で言い表す。

例題 1.1

円周率 $\pi = 3.14159265358979···$ を小数点以下6桁で4捨5入して $3.14159$ としたときの誤差、絶対誤差、相対誤差、4捨5入の誤差の限界はどれだけか。

［解］

誤差 = 3.14159 - 3.14159265358979… = -0.00000265358979…
絶対誤差 = 0.00000265358979…
相対誤差 = -0.0000008446638…
誤差の限界 = 0.000005

1.1.2 許容誤差

絶対誤差がある (小さな) 値 $\varepsilon_A$ より小さければ、近似値 $x$ を真の値 $a$ の代わりに用いてよいと判断できるとき、この $\varepsilon_A$ を許容絶対誤差 (tollerable absolute error)という。また相対誤差の絶対値 $ \left\vert e_R \right\vert $ がある (小さな) 正の値 $\varepsilon_R$ より小さければ、近似値 $x$ が真の値 $a$ の代わりに用いてよい判断できるとき、この $\varepsilon_R$ を許容相対誤差(tollerable relative error) という。これらの許容誤差は普通はあらかじめわかっている値である。こうして、近似値 $x$ が真の値 $a$ に近いと判断できる基準は

\[\left\vert x - a \right\vert < {\varepsilon}_{A} \quad \text{許容絶対誤差} \tag{1.6}\] \[\left\vert x - a \right\vert < {\varepsilon}_{R} \vert a \vert \quad \text{許容相対誤差} \tag{1.7}\]

と書けるが、 $\left\vert a \right\vert$ が0、もしくは0に近いときも遠いときも使える誤差の基準としては

\[\left\vert x - a \right\vert < \varepsilon_A + \varepsilon_R \left\vert a \right\vert \tag{1.8}\]

が実用的である。一般に $a$ はわかっていないから、 $\vert a \vert$ の代わりに $\vert x \vert$ を使って

\[\left\vert x - a \right\vert < \varepsilon_A + \varepsilon_R \left\vert x \right\vert \tag{1.9}\]

とすれば、右辺の値は計算できる。

(1.9)のよく用いられる応用例を示す。いま $a$ に収束する数列の中の連続する2つの数 $x_n$ と $x_{n+1}$ が

\[\begin{aligned} \left\vert x_n - x_{n+1} \right\vert &= \left\vert(x_n - a) - (x_{n+1} - a) \right\vert \cr &\leqq \left\vert x_n - a \right\vert + \left\vert x_{n+1} - a \right\vert \cr & < \varepsilon_A + \varepsilon_R \left( \left\vert x_n \right\vert + \left\vert x_{n+1} \right\vert \right) \end{aligned} \tag{1.10}\]

を満たすならば、$x_n$ も $x_{n+1}$ も、ともに $a$ のよい近似値であると見なす。

許容誤差を小さくすればするほど高い精度の数値を要求することになる。最高精度を要求した場合、許容絶対誤差 $\varepsilon_A$ は、(1.6) より

\[0 \leqq \vert x - a \vert < \varepsilon_A \quad \text{すなわち} \quad \varepsilon_A > 0\]

であるから、用いているコンピュ－タの記憶できる最小の正の数値である。この数値より小さい数値 $\varepsilon < \varepsilon_A$ を記憶させようとすると $\varepsilon$ = 0 とアンダーフローする。また、許容相対誤差 $\varepsilon_R$ の最小値 $\varepsilon_M$ は、(1.7) より

\[\left\vert a \right\vert \leqq \left\vert x - a \right\vert + \left\vert a \right\vert < \left( 1 + \varepsilon_M \right) \left\vert a \right\vert \quad \text{すなわち} \quad 1 + \varepsilon_M > 1\]

を満たす最小の正の数値である。この値はマシン・エプシロンと呼ばれる。 $\varepsilon_M$ より小さな数 $\varepsilon < \varepsilon_M$ は、1に加算されるときには数値的には無視され $1 + \varepsilon = 1$ となる。

倍精度計算の場合、指数部7ビットで仮数部16進14桁のIBM方式では $\varepsilon_A = 16^{-65} \simeq 5.3976\cdot10^{-79}$ 、 $\varepsilon_M = 16^{-(14-1)} \simeq 2.2204\cdot10^{-16}$ である。また、指数部11ビットで仮数部52ビットのIEEE方式では $\varepsilon_A = 2^{-1074} \simeq 4.9407\cdot10^{-324}$ 、 $\varepsilon_M = 2^{-52} \simeq 2.2204\cdot10^{-16}$ である ^$1)$ (問題 1-1 参照)。

$^{1)}$ 「付録 A 浮動小数点数の表現」

1.1.3 丸めの誤差

コンピュ－タが記憶する数値は、有限桁の2進数であるから、数値を記憶する際に有限桁の2進数に納まるように丸められる。このとき生じる誤差を、丸めの誤差(round off error) という。われわれが見る数値は10進数であるが、10進数を2進数に変換するときにも誤差が生じる。たとえば、10進数の0.1は、2進数では

\[0.0~0011~0011~0011~0011~0011......\]

と無限循環小数であり、コンピュ－タはこれを有限のビット数で丸めて記憶する。また、乗算の積は、いつも乗数と被乗数の桁数の合計の桁数が必要である。除算の商はほとんどの場合被除数の桁数より大きい。加減算の場合は、小数点の位置合わせが行われて、絶対値の大きな数の方が優先し、小さな数の末尾は捨てられてしまう。これらの演算の結果は一定の桁数に丸められ、丸めの誤差が生じる。

丸めの誤差自体は末尾の小さな値の誤差であるが、これが集積していき、計算法によっては増幅されて行くと、大きな誤差になり得る。

丸めの誤差を完全に避けることは困難である。丸めの誤差を小さくおさえるための一般的処方は、計算の精度を大きくするよりなかろう。原則として、数値計算では単精度ではなく、倍精度の計算をするように心がけることが必要である。

1.1.4 桁落ち

丸めの誤差は数値の末尾の小さな値が失われることによって生じる誤差である。ここで述べる桁落ちは、大きな誤差を生み出す現象で、数値計算でもっとも警戒しなければならない現象である。

桁落ち (loss of significant digits) は絶対値のほとんど等しい2つの数値の間の減算 (2つの数値が同符号の時) と加算 (異符号の時) のときに起こる。たとえば、有効数字6桁の2つの数

\[a = 1.23456 \quad \text{と} \quad b = 1.23455\]

との差 $a - b$ は、有効数字 $1$ 桁の $0.00001$ となり、有効数字は $6 - 1 = 5$ 桁失われる。すなわち $5$ 桁分桁落ちしてしまっている。しかも上位の桁が失われている。桁落ち以降の演算 (特に乗除算) の精度は落ちてしまい、結果は信用できなくなる。

桁落ちは是非とも避けなければならない。桁落ちを避けるには、桁落ちが起こる可能性のある演算は行わないようにする。

［例 1.1］

2次方程式 $ax^2 + 2bx + c = 0$ の解は、公式より
\[x_1 = \frac{-b + \sqrt{b^2 - ac}}{a} ,~x_2 = \frac{-b - \sqrt{b^2 - ac}}{a} \tag{1.11}\]
であるが、$b^2 \gg \left\vert ac \right\vert$ なら、$b > 0$ のとき $x_1$ の分子が、$b < 0$ のとき $x_2$ の分子が桁落ちする。したがって、 $b > 0$ のときは $x_2$ を求めて、$x_1$ は根と係数の関係

$\displaystyle x_1 x_2 = \frac{c}{a} $ より $\displaystyle x_1 = \frac{c}{ax_2}$

と求める。 $b < 0$ のときは $x_1$ を公式より求め、$x_2$ は
\[x_2 = \frac{c}{ax_1}\]
と求める。まとめて書けば、2次方程式の解は
\[x_1 = \frac{-b - \rm{sign}(b) \sqrt{b^2 - ac}}{a} ,~ x_2 = \frac{c}{ax_1} \tag{1.12}\]
とする。ここに、$\rm{sign}(b)$ は $b \geqq 0$ のとき $+1$、$b < 0$ のとき $-1$ である符号関数である。

1.1.5 打切り誤差

もともと無限大や無限小の極限で定義されている値を、有限のところで打ち切った値で計算するとき生ずる誤差を、打切り誤差 (truncation error) という。

［例 1.2］

導関数の値は
\[\frac{df}{dx} = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h} \tag{1.13}\]
で定義されているが、これを $h$ を有限のままにした値で近似したときの打切り誤差は
\[\frac{f(x + h) - f(x)}{h} - \frac{df}{dx} \tag{1.14}\]
である。数値計算では $h$ の無限小は不可能であるし、無限小でなくても $h$ が小さいと、関数値の差の桁落ちは避けられない。

同じようなことが、無限小の量を無限大個加える積分の場合にもいえる。

［例 1.3］

初等関数の値は、無限級数で定義される。 $\sin x $は
\[\sin x = \sum^{\infty}_{k=0} \frac{(-1)^k}{(2k + 1)!} x^{2k+1} \tag{1.15}\]
これを $n$ 項で打ち切って求めれば、やはり打切り誤差
\[\sum^{n}_{k=0} \frac{(-1)^k}{(2k + 1)!} x^{2k+1} - \sin x = - \sum^{\infty}_{k=n+1} \frac{(-1)^k}{(2k + 1)!} x^{2k+1} \tag{1.16}\]
が生じる。

1.2 数値計算の手順とPAD

本来は「PADについて」と書内でFORTRANコードを用いることが書かれている。表示が困難であるためこのリポジトリではPAD図の記載を省略する。また、FORTRANによるプログラムはJuliaによって置き換える。

1.2.1 構造化プログラミング

1.2.2 PADの制御構造

1.2.3 PADとFORTRAN

1.3 第1章の問題

問題 1-1

許容絶対誤差の最小値 $\varepsilon_A$ は数値的に $\varepsilon_A > 0$ を満足する最小正値、マシン・エプシロン $\varepsilon_M$ は、数値的に $1 + \varepsilon_M > 1$ を満足する最小正値である。 $\varepsilon_A$ および $\varepsilon_M$ を求める Julia プログラムを書け。

［解］

倍精度計算における $\varepsilon_A$(epsA) および $\varepsilon_M$ (epsM) を求めるプログラムの例は、それぞれ次の通りである。

function calc_εA() 
    εA = Float64(1.0)
    εX = Float64(0.5)
    while εX > 0
        εA = εX
        εX /= 2
    end
    return εA
end

function calc_εM()
    εM = Float64(1.0)
    while 1 + εM/2 > 1
        εM /= 2
    end
    return εM
end

問題 1-2

図 1.1 の簡単な PAD を Julia の文に翻訳せよ。

［解］

スカラー積

function scalar_product(X, Y)
    S = 0
    for (x, y) = zip(X, Y)
        S += x*y
    end
    return S
end

行列とベクトルの積

function product_mat_vec(A, X)
    N, M = size(A)
    b = zeros(N)
    for i = 1:N
        for j = 1:M
            b[i] += A[i, j] * X[j]
        end
    end
    return b
end

2つの行列の積

function matrix_multiplicate(A, B)
    N, K = size(A)
    K, M = size(B)
    C = zeros(N, M)
    for i = 1:N
        for j = 1:M
            for k = 1:K
                C[i, j] += A[i, k] * B[k, j]
            end
        end
    end
    return C
end

2次方程式の解

function completing_square(a, b, c)
    if a == 0
        if b == 0
            if c == 0
                println("xは任意")
            else
                println("xは不能")
            end
            return
        else
            x = - c / b
        end
    else
        d = b^2 - (4 * a * c)
        if d >= 0
            x1 = (-b - sqrt(d)) / 2a
        else
            x1 = (-b - sqrt(Complex(d))) / 2a
        end
        x2 = c / (a * x1)
        x = (x1, x2)
    end
    return x
end

第1章 数値計算