Edited at

数値計算の常識 概要(個人的メモ)

More than 3 years have passed since last update.


第1章:数の表現と誤差


16進切り捨てと2進切り捨て

16進切り捨てのとき$ε=6*10^{-8}〜10^{-6}$

2進切り捨てのとき$ε=3*10^{-8}〜6*10^{-6}$


ループ制御時の誤差の考慮

$0.1≒(0.CCCCCD)_{16}*2^3$を10回足すと、10進で約1.0000001192になって1.0よりわずかに大きくなる。

この場合、whileの条件をwhile(x<=1.0)にしていては条件を満たさない。

こんなときには、ループの制御に整数変数を持たせるとか、whileの条件に刻み幅0.1より小さい、かなりの"余裕"をもたせて"while x<=1.001 do"のようにするとかが定石となっている。


第2章:桁落ちに気をつけよう(その1)


桁落ちとは

絶対値がごく近い2数を足したり引いたりして結果の絶対値が小さくなるような計算をすると、絶対値が小さくなった分だけ相対誤差が大きくなる。すなわち有効数字が減る。このような現象を桁落ちと呼んでいる。

結果が0に近くなるような計算をする場合、必ず桁落ちが生じている。


式の変形による桁落ちの回避

例えば


  • 分母の有理化

  • 三角関数の半角の公式の利用


m重根の有効数字

m重根の有効数字の桁数は、計算桁数の1/mになると言われている


加減算と乗除算の誤差

数値計算の誤差の観点からすると、乗除算よりは加減算のほうがはるかに難物


第3章:桁落ちに気をつけよう(その2)


収束判定と誤差

$f(x)=x^3-15.70875x^2+61.69729x-0.04844725$

の改を求めたいとする。

この解は、0.0007853982、7.844763、7.863201の付近にある。

ここで、$x=0.0007853982$として10進7桁程度の有効数字で計算しているとすると、

各項では以下のように誤差が生じる。

右辺第一項:$±7.5*10^{-16}$

右辺第二項:$±1.0*10^{-11}$

右辺第三項:$±2.5*10^{-8}$

したがって、$f(x)$には$±2.5*10^{-8}$のオーダーくらいの誤差が混入する可能性がある。このことは$|f(x)|$が$10^{-8}$のオーダーになったら、それ以上細かくxの値を調整してもそれに伴う$f(x)$の値の変化には意味がないことを示している。

$f(x)≒0$の"$≒0$"ということの意味は$f(x)$の計算方法とxの値とに強く依存しているということである。

このことに常に注意していないと適切な算法設計は出来ない。

収束判定について、上の話を適用すると、

本来、上記のような"誤差に関する慎重な配慮"をしながら収束判定や計算打ち切りの事を決めなければならないのに、それを"ε"という量にしわ寄せしているケースが多い。


第4章:たった一回の計算なんて


数値計算の正しさ

数値計算の結果が正しいこと、あるいは、どの程度の許容誤差の範囲内では正しいかということを、かなりの説得力をもって示すこと(計算をしている本人自身に対しても、他人に対しても)は、計算作業自体と同様に重要なこと。

個々の問題ごとに、その背景にある数値解析の理論を十分理解した上で、それぞれの問題に則した適切なやり方で、計算結果の吟味を行わなければなれない。

その際、かなりの一般性をもっていえることは、

『ただ一回計算をしただけでは結果の正しさについてはほとんど何もわからない』

『系統的に計算方法を変えて2回以上計算すると非常に有用な情報が得られる』

ということ。


  • 刻みを変えて計算することの効用

  • 計算桁数を変えて計算することの効用


第5章:逆行列よさようなら


  • 逆行列は追放しLU分解を使うべき

  • LU分解を経由しても連立方程式を解く際の計算回数は変わらない

  • 逆にAが実対称(あるいは複素エルミート)行列であれば、計算回数は半分になる(対称ホレスキー分解)

実用的な観点からは、Gauss-Jordan法や逆行列とLU分解との間にはもっと重要な差異がある。

大規模な問題ではAはたいていゼロ要素の多い疎行列であるが、そのようなAに対しても逆行列はほとんどの場合全ての要素が非ゼロである。

それにたいして、LおよびUはかなりAの疎構造を継承して、ゼロ要素が多く、しかも、数値計算に入る前に"構造的に0"である要素があらかじめわかる。


第6章:単位と次元

