数値計算における誤差について【数値計算】

受け売りの言葉ですが，物理の問題は紙とペンで解くには複雑すぎることがしばしばです．

そんなときはコンピュータに任せます．

ただしコンピュータは有限桁の計算しかできません．

抽象的な理論としての数学だったら無限は扱えます．ただそれをコンピュータでやろうとすると，有限桁をもって無限を扱う必要があります．要はどうしても近似になってしまいます．

するとコンピュータによる計算，つまり数値計算で求めた結果は誤差を含みます．

今回は誤差の評価法として基本的なテイラー展開を用いた方法について書かせていただきます．

１．誤差を微小量としてテイラー展開

無限を有限で近似的に扱うにあたって生じる誤差を簡単に確認してみたいと思います．

例えば\(\sqrt{2}\)というのは無理数ですのでもしコンピュータ上で扱おうと思うと，有限桁で近似しなければならないです．そこでちょこっとだけpythonを使ってみます，

import numpy as np

x = np.sqrt(2)

print(f"x      = {x}")
print(f"x ** 2 = {x ** 2}")

とすると出力は

x      = 1.4142135623730951
x ** 2 = 2.0000000000000004

となります．\(\sqrt{2}^2 = 2\)ですから出力はぴったり\(2\)であるはずが，少しだけずれることになります．これが有限桁しか扱えないがために出てくる誤差です．

ここでこの誤差というのを少しだけ一般的に扱ってみます．

数学での厳密な値を\(x\)とし，それをコンピュータで有限桁で近似した値を\(\tilde{x}\)と表しましょう．このとき誤差\(h\)は\[h = \tilde{x} -x \]と表すことにしましょう．

例えば先ほどの例では

\begin{align*}
x &= \sqrt{2} \\
\tilde{x} &= 1.4142135623730951 \\
\end{align*}

となっていたわけです．コンピュータ内では\(x^2\)の代わりに\[\tilde{x}^2 = (x + h)^2 = x ^2 + 2 x h + h^2\]の計算をしたことになっていると見なせます．

よってこの「2乗する」という計算において誤差\(h\)の影響は\[\tilde{x}^2 -x^2 = \underline{2 x h + h^2}\]となるわけです．これは誤差\(h\)が2乗という操作に対して及ぼす誤差だと考えることができます．

ここでさきほどの具体的な値を思い出すと，\(h\)はおよそ\(10^{-15}\)ほどの非常に小さな値です．よって\(2xh\)と\(h^2\)の桁だけに注目しておおざっぱに比べると,

\begin{align*}
2xh &\simeq 10^{-15} \\
h^2 &\simeq 10^{-30}
\end{align*}

となりまして，\(2xh\)に比べて\(h^2\)ははるかに小さいです．

よって\[\tilde{x}^2 -x^2 \simeq 2 x h\]と見なすことができますから，数学的に厳密な値とコンピュータで計算した値との差は\(h\)の\(1\)乗で影響してくるということが分かります．

さてではもう一段階一般的にします．今まで「2乗する」という操作による誤差を考えましたが，一般の関数\(f\)で表されるとき誤差\(h\)による操作\(f\)に及ぼす誤差はどのように表現されるでしょうか．これは\[f(\tilde{x}) -f(x) = f(x + h) – f(x)\]と表現されることになります．

ここで関数\(f\)がテイラー展開できるとしましょう．すると\[ f(x + h) -f(x) = \frac{df(x)}{dx} h + \frac{1}{2} \frac{d^2 f(x)}{dx^2} h^2 + \cdots \]となりますから今\(h\)が十分に小さいのであれば，さきほどと同じように\(h\)の\(1\)次が最も影響が大きく，

\[f(x + h) -f(x) \simeq \frac{df(x)}{dx} h \]

として\(h\)の\(2\)次以上の項はすべて無視することができます．

ここまでの議論で大切なことは何かというと，コンピュータの有限性による誤差\(h\)を持った量に対して操作\(f\)を何の工夫もせずに行うと，\(h\)の\(1\)次の大きさの誤差が残るということです．

