あざらしとペンギンの問題

主に漫画、数値計算、幾何計算、TCS、一鰭旅、水族館、鰭脚類のことを書きます。

アクセルほど微分精度は加速していくよ!

今回のテーマは微分の計算です。やりたいことは、実数の範囲で定義されたなめらかな関数 f(x) と実数値 a が与えられたとき、fx = a での微分係数 f'(a) を求めることです。

もし f が初等関数の式として与えられているならば話は簡単です。f微分公式を適用していけば必ず導関数 f' が初等関数の式として得られるので、それに x = a を代入して計算すればいいというだけです。もちろんこの計算はコンピュータで実行できます。例えば、MathematicaフリーソフトMaxima のような数式処理システムを使えば導関数の計算から値の計算までをやってくれます。さらには、ソフトウェアをインストールしなくても、Web サービスの Wolfram Alpha を利用することができます。実際 Wolfram alpha は大変良いサービスで、大雑把な入力に対しても大抵欲しい答えを返してくれます。

では f の式が与えられていない場合はどうやって計算すればいいか?精度を高めるにはどうすれば?というのが今回の話の主題です。

単純な方法

定義式をそのまま使う

まず微分の定義からおさらいしましょう。実数の範囲で定義された関数 f について、極限

\displaystyle \lim_{h \to 0} \frac{f(a + h) - f(a)}{h}

が存在するならば、その極限値微分係数といい、f'(a) で表すのでした。

というわけで f'(a) を計算する最も簡単な方法は、小さい h を代入して

\displaystyle D_\mathrm{f}(h) := \frac{f(a + h) - f(a)}{h}

を計算することです。

ここれを関数 f(x) = e^xx = 1 での微分係数を求めるとして試してみましょう。ここで f'(x) = e^x なので真の値は f'(1) = e = 2.718281828459045\dots です。

hh = 1 / 2^3 = 0.125 から初めて 1/2 ずつにしていくとし、それぞれの h について D_\mathrm{f}(h) を計算します。この計算を統計解析ツールの R を使って行うことにします。世の中では汎用プログラム言語である Python が同様の目的で多く使われていますが、パッケージのインストール手段が公式で提供されておらず若干手間がかかるので、今回は R を使うことにします。R はベクトルへの関数適用を、大抵の場合は自動的に要素ごとに分配して、結果をベクトルで返してくれるのでかなり楽です(多くのプログラム言語に存在する map と同様のことを、単に式を書くだけでやってくれる。ループを回すと遅いので、可能な限りベクトルの計算で書くことが望ましい)。インストール方法などは公式サイトを参照してください

公式サイト: R: The R Project for Statistical Computing

R をインストールしたならばそのまま実行します。まずは準備として、次のコマンド群を1行ずつ打ってみてください。ただし、#はコメントなので無視してください。

# R には printf がないので定義しておく
printf = function(...) print(sprintf(...))

# 範囲(1から40まで。R の配列のインデックスは1から始まるので注意)
> R = 1:40

# h のリストを作る
> H = 2^(-(R + 2))

# 関数 f
> f = function(x) exp(x)

#相対誤差を計算する関数
> err = function(x) abs(x - f(1)) / f(1)

# D_f
> dif = function(h) (f(1 + h) - f(1)) / h

# 配列 A に結果を入れる
> A = dif(H)

# 結果を表示
> printf("%e  %.16e  %e", H, A, err(A))

これを実行した結果は次のようになりました。各列は左から順に hD_\mathrm{f}(h)、真の値との相対誤差 |D_\mathrm{f}(h) - e| / e です。

 [1] "1.250000e-01  2.8954801636718877e+00  6.518762e-02"
 [2] "6.250000e-02  2.8050258514034567e+00  3.191134e-02"
 [3] "3.125000e-02  2.7612008889018114e+00  1.578904e-02"
 [4] "1.562500e-02  2.7396294458276031e+00  7.853350e-03"
 [5] "7.812500e-03  2.7289278227360683e+00  3.916442e-03"
 [6] "3.906250e-03  2.7235978923596349e+00  1.955671e-03"
 [7] "1.953125e-03  2.7209381296383981e+00  9.771986e-04"
 [8] "9.765625e-04  2.7196095466729275e+00  4.884402e-04"
 [9] "4.882812e-04  2.7189455795114554e+00  2.441804e-04"
