pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Pythonで学ぶ数値計算のツボ(桁落ち)

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使って数値計算の桁落ちをみてみたいと思います。動作の確認はJupyterNotebookを用いた場合で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。



桁落ちの具体例

 \ \displaystyle{\sqrt{X+1}-\sqrt{X}}

このような計算で、 X が非常に大きな数の場合(ここでは  X=10^{15} としましょう)計算精度に誤差が生じる場合があります。というのは、 \sqrt{X+1} \sqrt{X} とがほとんど同じ値になってしまい

import numpy as np
s = np.sqrt((10**15)+1)
t = np.sqrt(10**15)
print('sqrt(X+1) =', s, ': sqrt(X) =', t)

から

 \ \displaystyle{\begin{aligned} \sqrt{X+1}=31622776.60168381\\ \sqrt{X}=31622776.60168379 \end{aligned}}

その数値的な違いがすでに有効数字と同じオーダーになってしまっていて、計算結果の精度が不明瞭になってしまうためです:

s-t

より

 \ \displaystyle{\sqrt{X+1}-\sqrt{X}\overset{?}{=}1.862645149230957\times 10^{-8}}

となります。そこで、大きい数同士の差をとることに誤差が生じる危険性があるので、分子を有理化することによって、差をとるという操作をなくしてしまうことができます:

 \ \displaystyle{\sqrt{X+1}-\sqrt{X}=\frac{\left(\sqrt{X+1}-\sqrt{X}\right)\left(\sqrt{X+1}+\sqrt{X}\right)}{\sqrt{X+1}+\sqrt{X}}=\frac{1}{\sqrt{X+1}+\sqrt{X}}}

このようにして計算すると

1/(s+t)

より

 \ \displaystyle{\sqrt{X+1}-\sqrt{X}=1.5811388300841893\times 10^{-8}}

と計算されます。このようにすると、有効数字の桁数を保ったまま計算できます。



具体例:2次方程式

上のような具体例が現実的に起こりうるのは、2次方程式

 \ \displaystyle{ax^2+bx+c=0}

で、解の公式から

 \ \displaystyle{x=\frac{-b+\sqrt{b^2-4ac}}{2a} \qquad (b>0,\ a\neq0)}

 b がじゅうぶん大きな値だった場合です。

たとえば  a=1/2,\ c=-1/2 とすると

 \ \displaystyle{x=-b+\sqrt{b^2+1}}

で、まさに上の例で  X=b^2 とした場合になります。


2次方程式の数値解

ここで、Pythonの代数方程式の数値解を出力する組み込み関数を使うと(詳細は以下の記事を参照してください)

pianofisica.hatenablog.com

2次方程式

 \ \displaystyle{\frac{1}{2}x^2+\sqrt{10^{15}}x-\frac{1}{2}=0}

では

np.roots([1/2, sqrt(10**15), -1/2])

より、(2つのうちの1つの)解が

 \ \displaystyle{-b+\sqrt{b^2+1}=1.58113883\times 10^{-8}\qquad b=\sqrt{10^{15}}}

であると計算されます。これは上の例で分子を有理化した場合に得られた数値と一致しており、有理化せずに差をとった場合の結果が誤りだったことを示しています。

以上のことから、数値計算で大きい数同士の差をとるときには、有効数字の桁数にじゅうぶん注意する必要があることがわかります。



キーワードPython数値計算、桁落ち