読者です 読者をやめる 読者になる 読者になる

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

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

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

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

もし f が初等関数の式として与えられているならば話は簡単です。f微分公式を適用していけば必ず導関数 f' が初等関数の式として得られるので、それに x = a を代入して計算すればいいというだけです。もちろんこの計算はコンピュータで実行できます。例えば、MathematicaフリーソフトMaxima のような数式処理システムを使えば導関数の計算から値の計算までをやってくれます。さらには、ソフトウェアをインストールしなくても、Web サービスの 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 での微分係数を求めるとして試してみましょう。*1ここで f'(x) = e^x なので真の値は f'(1) = e = 2.718281828459045\dots です。

hh = 1 / 2^3 = 0.125 から初めて 1/2 ずつにしていくとし、それぞれの h について D_\mathrm{f}(h) を計算します。この計算を Julia というこの記事が書かれた時点ではまだ新しいプログラム言語で書いたコードを次に示します。

公式サイト: http://julialang.org/

Julia は科学技術計算での使用を第一に考慮して作られたプログラム言語で、MATLAB に似た簡潔な行列記法と豊富な数値計算関数が標準で入っているのが特徴です。言語としては PythonRuby に近い動的言語であり、簡潔な記法でプログラミングを行うことが可能であると同時に、それらの言語になかった高度なパフォーマンスの向上が謳われています。筆者も最近使い始めたばかりなので、自分の練習も兼ねて Julia で書いてみることにしました。Julia のコマンドプロンプトで次のコード(#はコメントなので除く)を1行ずつ打っていきます。

#範囲 (Julia ではインデックスは基本的に1から始まり範囲は終端を含むので注意)
R = 1:40

#h のリスト
h = [2.0^(-k) for k in R + 2]

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

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

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

#配列 a に結果を入れる
a = map(dif, h)

for k in R
   @printf("%e\t%.16e\t%e\n", h[k], a[k], err(a[k]))
end

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

1.250000e-01    2.8954801636718877e+00  6.518762e-02
6.250000e-02    2.8050258514034567e+00  3.191134e-02
3.125000e-02    2.7612008889018114e+00  1.578904e-02
1.562500e-02    2.7396294458276031e+00  7.853350e-03
7.812500e-03    2.7289278227360683e+00  3.916442e-03
3.906250e-03    2.7235978923596349e+00  1.955671e-03
1.953125e-03    2.7209381296383981e+00  9.771986e-04
9.765625e-04    2.7196095466729275e+00  4.884402e-04
4.882813e-04    2.7189455795114554e+00  2.441804e-04
2.441406e-04    2.7186136769778386e+00  1.220802e-04
1.220703e-04    2.7184477459668415e+00  6.103764e-05
6.103516e-05    2.7183647855272284e+00  3.051820e-05
3.051758e-05    2.7183233065734385e+00  1.525895e-05
1.525879e-05    2.7183025674312375e+00  7.629442e-06
7.629395e-06    2.7182921979110688e+00  3.814708e-06
3.814697e-06    2.7182870132382959e+00  1.907374e-06
1.907349e-06    2.7182844209019095e+00  9.537064e-07
9.536743e-07    2.7182831247337162e+00  4.768728e-07
4.768372e-07    2.7182824769988656e+00  2.385845e-07
2.384186e-07    2.7182821538299322e+00  1.196973e-07
1.192093e-07    2.7182819917798042e+00  6.008235e-08
5.960464e-08    2.7182819098234177e+00  2.993228e-08
2.980232e-08    2.7182818800210953e+00  1.896862e-08
1.490116e-08    2.7182818651199341e+00  1.348679e-08
7.450581e-09    2.7182818651199341e+00  1.348679e-08
3.725290e-09    2.7182818651199341e+00  1.348679e-08
1.862645e-09    2.7182819843292236e+00  5.734143e-08
9.313226e-10    2.7182822227478027e+00  1.450507e-07
4.656613e-10    2.7182817459106445e+00  3.036786e-08
2.328306e-10    2.7182826995849609e+00  3.204693e-07
1.164153e-10    2.7182846069335938e+00  1.022144e-06
5.820766e-11    2.7182846069335938e+00  1.022144e-06
2.910383e-11    2.7182922363281250e+00  3.828841e-06
1.455192e-11    2.7182922363281250e+00  3.828841e-06
7.275958e-12    2.7183227539062500e+00  1.505563e-05
3.637979e-12    2.7182617187500000e+00  7.397949e-06
1.818989e-12    2.7182617187500000e+00  7.397949e-06
9.094947e-13    2.7182617187500000e+00  7.397949e-06
4.547474e-13    2.7187500000000000e+00  1.722307e-04
2.273737e-13    2.7187500000000000e+00  1.722307e-04

上から見ると h が小さくなるほど誤差が小さくなっているのが分かります。しかし、下の方はむしろ誤差が大きくなっているのが分かるでしょうか?これは桁落ちという現象によるもので、精度を高めようと h をあまり小さくし過ぎると、かえって誤差が大きくなってしまうというジレンマがあるのです。そのため、上の表では小数第7位くらいの精度しか得られていないのです。Julia が基本的に採用している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(h) = (f(1 + h) - f(1 - h)) / 2h
end

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

1.250000e-01    2.7253662198037318e+00  2.606202e-03    6.518762e-02
1.250000e-01    2.7253662198037318e+00  2.606202e-03    6.518762e-02
6.250000e-02    2.7200518888706746e+00  6.511688e-04    3.191134e-02
3.125000e-02    2.7187242787455261e+00  1.627684e-04    1.578904e-02
1.562500e-02    2.7183924369799826e+00  4.069060e-05    7.853350e-03
7.812500e-03    2.7183094803361314e+00  1.017256e-05    3.916442e-03
3.906250e-03    2.7182887414124934e+00  2.543133e-06    1.955671e-03
1.953125e-03    2.7182835566964059e+00  6.357830e-07    9.771986e-04
9.765625e-04    2.7182822605182082e+00  1.589457e-07    4.884402e-04
4.882813e-04    2.7182819364734314e+00  3.973627e-08    2.441804e-04
2.441406e-04    2.7182818554629193e+00  9.934170e-09    1.220802e-04
1.220703e-04    2.7182818352102913e+00  2.483645e-09    6.103764e-05
6.103516e-05    2.7182818301480438e+00  6.213479e-10    3.051820e-05
3.051758e-05    2.7182818288783892e+00  1.542681e-10    1.525895e-05
1.525879e-05    2.7182818285655230e+00  3.917104e-11    7.629442e-06
7.629395e-06    2.7182818284782115e+00  7.050942e-12    3.814708e-06
3.814697e-06    2.7182818284491077e+00  3.655759e-12    1.907374e-06
1.907349e-06    2.7182818284491077e+00  3.655759e-12    9.537064e-07
9.536743e-07    2.7182818283326924e+00  4.648256e-11    4.768728e-07
4.768372e-07    2.7182818283326924e+00  4.648256e-11    2.385845e-07
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世紀初期にこの方法を発明したイギリスの数学者リチャードソン*2から取られていますが、それより200年ほど前に和算家の建部賢弘*3が同様の方法を発明していたことが知られています。

加速法は繰り返し適用することによってさらに収束を早めることができます。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) のリストを計算する関数を次に示します。