[10] "2.441406e-04  2.7186136769778386e+00  1.220802e-04"
[11] "1.220703e-04  2.7184477459668415e+00  6.103764e-05"
[12] "6.103516e-05  2.7183647855272284e+00  3.051820e-05"
[13] "3.051758e-05  2.7183233065734385e+00  1.525895e-05"
[14] "1.525879e-05  2.7183025674312375e+00  7.629442e-06"
[15] "7.629395e-06  2.7182921979110688e+00  3.814708e-06"
[16] "3.814697e-06  2.7182870132382959e+00  1.907374e-06"
[17] "1.907349e-06  2.7182844209019095e+00  9.537064e-07"
[18] "9.536743e-07  2.7182831247337162e+00  4.768728e-07"
[19] "4.768372e-07  2.7182824769988656e+00  2.385845e-07"
[20] "2.384186e-07  2.7182821538299322e+00  1.196973e-07"
[21] "1.192093e-07  2.7182819917798042e+00  6.008235e-08"
[22] "5.960464e-08  2.7182819098234177e+00  2.993228e-08"
[23] "2.980232e-08  2.7182818800210953e+00  1.896862e-08"
[24] "1.490116e-08  2.7182818651199341e+00  1.348679e-08"
[25] "7.450581e-09  2.7182818651199341e+00  1.348679e-08"
[26] "3.725290e-09  2.7182818651199341e+00  1.348679e-08"
[27] "1.862645e-09  2.7182819843292236e+00  5.734143e-08"
[28] "9.313226e-10  2.7182822227478027e+00  1.450507e-07"
[29] "4.656613e-10  2.7182817459106445e+00  3.036786e-08"
[30] "2.328306e-10  2.7182826995849609e+00  3.204693e-07"
[31] "1.164153e-10  2.7182846069335938e+00  1.022144e-06"
[32] "5.820766e-11  2.7182846069335938e+00  1.022144e-06"
[33] "2.910383e-11  2.7182922363281250e+00  3.828841e-06"
[34] "1.455192e-11  2.7182922363281250e+00  3.828841e-06"
[35] "7.275958e-12  2.7183227539062500e+00  1.505563e-05"
[36] "3.637979e-12  2.7182617187500000e+00  7.397949e-06"
[37] "1.818989e-12  2.7182617187500000e+00  7.397949e-06"
[38] "9.094947e-13  2.7182617187500000e+00  7.397949e-06"
[39] "4.547474e-13  2.7187500000000000e+00  1.722307e-04"
[40] "2.273737e-13  2.7187500000000000e+00  1.722307e-04"

上から見ると h が小さくなるほど誤差が小さくなっているのが分かります。しかし、下の方はむしろ誤差が大きくなっているのが分かるでしょうか?これは桁落ちという現象によるもので、精度を高めようと h をあまり小さくし過ぎると、かえって誤差が大きくなってしまうというジレンマがあるのです。そのため、上の表では小数第7位くらいの精度しか得られていないのです。筆者の環境で R が基本的に採用している64ビット浮動小数点型(C言語でいうところの double)の10進有効精度はおよそ16桁なので、望みうる相対誤差は 1e-16 (= 10^{-16}) 程度です。それと比較すると半分くらいの桁数の精度しか得られていないということになります。

この結果から、微分を精度よく計算するには、h がそれほど小さくないうちから真の値に素早く近づくようにすることが必要だと分かりました。

中点法

微分の精度を高める方法として古くから知られているものに中点法があります。これは、上の D_\mathrm{f}(h) の代わりに

\displaystyle D_\mathrm{m}(h) := \frac{f(a + h) - f(a  - h)}{2h}

を計算するというものです。こうする方が良さそうだということは、直感的には次の図から分かるでしょう。

f:id:azapen6:20151222103901p:plain

それでは、先と同じ例 f(x) = e^xx = 1 で計算してみましょう。先ほどのプログラムの関数 dif を次のように書き換えます。

# D_m
> dif = function(h) (f(1 + h) - f(1 - h)) / (2 * h)