物理法則などで単位を持った式の数値計算をする場合(オイラーの公式etc)、パラメタが次元を持ってしまう。

つまり、単位が変わるとパラメタのオーダーも変更しなければならなくなってしまう。

そうならないように、数値計算に入る前に、無次元量の関係式にして、パラメタをなるべく少ししか含まない形に変形しておくべきである。


第7章:数値積分法-台形則を使いこなすには


基本的な性質

式で書いてある関数はよほどの事がない限り式で微分できるのに対して、積分の方は極限られた形のものしか式の形で積分することはできない。

そこで、定積分の数値計算方法である"数値積分法"が重要視されるわけである。

区間数Nを倍(刻み幅hを半分)にすると台形則の誤差は約1/4になる。この事実は本性のよくわかっていない関数f(x)を相手にしているとき、あるいは計算上の丸め誤差の大きさの見当がつけ難いときなどに、"台形則が台形則らしく"働いているかどうかを確認するのに役に立つ(実際には$I_N - I_{N/2}$がNを倍にした時に約1/4になっている事を確認する)

しかし、丸め誤差の観点からは、やたらにNを大きくすればよいというものではない事にくれぐれも注意されたい。


高精度公式の導出


  • Richardson補外

  • Simpson則


不連続点があればそこで分割せよ


端点の特異性は変数変換で除け


  • Aitken加速


周期関数に対する台形則

解析的な周期関数$f(x)$の全周期[a,b]にわたる積分を台形則で計算したときの誤差は、$f^{(m)}(a)=f^{(m)}(b)$(すべてのmに対して)であるから、「どんなpに対しても$I_N - I=O(h^{(2p+2)})$となる」ことになる。すなわち台形則はきわめて良い性能を持つことが期待される。


無限積分に対する台形則

"無限和"が実際には非常に僅かな項数の"有限和"ですんでいることにも注目すべきである。


万能変数変換

たいていの積分を、前節で示したような、台形則がきわめて早い収束性をもつ形の無限積分に変換してしまう"2重指数変換"とよばれる"万能型"の変数変換が知られている。


多重積分

数値積分をするにあたり、特異点や端点について事前に処理しておく必要があるが、これらは、次元の増加とともに耐え難いものとなる。

多重積分に対しては、結局モンテ・カルロ法等に頼らざるをえない場合が多い。 


第8章:数値微分法−打ち切り誤差と丸め誤差の闘い


数値微分の"刻み"

$f_1(x,h) = (f(x+h) - f(x))/h$

刻み幅hを小さくすると、初めのうちは誤差が減少するが、あるところから不規則な挙動をはじめ、そしてそれからはかえって誤差が大きくなる傾向がみられる。

初めのうちの傾向は

「$h$->0のとき$f_1(x,h)$->$f'(x)$」

という"理論"の反映であり、途中から先は

「$f(x+h)-f(x)$の計算の際の桁落ち」

の影響である。

丸め誤差と打ち切り誤差の両者の和が数値微分の誤差だから、hは大きすぎても小さすぎてもいけない。

ごく荒っぽくいって、$f_1(x,h)$を用いるときには

「hを演算桁数の半分くらいの桁のところに選ぶのがよく、そのとき数値微分の結果の有効数字は約半分になっている」

ということになる。

$f_2(x,h) = (f(x+h) - f(x-h))/2h$

を使うと、

「hを計算桁数の1/3くらいの桁のところに選ぶのがよく、そのとき結果の有効数字は計算桁数の約2/3になっている」


数値微分の加速


  • Romberg1段公式

  • Aitkenの公式

この考え方では「より高精度な公式と、より大きなhを使うのがよい」となりそうだが、もちろんそれには限度がある。hは展開式

$f(x+h) = f(x) + hf'(x) +$ ・・・

の収束が速いような大きさ("小ささ")でなければならない。

実践的には、 $f_1(x,h)$と$f_1(x,2h)$との"差"を観察して、

「hを半分にしたとき"差"が約半分になる傾向が保たれる範囲で最も小さなh」

を選ぶのがよい。$f_2(x,h)$に対しては、

「hを半分にしたとき1/4になる傾向が保たれる範囲」

とする。


第9章:Newton法


どうやって計算を終了するか

「$x=\alpha$の近くで関数値$f(x)$を計算する際に予想される丸め誤差などの計算誤差の大体の大きさ$\delta$の見当を何らかの方法で付けておいて、

 $|f(x^{(\nu)})|<\delta$