accel(a, p) = [(2^p * a[k + 1] -  a[k]) / (2^p - 1) for k in 1:(length(a) - 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 として採用すればよいでしょう。上の値を表示する関数を次に示します。

function order(a)
    for k in 1:(length(a) - 2)
        println(log2((a[k + 2] - a[k + 1]) / (a[k + 1] - a[k])))
    end
end

敢えて値を表示するようにしているのは、あまりに差が小さくなると桁落ちや0除算で上手くいかないので、ちゃんと値が出ているところを見極めるためです。

以下で使うために今までのものを acceldif.jl というファイルに保存しておきます。結果を表示するための関数も print_results として用意しました。

微分を加速する

それでは本題の微分を加速しましょう。目標は 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.py をインポートした段階では、変数 a に D_\mathrm{m}(h) のリストが入っています。

include("acceldif.jl")

order(a)
-2.0010565324454466
-2.000264157259162
-2.0000660412339974
-2.0000165093134608
-2.000004124330373
-2.000000917398196
-1.9999989876985755
-2.000012147663972
-2.0
-2.000259174534226
-1.9953419671219552
-2.020818471507554
-1.8413022539809418
-1.5849625007211563
-Inf
#エラー

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

#a を加速
a1 = accel(a, 2)

print_results(a1)
1.250000e-01    2.7182804452263221e+00  5.088629e-07
6.250000e-02    2.7182817420371435e+00  3.179284e-08
3.125000e-02    2.7182818230581347e+00  1.986884e-09
1.562500e-02    2.7182818281215142e+00  1.241707e-10
7.812500e-03    2.7182818284379473e+00  7.761443e-12
3.906250e-03    2.7182818284577102e+00  4.910941e-13
1.953125e-03    2.7182818284588088e+00  8.691353e-14
9.765625e-04    2.7182818284585060e+00  1.983327e-13
4.882813e-04    2.7182818284594155e+00  1.362517e-13
2.441406e-04    2.7182818284594155e+00  1.362517e-13
1.220703e-04    2.7182818284606278e+00  5.822553e-13
6.103516e-05    2.7182818284551709e+00  1.425251e-12
3.051758e-05    2.7182818284612345e+00  8.054205e-13
1.525879e-05    2.7182818284491077e+00  3.655759e-12
7.629395e-06    2.7182818284394066e+00  7.224605e-12
3.814697e-06    2.7182818284491077e+00  3.655759e-12
1.907349e-06    2.7182818282938874e+00  6.075811e-11
9.536743e-07    2.7182818283326924e+00  4.648256e-11
4.768372e-07    2.7182818289535740e+00  1.819270e-10

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

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

b = accel(a, 1)

print_results(b)
1.250000e-01    2.7192926987373807e+00  3.718784e-04
6.250000e-02    2.7185346201562193e+00  9.299687e-05
3.125000e-02    2.7183450310134765e+00  2.325092e-05
1.562500e-02    2.7182976293870098e+00  5.812837e-06
7.812500e-03    2.7182857787091166e+00  1.453216e-06
3.906250e-03    2.7182828160226791e+00  3.633044e-07
1.953125e-03    2.7182820753498942e+00  9.082607e-08
9.765625e-04    2.7182818901813204e+00  2.270636e-08
4.882813e-04    2.7182818438899892e+00  5.676727e-09
2.441406e-04    2.7182818323170586e+00  1.419284e-09
1.220703e-04    2.7182818294248654e+00  3.553054e-10
6.103516e-05    2.7182818286970098e+00  8.754234e-11
3.051758e-05    2.7182818285208277e+00  2.272854e-11
1.525879e-05    2.7182818284657384e+00  2.462332e-12
7.629395e-06    2.7182818284449501e+00  5.185241e-12
3.814697e-06    2.7182818284491077e+00  3.655759e-12
1.907349e-06    2.7182818283160617e+00  5.260065e-11
9.536743e-07    2.7182818283326924e+00  4.648256e-11
4.768372e-07    2.7182818288648769e+00  1.492972e-10

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

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

order(a1)
-4.000528488784751
-4.00012321429058
-4.000128183764117
-4.001037024199116
-4.168952769421632
#エラー

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

a2 = accel(a1, 4)

print_results(a2)
1.250000e-01    2.7182818284911985e+00  1.182857e-11
6.250000e-02    2.7182818284595340e+00  1.798718e-13
3.125000e-02    2.7182818284590726e+00  1.012902e-14
1.562500e-02    2.7182818284590429e+00  8.168565e-16 #真の値は 2.718281828459045...
7.812500e-03    2.7182818284590278e+00  6.371480e-15
3.906250e-03    2.7182818284588821e+00  5.995726e-14
1.953125e-03    2.7182818284584860e+00  2.056845e-13
9.765625e-04    2.7182818284594759e+00  1.584702e-13
4.882813e-04    2.7182818284594155e+00  1.362517e-13
2.441406e-04    2.7182818284607086e+00  6.119889e-13
1.220703e-04    2.7182818284548071e+00  1.559052e-12

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

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

order(a2)
-6.100674985050926
-3.9548907484464393
-0.978626349207433
3.2700891633677442
1.4433478953022212
#エラー

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

a3 = accel(a2, 6)

print_results(a3)
1.250000e-01    2.7182818284590313e+00  5.064510e-15
6.250000e-02    2.7182818284590651e+00  7.351708e-15
3.125000e-02    2.7182818284590424e+00  9.802277e-16
1.562500e-02    2.7182818284590273e+00  6.534852e-15
7.812500e-03    2.7182818284588799e+00  6.077412e-14
3.906250e-03    2.7182818284584798e+00  2.079717e-13
1.953125e-03    2.7182818284594914e+00  1.641881e-13
9.765625e-04    2.7182818284594146e+00  1.359249e-13
4.882813e-04    2.7182818284607291e+00  6.195039e-13

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

order(a3)
#エラー

はい、終わり!

改めて a, a1, a2, a3 の誤差だけを横に並べてみます。

R = 1:10
map(err, [a[R] a1[R] a2[R] a3[R]])
10x4 Array{Float64,2}:
 0.0026062    5.08863e-7   1.18286e-11  5.06451e-15
 0.000651169  3.17928e-8   1.79872e-13  7.35171e-15
 0.000162768  1.98688e-9   1.0129e-14   9.80228e-16
 4.06906e-5   1.24171e-10  8.16856e-16  6.53485e-15
 1.01726e-5   7.76144e-12  6.37148e-15  6.07741e-14
 2.54313e-6   4.91094e-13  5.99573e-14  2.07972e-13
 6.35783e-7   8.69135e-14  2.05684e-13  1.64188e-13
 1.58946e-7   1.98333e-13  1.5847e-13   1.35925e-13
 3.97363e-8   1.36252e-13  1.36252e-13  6.19504e-13
 9.93417e-9   1.36252e-13  6.11989e-13  1.59352e-12

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

まとめ

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

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

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

ついでにプログラム言語 Julia 筆者は結構気に入りました。ただし、環境やバージョンによって精度に違いが出るという問題もあるようです。元々環境に最適化して高度なパフォーマンスを発揮することを目指しているとのことなので、ビルドに使うコンパイラやFPUの拡張精度の使用およびSIMDの使用などで違いが出るのかもしれません。その辺りはもっと使ってみてから判断すればいいでしょう。では。

*1:今回の例に多項式関数を用いなかったのは、多項式関数では1回加速しただけで厳密に真の値が求まってしまうためです。

*2:Lewis Fry Richardson : 1881-1953

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