この式を使って R = 1:20 で計算した結果は次のようになりました。各列は左から順に hD_\mathrm{f}(h)、真の値との相対誤差 |D_\mathrm{m}(h) - e| / e で、比較のために4列目に |D_\mathrm{f}(h) - e| / e を記しました。

 [1] "1.250000e-01  2.7253662198037318e+00  2.606202e-03  6.518762e-02"
 [2] "6.250000e-02  2.7200518888706746e+00  6.511688e-04  3.191134e-02"
 [3] "3.125000e-02  2.7187242787455261e+00  1.627684e-04  1.578904e-02"
 [4] "1.562500e-02  2.7183924369799826e+00  4.069060e-05  7.853350e-03"
 [5] "7.812500e-03  2.7183094803361314e+00  1.017256e-05  3.916442e-03"
 [6] "3.906250e-03  2.7182887414124934e+00  2.543133e-06  1.955671e-03"
 [7] "1.953125e-03  2.7182835566964059e+00  6.357830e-07  9.771986e-04"
 [8] "9.765625e-04  2.7182822605182082e+00  1.589457e-07  4.884402e-04"
 [9] "4.882812e-04  2.7182819364734314e+00  3.973627e-08  2.441804e-04"
[10] "2.441406e-04  2.7182818554629193e+00  9.934170e-09  1.220802e-04"
[11] "1.220703e-04  2.7182818352102913e+00  2.483645e-09  6.103764e-05"
[12] "6.103516e-05  2.7182818301480438e+00  6.213479e-10  3.051820e-05"
[13] "3.051758e-05  2.7182818288783892e+00  1.542681e-10  1.525895e-05"
[14] "1.525879e-05  2.7182818285655230e+00  3.917104e-11  7.629442e-06"
[15] "7.629395e-06  2.7182818284782115e+00  7.050942e-12  3.814708e-06"
[16] "3.814697e-06  2.7182818284491077e+00  3.655759e-12  1.907374e-06"
[17] "1.907349e-06  2.7182818284491077e+00  3.655759e-12  9.537064e-07"
[18] "9.536743e-07  2.7182818283326924e+00  4.648256e-11  4.768728e-07"
[19] "4.768372e-07  2.7182818283326924e+00  4.648256e-11  2.385845e-07"
[20] "2.384186e-07  2.7182818287983537e+00  1.248247e-10  1.196973e-07"

先ほどより早く真の値に近づいていることが見て取れます。精度も8桁程度から12桁程度まで改善されました。上の方を見ると誤差の指数の減り方が D_\mathrm{m}(h) の方が2倍ほど速くなっています。下の方では桁落ちによって誤差が大きくなっていますが、それでも D_\mathrm{m}(h) の方が誤差が小さくなっています。

では、どうして中点法の方が定義式にそのまま入れるより精度がよくなるのでしょうか?それを見るには少しばかり計算が必要です。

仮定として f は3回微分可能だとします。このとき、テイラーの定理より定数 A, B が存在して、h \to 0

\displaystyle f(a + h) = f(a) + f'(a) h + A h^2 + B h^3 + o(h^3)

が成り立ちます。これより

\displaystyle D_\mathrm{f}(h) = f'(a) + A h + B h^2 + o(h^2)

となります。すなわち、誤差のオーダーは h の1次です。

一方、中点法では、

\displaystyle f(a + h) = f(a) + f'(a) h + A h^2 + B h^3 + o(h^3)
\displaystyle f(a - h) = f(a) - f'(a) h + A h^2 - B h^3 + o(h^3)
\displaystyle f(a + h) - f(a - h) = 2 f'(a) h + 2 B h^3 + o(h^3)

より

\displaystyle D_\mathrm{m}(h) = f'(a) + B h^2 + o(h^2)

となり、誤差のオーダーは h の2次です。よって、こちらの方が h を小さくしていくと誤差が早く0に近づきます。

これで中点法が定義式をそのまま計算するより理論的にも優れた方法であることが分かりました。

加速法

今までは前フリで本題はここからです。以上の単純な方法よりも高い精度を得たいときどうするかを考えましょう。

一般的な設定から導出

まず、一般的な設定として、関数 A(h)h \to 0

\displaystyle A(h) = A + a_1 h^{p_1} + a_2 h^{p_2} + \dots + a_m h^{p_m} + o(h^{p_m})
ただし p_1 < p_2 < \dots < p_m

という展開を持つとします。このとき、極限 \lim_{h \to 0} A(h) = A を求めたいとします。(ただし、h = 0 で評価することはできないとします。)

