数値微分の高精度化（上）

テーラー展開と数値微分

関数 $f(x)$ のテーラー展開は次のようになる。

$$f(x+h) = f(x) + \frac{f^{\prime}(x)}{1!} h+ \frac{f^{\prime\prime}(x)}{2!}h^2
+ \frac{f^{\prime\prime\prime}(x)}{2!}h^3 + \cdots$$

これを並べ替えて整理すると、

$$ f^{\prime}(x) = \frac{f(x+h) – f(x)}{h} – \frac{f^{\prime\prime}(x)}{2!}h
– \frac{f^{\prime\prime\prime}(x)}{2!}h^2 – \cdots$$

$$ f^{\prime}(x) \approx \frac{f(x+h) – f(x)}{h}$$

というお馴染みの1階微分の近似式が得られ、これに基づいて数値微分をやってみよう。
この数値微分の方法を前進差分という。

今回は、数値計算にも非常に優れているプログラミング言語Juliaを用いてみる。
関数$f(x)$として$sin(x^2)$として、$x=\pi/2$の時の値を実験してみる。

計算には、デフォルトの倍精度浮動小数点型のFloat64が使われる。

間隔$h$を $h \to 0$ とすることで精度が上がっていくのだが、Float64の精度の限界が$2.2e^{-16}$なので、間隔を0.1から順に1/10にし、$1.0^{e-17}$　より大きい間だけ繰り返している。

using Printf

diff_forward(f,x; h=1e-20) = (f(x+h) - f(x))/h

f(x)  = sin(x^2)
df(x) = 2*x*cos(x^2)
pi_2  = π/2
exact_value = df(pi_2)
println( "f'(π/2)=",exact_value )
println( "            h                     diff                   error")
delta = 0.1
while delta > 1.0e-17
    cdf = diff_forward(f,pi_2,h=delta)
    @printf( "%20.16f\t%20.16f\t%20.16f\n", delta, cdf, cdf-exact_value )
    delta /= 10.0
end

f'(π/2)=-2.4542495411512917
            h                     diff                   error
  0.1000000000000000	 -2.8133781493584644	 -0.3591286082071727
  0.0100000000000000	 -2.4926542600235968	 -0.0384047188723051
  0.0010000000000000	 -2.4581093000713983	 -0.0038597589201066
  0.0001000000000000	 -2.4546357044752387	 -0.0003861633239470
  0.0000100000000000	 -2.4542881593814236	 -0.0000386182301320
  0.0000010000000000	 -2.4542534027283121	 -0.0000038615770204
  0.0000001000000000	 -2.4542499288404680	 -0.0000003876891763
  0.0000000100000000	 -2.4542495591362008	 -0.0000000179849091
  0.0000000010000000	 -2.4542496923629638	 -0.0000001512116721
  0.0000000001000000	 -2.4542501364521732	 -0.0000005953008815
  0.0000000000100000	 -2.4542368137758781	  0.0000127273754136
  0.0000000000010000	 -2.4545920851437582	 -0.0003425439924665
  0.0000000000001000	 -2.4535928844215955	  0.0006566567296962
  0.0000000000000100	 -2.4646951146678471	 -0.0104455735165554
  0.0000000000000010	 -2.7755575615628914	 -0.3213080204115997
  0.0000000000000001	  0.0000000000000000	  2.4542495411512917

右のerrorの欄が正確な値 $f\prime(\pi/2) = 2x cos(\pi/2)$ との誤差である。
$h=1.0e^{-8}$ までは誤差が減少しているが、それ以上になると、誤差が急増しているので、間隔をやたらに小さくするのは不味い。

中心差分

これを改良する方法として良く紹介されているのが中心差分による方法である。

$$ f^{\prime}(x) \approx \frac{f(x+h/2) – f(x-h/2)}{h}$$

微分したい$x$の値に対して、$x$を中心として両側に$h/2$だけ離れた点の関数値を求めることで、誤差が減ることが期待される。

using Printf

diff_central(f,x; h=1e-20) = (f(x+h/2) - f(x-h/2))/h

f(x)  = sin(x^2)
df(x) = 2*x*cos(x^2)
pi_2  = π/2
exact_value = df(pi_2)
println( "f'(π/2)=",exact_value )
println( "            h                     diff                   error")
delta = 0.1
while delta > 1.0e-17
    cdf = diff_central(f,pi_2,h=delta)
    @printf( "%20.16f\t%20.16f\t%20.16f\n", delta, cdf, cdf-exact_value )
    delta /= 10.0
end

f'(π/2)=-2.4542495411512917
            h                     diff                   error
  0.1000000000000000	 -2.4490444809132828	  0.0052050602380089
  0.0100000000000000	 -2.4541976423794742	  0.0000518987718174
  0.0010000000000000	 -2.4542490221782787	  0.0000005189730130
  0.0001000000000000	 -2.4542495359680672	  0.0000000051832245
  0.0000100000000000	 -2.4542495410950771	  0.0000000000562146
  0.0000010000000000	 -2.4542495415946770	 -0.0000000004433853
  0.0000001000000000	 -2.4542495369317403	  0.0000000042195514
  0.0000000100000000	 -2.4542495591362008	 -0.0000000179849091
  0.0000000010000000	 -2.4542495813406613	 -0.0000000401893696
  0.0000000001000000	 -2.4542501364521732	 -0.0000005953008815
  0.0000000000100000	 -2.4542368137758781	  0.0000127273754136
  0.0000000000010000	 -2.4544810628412956	 -0.0002315216900040
  0.0000000000001000	 -2.4524826613969704	  0.0017668797543213
  0.0000000000000100	 -2.4980018054066018	 -0.0437522642553101
  0.0000000000000010	 -2.1094237467877974	  0.3448257943634943
  0.0000000000000001	  0.0000000000000000	  2.4542495411512917

前進差分と比較すると、hが1/10、誤差が1/100のオーダーで良くなったみたいだ。

しかし、間隔を狭めると誤差の方が目立つ傾向に変わりがない。

ここまでの説明は、数値計算に関する教科書やネット上でいっぱい説明されていることである。
間隔を狭め過ぎた場合に誤差がどのように増えていくか。誤差が最小になるのは、間隔がどのあたりかを求める方法なども紹介されているので参考にされたい。

計算誤差との闘い

数値計算をやっていると、必ず計算誤差に悩まされる。
計算結果が信用できるか、精度はどこまで保証されているかが常に問題になる。

差分はよく利用するが、非常に怖い特徴をもっている。

$f(x+h/2) – f(x-h/2)$ の部分は、ほとんど等しい値である$\frac{f(x+h/2)$ と $\frac{f(x+h/2)$の差を求めているが、ほとんど等しい浮動小数点数の差を求めるとき、桁落ちが起きる。桁落ちの桁数分だけ有効桁数が減ってしまい、本来16桁の精度があっても、差の精度は半分以下になることがある。

上記の計算結果でも、間隔をどんどん少なくすると、diffが0になっているが、これは分子が０になっているからである。

たとえば、

julia> π/2　+ 1.0e-16 == π/2
true

のように、π/2 に非常に小さな値を足しても、完全に無視されて結果は元のπ/2のまま（完全一致）になる。

差を利用している限り、このジレンマからは抜けられそうにない。何か良い方法はないだろか。