となったら計算を終える」


どうやって計算を開始するか

出発値$x^{(0)}$の定め方も難しい。下手なところから出発すると、収束どころか、とんでもない発散・振動を起こすことさえある。

最も基本的な"収束性向上"のための手法として"減速"がある。

それは、

$x^{(\nu+1)}=x^{(\nu)}-\frac{f(x^{(\nu)})}{f'(x^{(\nu)})}$

の代わりに

$x^{(\nu+1)}=x^{(\nu)}-\mu^{(\nu)}\frac{f(x^{(\nu)})}{f'(x^{(\nu)})}$

を用いるということである。


第10章:複素数の計算


複素数の表現

複素数$z$の表現としては、実部$x=Re$ $z$と虚数部$y=Im$ $z$とで表す"直角座標表示"$z=x+iy$と、

絶対値$r=|z|$と偏角$\theta=arg$ $z$とで表す"極形式"$z=rexp(i\theta)$

とを考える。

両表現の間には以下の関係がある。

$x=rcos\theta$, $y=rsin\theta$

$r=\sqrt{x^2+y^2}$, $\theta=arctan(\frac{y}{x})$

もし、rがどの程度の大きさか知りたいのであれば、$r=|z|$の代わりに

$max(|x|, |y|)$ あるいは $|x|+|y|$

で十分なことが多い。

rの精密な値が欲しいときには

$|x| (y=0のとき)$

$|y| (x=0のとき)$

$|x|・\sqrt{1+(y/x)^2} (|x|≧|y|>0のとき)$

$|y|・\sqrt{1+(x/y)^2} (|y|>|x|>0のとき)$

とする。

xやyを浮動小数点表示としたとき、その指数部の絶対値が許される最大値の半分以上になっていると、rそのものはまだ"指数部あふれ"を起こさないのに、$x^2$や$y^2$は指数部あふれを起こしてしまう。

したがって、上のような計算のほうがよい。

複素数$z,z_1,z_2$について$(z_1/|z|)・(z_2/|z|)$というような計算をすることがあるが、このようなときはもちろん$z_1・z_2/|z|^2$としたほうがよい。

また、極形式直行座標表示に戻すときは、できるだけ$n\pi/2$の形の角の部分を"式の変形"で処理してしまうのが良い。


第11章:代数方程式

ここでいう代数方程式は「(一般には複素数の係数をもつ)一変数の多項式を0に等しいと置いた方程式」という狭義にとる。

すなわち、"代数方程式を解く"というのは"多項式の零点を求める"ことであるとする。


誤差の見積もり

Smithの法則


第12章:常微分方程式の初期値問題

常微分方程式の初期値問題の数値解法というのは、数値積分法の親玉みたいなもの。


Euler法

刻み幅を変えて何回か計算しましょう。

高々倍の手間をかけることを厭わなことによって、単なる"計算のやりっぱなし"ではなく、結果に対する自信と保証とが得られる。微分方程式の数値解にときどき現れる"不安定現象"とか"幻影解"とかも、このようにな眺めればすぐ発見できる。


Euler法で満足するな

例えば小数点以下2桁目まで正しい答えを出す刻み幅hが0.05だったとする。

誤差がhに比例するというEuler法の性質からして、もう2桁正しく出すhは0.0005となる。

これでは複雑な問題には計算時間がかかって仕方ない。また丸め誤差の集積という伏兵も頭をもたげてくる。


種々の高精度公式

Runge-Kutta法、多段階法を知っておけば十分


第13章:数値計算の手間の理論と実際

単純に指数計算する場合も部品みわけて作って、組み合わせた方が良い


第14章:数値計算の中の非数値計算

スパース行列の効率的な計算をしましょう


第15章:補間(内挿)

性質の良い補間関数を使いなさい


第19章:数列の収束の速度

数値計算の立場からは、単に"収束する"かどうかが問題なのではなくて、"どのくらい速く"収束するかが問題なのである。


第20章:数列の収束の加速


  • Romberg-Richardson法

  • Aitken法

  • $\epsilon$算法


第23章:モンテ・カルロ法


  • 誤差は$1/\sqrt{N}$で減少

  • 誤差の目安は2sd

  • 次元nが増えると計算すべき点は次元数の指数関数的に増大