先ほどと同じように h を半分ずつにしていくことを考えましょう。このとき、

\displaystyle A(h / 2) = A + \frac{a_1 h^{p_1}}{2^{p_1}} +  \frac{a_2 h^{p_2}}{2^{p_2}} + \dots + \frac{a_m h^{p_m}}{2^{p_m}} + o(h^{p_m})

となります。ここで次の計算をすると h^{p_1} の項を消すことができます。

\displaystyle 2^{p_1} A(h / 2) -  A(h) = (2^{p_1} - 1) A + a'_2 h^{p_2} + \dots + a'_m h^{p_m} + o(h^{p_m})
ただし a'_i \; (i=2, 3, \dots, m) はこの計算における h^{p_i} の係数

さらに A を求めるために変形すると、

\displaystyle A^{(1)}(h) := \frac{2^{p_1} A(h / 2) - A(h)}{2^{p_1} - 1}
\displaystyle = A + a''_2 h^{p_2} + \dots + a''_m h^{p_m} + o(h^{p_m})
ただし a''_i \; (i=2, 3, \dots, m) はこの計算における h^{p_i} の係数

となります。上で定義される関数 A^{(1)}(h)A(h) と同じく h \to 0A に収束します。ここで A(h) では誤差が hp_1 次だったのに対し、A^{(1)}(h) では hp_2 次となるため、より速く A に収束することが期待されます。

このようにして収束を早める方法を加速法といいます。ここで説明した方法はリチャードソン補外法(Richardson extrapolation)と呼ばれています。名前は20世紀初期にこの方法を発明したイギリスの数学者リチャードソン*1から取られていますが、それより200年ほど前に和算家の建部賢弘*2が同様の方法を発明していたことが知られています。

加速法は繰り返し適用することによってさらに収束を早めることができます。A^{(1)}(h) の係数を改めて

\displaystyle A^{(1)}(h) = A + a_2 h^{p_2} + a_3 h^{p_3} + \dots + a_m h^{p_m} + o(h^{p_m})

とします。ここで先と同じように

\displaystyle A^{(2)}(h) := \frac{2^{p_2} A^{(1)}(h / 2) - A^{(1)}(h)}{2^{p_2} - 1}

を定義します。すると

\displaystyle A^{(2)}(h) = A + a'_3 h^{p_3} + \dots + a'_m h^{p_m} + o(h^{p_m})
ただし a'_i \; (i=3, 4, \dots, m) はこの計算における h^{p_i} の係数

となります。誤差は hp_3 次になりました。

さらに

\displaystyle A^{(3)}(h) := \frac{2^{p_3} A^{(2)}(h / 2) - A^{(2)}(h)}{2^{p_3} - 1}

を計算すると、誤差のオーダーはさらに小さくなると考えられます。

計算方法

それでは A^{(j)}(h) を計算する方法を考えましょう。記号の統一のため A^{(0)}(h) := A(h) とします。h は適当に小さすぎない h_0 (先の例では h_0 = 2^{-3} = 0.125)から出発して半分ずつにしていくとします。すなわち h_{k + 1} = h_k / 2 です。このとき、

\displaystyle a^{(j)}_k := A^{(j)}(h_k)

と定義すると、

\displaystyle a^{(j + 1)}_k = A^{(j + 1)}(h_k) = \frac{2^{p_j} A^{(j)}(h_k / 2) - A^{(j)}(h)}{2^{p_j} - 1}
\displaystyle = \frac{2^{p_j} A^{(j)}(h_{k + 1}) - A^{(j)}(h)}{2^{p_j} - 1} = \frac{2^{p_j} a^{(j)}_{k+1} - a^{(j)}_k}{2^{p_j} - 1}

すなわち

\displaystyle a^{(j + 1)}_k = \frac{2^{p_j} a^{(j)}_{k+1} - a^{(j)}_k}{2^{p_j} - 1}

で、a^{(j + 1)}_k を計算するには a^{(j)}_{k+1}, a^{(j)}_k, p_j があればいいということになります。

a^{(j)}_k \; (k = 1, 2, \dots, n) のリストと p = p_j を受け取って a^{(j + 1)}_k \; (k = 1, 2, \dots, n - 1) のリストを計算する関数を次に示します。