つまり数値計算の誤差を抑えるためには工夫が必要ってことですね．そりゃそうですよね．

２．例として数値微分

ここでは誤差\(h\)が操作に及ぼす影響を小さくする工夫について書きたいと思います．

ここでは操作\(f\)として微分をとってみましょう．

微分の定義式は\[\frac{df(x)}{dx} = \lim_{h \to 0} \frac{f(x + h) -f(x)}{h}\]ですね．これをコンピュータで数値計算しようということを考えます．

２－１．素直な方法

問題となるのは\(h\)ですね．有限の値しか扱えませんから，可能な限り小さな\(h\)をとって\[\frac{df(x)}{dx} \simeq \frac{f(x + h) -f(x)}{h}\]と素直に定義式そのままの形で近似することにします．

ここで，さきほどと類似の記号をとりましょう．すなわち数学的に厳密な微分係数を\[k = \frac{df(x)}{dx}\]とし，コンピュータで計算した値を\[\tilde{k}_1 = \frac{f(x + h) -f(x)}{h}\]とします．

さて繰り返しになりますが，コンピュータ上での\(h\)は\(0\)への極限がとれるわけではなく，ごく小さい有限の値です．ですから今回の数値微分においてはコンピュータ上の有限な\(h\)と数学的に厳密な\(h \to 0\)との差が微分にどの程度影響を与えるかというのを評価することになります．

ということで\(\tilde{k}_1 -k\)を計算してみます．テイラー展開すると

\begin{align*}
\tilde{k}_1 -k &= \frac{1}{h} (\underset{f(x+h)のテイラー展開}{\underline{f(x) + \frac{df(x)}{dx} h + \frac{1}{2} \frac{d^2 f(x)}{dx^2} h^2 + \cdots}} -f(x)) -\frac{df(x)}{dx} \\
&= \frac{d^2 f(x)}{dx^2} h + \cdots
\end{align*}

ということで確かに\(h\)の\(1\)乗の誤差が残っていることが分かります．

２－２．ちょっと工夫した方法

対称性を利用すると\(h\)の\(1\)乗の誤差をキャンセルすることができます．それを実際に書いてみます．つまり中点をつかって誤差\(h\)が均等に登場するようにすれば打ち消しあってくれるんじゃないかという目論見です．

\(f(x + h)\)と\(f(x -h)\)のテイラー展開をならべてみますと，

\begin{align*}
f(x + h) &= f(x) + \frac{df(x)}{dx} h + \frac{1}{2} \frac{d^2 f(x)}{dx^2} h^2 + \frac{1}{6} \frac{d^3 f(x)}{dx^3} h^3 + \cdots \\
f(x -h) &= f(x) -\frac{df(x)}{dx} h + \frac{1}{2} \frac{d^2 f(x)}{dx^2} h^2 -\frac{1}{6} \frac{d^3 f(x)}{dx^3} h^3 + \cdots \\
\end{align*}

引き算してみますと，

\begin{align*}
f(x + h) -f(x-h) &= \frac{df(x)}{dx} 2h + \frac{1}{3} \frac{d^3 f(x)}{dx^3} h^3 + \cdots
\end{align*}

となって\(h\)の\(2\)次の項が消えました．こうすることで微分の近似式として，\[\tilde{k}_2 = \frac{f(x+h)-f(x-h)}{2h}\]を採用すると誤差\(h\)の影響が\(2\)次に小さくなっていることがわかります．実際に書いてみると，

\begin{align*}
\tilde{k}_2 -k &= \frac{f(x+h)-f(x-h)}{2h} -\frac{df(x)}{dx} \\
&= \underset{f(x+h)-f(x-h)のテイラー展開}{\underline{\frac{df(x)}{dx} + \frac{1}{6} \frac{d^3 f(x)}{dx^3} h^2 + \cdots}} -\frac{df(x)}{dx} \\
&= \frac{1}{6} \frac{d^3 f(x)}{dx^3} h^2 + \cdots
\end{align*}

ということで誤差\(h\)の影響が\(2\)次に落ちていることが確認できました．

このように工夫をすることで誤差の影響を小さく抑えることができたりできなかったりします．