計算物理学:No.1

1 はじめに

この講義では、既にある程度のプログラミングについて慣れてきた人を対象に収束加速法、常微分方程式などを学習し、物理現象に使えるまで馴染んでいく。

2 数列の収束

3 年前期の数値計算法概論で扱った数値積分、ニュートン法で経験したが、一般に有限回の数値計算では近似値が得られる。勿論演算量を大きくすれば近似が良くなるが、少ない演算量から良い近似が推定できれば効率的であるし、桁落ち、丸め誤差などの悪影響も受け難くなる。

そこで問題を一般化し、数列 a_{1}, a_{2}, a_{3}, ⋯ が極限値 a に収束する

lim_{n→∞} a_{n} = a  (1)

として、その収束の早さの分類と、収束加速法 (数値積分で扱った Richardson 加速が 1 例) について議論する。

2.1 分類

2.1.1 1 次収束

|a_{n+1} − a| ≤ M|a_{n} − a| (0 < M < 1)  (2)

となるようなもの。これは n が大きいところで

a_{n+1} − a ≈ λ(a_{n} − a) (0 < |λ| < 1)  (3)

となっていることが多い。この時、 λ, b を定数として

a_{n} = a + bλ^{n} + ⋯  (4)

と表される。

収束の早さは λ の絶対値により、 |λ| が小さければ収束が早く、 |λ| が１に近いと収束が遅い。収束加速法の良い対象となる。

例としては重根近くでのニュートン法での振舞いなど。

2.1.2 2 次収束

|a_{n+1} − a| ≤ M|a_{n} − a|^{2} (0 < M < 1)  (5)

となるもの。

例としてはニュートン法がうまく働く時で、解の近くでの振舞い。この場合は、大体 n が 1 つ増すごとに正しい桁数が倍になるので、収束が十分早く、特に加速法の必要がない。

2.1.3 p 次収束

上記 2 つを一般化し

|a_{n+1} − a| ≤ M|a_{n} − a|^{p} (0 < M < 1)  (6)

となるもの。 p が整数でない場合もある。

2.1.4 O(n^{−α}) の収束

|a_{n} − a| ≤ cn^{−α} (c > 0, α > 0)  (7)

となるようなもの。しばしば、

a_{n} = a + b_{1}n^{−α} + ⋯  (8)

のように表される。

α が大きくなるにつれ、収束は早くなるが、十分早いとはいえない。 α が 0 に近付くと、収束が大変遅くなる。

例としては台形公式による数値積分で関数が滑らかな場合 (α = 2)。

この型のものは n が等比数列をなすような a_{n} の部分列を選ぶと 1 次収束になる。例えば、台形公式による数値積分で n = 2^{m}(m = 0, 1, 2, 3, ⋯) として、 a_{1}, a_{2}, a_{4}, a_{8}, ⋯ ととると、 (4) の型になり、 Richardson 補外が応用できる。

この型のもので、収束が交代的であるもの

a_{n} = a + (−1)^{n}(b_{1}n^{−α} + ⋯)  (9)

もある。 2.1.5 混合型

1 次収束と、 O(n^{−α}) の両者が混合したケース

a_{n} = a + λ^{n}(c_{1}n^{−α_{1}} + c_{2}n^{−α_{2}} + ⋯)  (10)

(ここで |λ| ≤ 1, λ ≠ 1, α_{1} < α_{2} < ⋯ である。 |λ| = 1 の場合はさらに 0 < α_{1} < α_{2} < ⋯ とする) の数列も考えられる。交替収束 (9) は λ = −1 のケースでもある。

2.2 収束加速法

2.2.1 Richardson 補外

1 次収束で、収束率

λ = lim_{n→∞} \frac{a_{n+1} − a}{a_{n} − a}  (11)

が理論的考察などで分かっている場合は、収束値に近付くと (4) 式から

a_{n+1} − a ≈ λ(a_{n} − a)  (12)

になるので、これを解くと改善した推定値として

a ≈ a_{n}^{(1)} = \frac{a_{n+1} − λa_{n}}{1 − λ} = a_{n} + \frac{a_{n+1} − a_{n}}{1 − λ}  (13)

が得られる (例:台形公式に対する Richardson 補外では λ = 1/4)。さらに、高次の補正が

a_{n} = a + b_{1}λ_{1}^{n} + b_{2}λ_{2}^{n} + ⋯  (14)

であるとき、得られた数列 a_{n}^{(1)} に対し、 Richardson 補外を繰り返し a_{n}^{(2)}, a_{n}^{(3)}, ⋯ 収束性を向上できる。

2.2.2 Aitken 加速

実際の問題では収束率 λ が未知であることが多い。そのような場合、 (4) 式から、

a_{n+1} − a_{n} ≈ bλ^{n+1} − bλ^{n} = bλ^{n}(λ − 1)  (15)

であるので、 a_{n+1} − a_{n} と a_{n+2} − a_{n+1} の比から

λ = lim_{n→∞} \frac{a_{n+2} − a_{n+1}}{a_{n+1} − a_{n}}  (16)

で λ を推定を推定することができる。

この λ の推定式を Richardson 補外で導いた式 (13) に代入すると、

a ≈ a_{n}^{(1)} = a_{n} − \frac{(a_{n+1} − a_{n})^{2}}{a_{n+2} − 2a_{n+1} + a_{n}} = a_{n+2} − \frac{(a_{n+2} − a_{n+1})^{2}}{a_{n+2} − 2a_{n+1} + a_{n}}  (17)

が得られる。

Aitken 加速は 1 次収束に対してだけでなく、混合型 (10) や交代的な収束性を示すもの (9) にも有効である (説明は次回以降)。しかし、単純な冪乗的収束を示すもの (8) には、あまり効果がない。

得られた数列 a_{n}^{(1)} に対し、 Aitken 加速をさらに繰り返し、 a_{n}^{(2)}, a_{n}^{(3)}, ⋯ と収束を向上できる。

\underline{課題：}

与えられた数列に対し、 Aitken 加速された数列を返すプログラムを（できればサブルーチンの形で）作れ。倍精度実数を使うこと。

3 一般的注意

3.1 連絡先

質問などは、講義中に積極的にすること、また、講義後に疑問がある場合は、直接以下の部屋に訪れても良いし、 e-mail を使っての質問も受け付ける。

野村 2641 号室 (knomura@stat.phys.kyushu-u.ac.jp)

3.2 講義のプリント

講義のプリントを次の URL のホームページに置いておいた。 http://maya.phys.kyushu-u.ac.jp/˜ knomura/index.html.ja この中の

講義

計算物理学

から参照すること。また、今後レポート提出者の一覧、レポート再提出の掲示もこの URL を使う予定である。