# R では配列外の値は NA として扱われる
# 警告は出るかもしれないが落ちることはない、はず
> accel = function(a,p) (2^p * a[R + 1] - a) / (2^p - 1)

未知となっている p_j を求めるには、次の計算をします。

\displaystyle \frac{a^{(j)}_{k + 2} - a^{(j)}_{k + 1}}{a^{(j)}_{k + 1} - a^{(j)}_k} = \frac{(2^{-2p_j} - 2^{-p_j}) a_j h^{p_j} + o(h^{p_j})}{(2^{-p_j} - 1) a_j h^{p^j} + o(h^{p_j})}
\displaystyle \approx \frac{2^{-2p_j} - 2^{-p_j}}{2^{-p_j} - 1} = 2^{-p_j}

これより

\displaystyle -\log_2 \frac{a^{(j)}_{k + 2} - a^{(j)}_{k + 1}}{a^{(j)}_{k + 1} - a^{(j)}_k}

を最も近い整数に丸めたものを p_j として採用すればよいでしょう。上の値を表示する関数を次に示します。

> order = function(a) -log2((a[R + 2] - a[R + 1]) / (a[R + 1] - a))

以下で使うために今までのものを acceldif.R というファイルに保存しておきます。

微分を加速する

それでは本題の微分を加速しましょう。目標は f(x) = e^x について f'(1) を求めることです。加速する元の数列は a^{(0)}_k = D_\mathrm{m}(h_k)h_k = 2^{-k - 3} です。

中点法の議論と同様にして、f2m + 1微分可能なとき

\displaystyle D_\mathrm{m}(h) = f'(a) + a_2 h^2 + a_4 h^4 + \dots + a_{2m} h^{2m} + o(h^{2m})

となります。これより、p_1 = 2 となることが予想されます。これを先の関数 order で確かめてみましょう。acceldif.R をインポートした段階では、変数 A0 に D_\mathrm{m}(h) のリストが入っています。

> source("acceldif.R")

> order(A0)
 [1] 2.001057 2.000264 2.000066 2.000017 2.000004 2.000001 1.999999 2.000012
 [9] 2.000000 2.000259 1.995342 2.020818 1.841302 1.584963      Inf      NaN
[17]      Inf     -Inf       NA       NA

これより、最初の方ではちゃんと期待通りの値が出ていることが分かります。よって、まずは p_1 = 2 として accel を呼び、a^{(1)}_k のリストを計算します。

# A0 を加速
> A1 = accel(A0, 2)

> print_results(A1)
 [1] "1.250000e-01  2.7182804452263221e+00  5.088629e-07"
 [2] "6.250000e-02  2.7182817420371435e+00  3.179284e-08"
 [3] "3.125000e-02  2.7182818230581347e+00  1.986884e-09"
 [4] "1.562500e-02  2.7182818281215142e+00  1.241707e-10"
 [5] "7.812500e-03  2.7182818284379473e+00  7.761443e-12"
 [6] "3.906250e-03  2.7182818284577102e+00  4.910941e-13"
 [7] "1.953125e-03  2.7182818284588088e+00  8.691353e-14"
 [8] "9.765625e-04  2.7182818284585060e+00  1.983327e-13"
 [9] "4.882812e-04  2.7182818284594155e+00  1.362517e-13"
[10] "2.441406e-04  2.7182818284594155e+00  1.362517e-13"
[11] "1.220703e-04  2.7182818284606278e+00  5.822553e-13"
[12] "6.103516e-05  2.7182818284551709e+00  1.425251e-12"
[13] "3.051758e-05  2.7182818284612345e+00  8.054205e-13"
[14] "1.525879e-05  2.7182818284491077e+00  3.655759e-12"
[15] "7.629395e-06  2.7182818284394066e+00  7.224605e-12"
[16] "3.814697e-06  2.7182818284491077e+00  3.655759e-12"
[17] "1.907349e-06  2.7182818282938874e+00  6.075811e-11"
[18] "9.536743e-07  2.7182818283326924e+00  4.648256e-11"
[19] "4.768372e-07  2.7182818289535740e+00  1.819270e-10"

3列目の誤差を見ると、この時点でかなり速く真の値に近づいていることが分かります。

今選んだ p_1 が適切であることを見るために、異なる p_1 についてもやってみます。

> B = accel(A0, 1)

> print_results(B)
 [1] "1.250000e-01  2.7147375579376174e+00  1.303864e-03"
 [2] "6.250000e-02  2.7173966686203777e+00  3.256321e-04"
 [3] "3.125000e-02  2.7180605952144390e+00  8.138716e-05"
 [4] "1.562500e-02  2.7182265236922802e+00  2.034549e-05"
 [5] "7.812500e-03  2.7182680024888555e+00  5.086290e-06"
 [6] "3.906250e-03  2.7182783719803183e+00  1.271567e-06"
 [7] "1.953125e-03  2.7182809643400105e+00  3.178916e-07"
 [8] "9.765625e-04  2.7182816124286546e+00  7.947314e-08"
 [9] "4.882812e-04  2.7182817744524073e+00  1.986793e-08"
[10] "2.441406e-04  2.7182818149576633e+00  4.966881e-09"
[11] "1.220703e-04  2.7182818250857963e+00  1.240949e-09"
[12] "6.103516e-05  2.7182818276087346e+00  3.128117e-10"
[13] "3.051758e-05  2.7182818282526569e+00  7.592599e-11"
[14] "1.525879e-05  2.7182818283909000e+00  2.506916e-11"
[15] "7.629395e-06  2.7182818284200039e+00  1.436246e-11"
[16] "3.814697e-06  2.7182818284491077e+00  3.655759e-12"
[17] "1.907349e-06  2.7182818282162771e+00  8.930937e-11"
[18] "9.536743e-07  2.7182818283326924e+00  4.648256e-11"
[19] "4.768372e-07  2.7182818292640150e+00  2.961319e-10"

誤差を見ると今ひとつ収束が悪いのが分かるでしょう。

話を戻して、先程求めた A1 に対して再び order を計算し、p_2 の値を決めます。理論的には p_2 = 4 となるはずです。

> order(A1)
 [1]  4.0005285  4.0001232  4.0001282  4.0010370  4.1689528        NaN

早々とエラーになったのは誤差が非常に近接しているためです。ともかく p_2 = 4 が決まりました。先ほどと同じように accel を呼び、a^{(2)}_k のリストを計算します。以下では見やすさのため最初の10個だけ表示します。

> A2 = accel(A1, 4)

> print_results(A2)
 [1] "1.250000e-01  2.7182818284911985e+00  1.182857e-11"
 [2] "6.250000e-02  2.7182818284595340e+00  1.798718e-13"
 [3] "3.125000e-02  2.7182818284590726e+00  1.012902e-14"
 [4] "1.562500e-02  2.7182818284590429e+00  8.168565e-16" # 真の値は 2.718281828459045...
 [5] "7.812500e-03  2.7182818284590278e+00  6.371480e-15"
 [6] "3.906250e-03  2.7182818284588821e+00  5.995726e-14"
 [7] "1.953125e-03  2.7182818284584860e+00  2.056845e-13"
 [8] "9.765625e-04  2.7182818284594759e+00  1.584702e-13"
 [9] "4.882812e-04  2.7182818284594155e+00  1.362517e-13"
[10] "2.441406e-04  2.7182818284607086e+00  6.119889e-13"

2番目の時点で既に中点法で得られた最良の値よりも誤差が小さくなっています。さらに、4番目は64ビット浮動小数点型でで到達できる精度の限界にかなり迫っています。一方で、加速するほど桁落ちのために列を先に進めるとすぐ精度が落ちてしまうという問題も見えてきます。

さらに加速を試みて order を計算すると、

> order(A2)
 [1]  6.1006750  3.9548907  0.9786263 -3.2700892 -1.4433479        NaN

となり、理論上は p_3 = 6 となるところが結果が安定しません。というのも、A2 の収束が非常に速いため、初期で著しい桁落ちが生じてしまうのです。ここまで来るとそろそろやめた方がいいのですが、とりあえずもう一度くらい加速してみましょう。

> A3 = accel(A2, 6)

> print_results(A3)
 [1] "1.250000e-01  2.7182818284590313e+00  5.064510e-15"
 [2] "6.250000e-02  2.7182818284590651e+00  7.351708e-15"
 [3] "3.125000e-02  2.7182818284590424e+00  9.802277e-16"
 [4] "1.562500e-02  2.7182818284590273e+00  6.534852e-15"
 [5] "7.812500e-03  2.7182818284588799e+00  6.077412e-14"
 [6] "3.906250e-03  2.7182818284584798e+00  2.079717e-13"
 [7] "1.953125e-03  2.7182818284594914e+00  1.641881e-13"
 [8] "9.765625e-04  2.7182818284594146e+00  1.359249e-13"
 [9] "4.882812e-04  2.7182818284607291e+00  6.195039e-13"
[10] "2.441406e-04  2.7182818284547134e+00  1.593524e-12"
[11] "1.220703e-04  2.7182818284617474e+00  9.941143e-13"
[12] "6.103516e-05  2.7182818284480872e+00  4.031187e-12"
[13] "3.051758e-05  2.7182818284386085e+00  7.518183e-12"
[14] "1.525879e-05  2.7182818284499288e+00  3.353686e-12"
[15] "7.629395e-06  2.7182818282809009e+00  6.553558e-11"
[16] "3.814697e-06  2.7182818283361003e+00  4.522885e-11"
[17] "1.907349e-06  2.7182818290054374e+00  2.010065e-10"

なんだかんだで最初の時点でもうほとんど有効精度に近い値が出ています。これでさらに order を計算すると、

> order(A3)
 [1]        NaN  0.5849625 -3.2875766 -1.4403439        NaN        NaN

はい、終わり!

改めて A0, A1, A2, A3 の誤差だけを横に並べてみます。

> R = 1:10
> err(cbind(A0[R], A1[R], A2[R], A3[R]))
 [1,] 2.606202e-03 5.088629e-07 1.182857e-11 5.064510e-15
 [2,] 6.511688e-04 3.179284e-08 1.798718e-13 7.351708e-15
 [3,] 1.627684e-04 1.986884e-09 1.012902e-14 9.802277e-16
 [4,] 4.069060e-05 1.241707e-10 8.168565e-16 6.534852e-15
 [5,] 1.017256e-05 7.761443e-12 6.371480e-15 6.077412e-14
 [6,] 2.543133e-06 4.910941e-13 5.995726e-14 2.079717e-13
 [7,] 6.357830e-07 8.691353e-14 2.056845e-13 1.641881e-13
 [8,] 1.589457e-07 1.983327e-13 1.584702e-13 1.359249e-13
 [9,] 3.973627e-08 1.362517e-13 1.362517e-13 6.195039e-13
[10,] 9.934170e-09 1.362517e-13 6.119889e-13 1.593524e-12

最初の5行を見れば加速するほど真の値に近いことが分かります。最も精度が高い A2[4] は A1[4], A1[5] に依存し、さらには A0[4], A0[5], A0[6] に依存しているので、たった3個の精度が低い値からこれだけの精度が出ているのです。これで加速法の凄さが分かって頂けたでしょうか?

まとめ

今回は微分を精度よく計算する方法について考えました。定義式をそのまま計算するよりも中点法を使った法が精度よく計算することを見ました。さらに、リチャードソン補外法によって収束を加速する方法について見ました。加速法を用いると、元の数列の最初の方を見るだけでも高い精度が得られることが分かりました。一方で、桁落ちのために加速した数列を進めることは必ずしもよくないことも分かりました。

リチャードソン補外法を用いて微分の精度を上げる方法はロンバーグ微分公式と呼ばれます。また、リチャードソン補外法を積分の台形公式に適用したロンバーグ積分公式もあります。

今回述べたような加速法は関数の滑らかさや数列の性質に依存して得意不得意が決まるので、万能の処方箋ではありません。実際に使うときはそのことに気をつけてください。

R は筆者がかなり前から使っているツールですが、Python と違ってそもそも統計解析向けに作られているので、ベクトル・行列計算用電卓として使う程度の計算はパッケージのインストールなしに行うことができます。パッケージのインストールが R の上でコマンドを叩くだけで行うことができるのも楽です。同様の機能を提供する比較的新しいプログラム言語として Julia というものもあります。とりあえずいくつか使ってみた限りでは、筆者に馴染みのある R の方がベクトル・行列計算用電卓としては使い勝手が良いように思われました。

というわけで、加速法の話しでした。では。

*1:Lewis Fry Richardson : 1881-1953

*2:建部賢弘(たけべ かたひろ/けんこう):1664-1739