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

2023年 12月 14日

複素数を利用してみよう

間隔が小さくすると、差分はほとんど同じ値の差を計算するので、誤差が出てくる。

もう一度、テーラー展開を見てみよう。

$$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$$

ここで突然だが、間隔を$h$の代わりに$ih$にしてみよう。$f(x)$は実関数と想定していたが、複素数にも対応する複素関数だとするとどうなるだろうか。

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

この式を実部と虚部に分けると、２つの等式に分けられる。

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

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

$f^{\prime}(x)$は虚部の等式にだけ含まれている。これから、$h^3$以上の項を無視すると

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

が得られる。

Juliaのプログラムは、以前の実数版を少し変更するだけで得られる。

using Printf

diff_complex(f,x; h=1e-20)= imag(f(x+h*im))/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_complex(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.4747595191387584	 -0.0205099779874667
  0.0100000000000000	 -2.4544571054739501	 -0.0002075643226584
  0.0010000000000000	 -2.4542516170381785	 -0.0000020758868868
  0.0001000000000000	 -2.4542495619101850	 -0.0000000207588933
  0.0000100000000000	 -2.4542495413588807	 -0.0000000002075891
  0.0000010000000000	 -2.4542495411533674	 -0.0000000000020757
  0.0000001000000000	 -2.4542495411513117	 -0.0000000000000200
  0.0000000100000000	 -2.4542495411512921	 -0.0000000000000004
  0.0000000010000000	 -2.4542495411512917	  0.0000000000000000
  0.0000000001000000	 -2.4542495411512917	  0.0000000000000000
  0.0000000000100000	 -2.4542495411512912	  0.0000000000000004
  0.0000000000010000	 -2.4542495411512917	  0.0000000000000000
  0.0000000000001000	 -2.4542495411512917	  0.0000000000000000
  0.0000000000000100	 -2.4542495411512912	  0.0000000000000004
  0.0000000000000010	 -2.4542495411512917	  0.0000000000000000
  0.0000000000000001	 -2.4542495411512917	  0.0000000000000000

間隔$h$は今までと同じであるが、順調に減り、浮動小数点数の精度の限界まで正確に計算できているようだ。

この方法を複素ステップ法(complex step method)と呼ぶ。関数評価は1回だけで差分をとっていないので、無駄な桁落ちが回避されている。

実数の世界で考えても上手く行かない場合、複素数に拡張して考えると解決することがある。
実関数の積分が上手くできないとき、複素関数の一周積分に変換して解く方法は、複素関数の本でほぼ確実に紹介されている。
今回紹介した方法は、複素数への拡張の非常に簡単な例である。

超高精度計算

Juliaには、非常に高精度な浮動小数点数としてBigFloat型が用意されている。

精度と、BigFloat型の数を作る方法を以下に示す。

julia> eps(BigFloat)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77

julia> big(π)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

BigFloatのデフォルトの精度は、有効桁数が77である。

今度は、BigFloatを用いて数値微分をやってみよう。

using Printf

diff_complex(f,x; h=1e-20)= imag(f(x+h*im))/h

f(x)  = sin(x^2)
df(x) = 2*x*cos(x^2)
pi_2  = big(π)/2
exact_value = df(pi_2)
println( "f'(π/2)=",exact_value )
println( "  h                                        diff                                                 error")
delta = 0.1
while delta > 1.0e-42
    cdf = diff_complex(f,pi_2,h=delta)
    @printf( "%-7.6g  %-83.80g  %-10.6g\n", delta, cdf, cdf-exact_value )
    delta /= 10.0
end

f'(π/2)=-2.454249541151291862598269654096490678510415782160307138398235372867426448612352
  h                                        diff                                                 error
0.1      -2.474759519138758286790777070616245166951945746783927951166869038241366754018327    -0.02051  
0.01     -2.4544571054739506471593614954667510548849612152921732548875673348131641874160557   -0.000207564
0.001    -2.4542516170381790136213679155693108054057279420359356021244924006031127093114253   -2.07589e-06
0.0001   -2.4542495619101850973776842098786972576451011834973246579325757128514922388872529   -2.07589e-08
1e-05    -2.4542495413588807973823879448965043751497799631185658194185427982574999151330794   -2.07589e-10
1e-06    -2.4542495411533677519463544697664906231659017678225648892726457303563292942749321   -2.07589e-12
1e-07    -2.4542495411513126214917505266186296415580806723192562650358728588165341624036285   -2.07589e-14
1e-08    -2.4542495411512920701872044628241346551605353248906687610443566131772795331326463   -2.07589e-16
1e-09    -2.4542495411512918646741590021837677053364117885363977801080920527026940664759687   -2.07589e-18
1e-10    -2.4542495411512918626190285475773634498762486223437618212913646044164720427386517   -2.07589e-20
1e-11    -2.4542495411512918625984772430312994062106614263976264376842649578076579704662954   -2.07589e-22
1e-12    -2.4542495411512918625982717299858387657870828608340972849435203735412080590634124   -2.07589e-24
1e-13    -2.4542495411512918625982696748553841593831814048907812298024305876642253551546366   -2.07589e-26
1e-14    -2.4542495411512918625982696543040796133191434383876118817486014280386046066591822   -2.07589e-28
1e-15    -2.4542495411512918625982696540985665678585030583950625957778257179230491800517256   -2.07589e-30
1e-16    -2.4542495411512918625982696540965114374038966549236781854409505411156770380977443   -2.07589e-32
1e-17    -2.4542495411512918625982696540964908860993505908879280552128010340469540354568647   -2.07589e-34
1e-18    -2.4542495411512918625982696540964906805863051302475833475663810294792219681132016   -2.07589e-36
1e-19    -2.4542495411512918625982696540964906785311746756411799014894211936124901341538961   -2.07589e-38
1e-20    -2.4542495411512918625982696540964906785106233710951158660291472310748770854677353   -2.07589e-40
1e-21    -2.4542495411512918625982696540964906785104178580496552256745444914495009549809058   -2.07589e-42
1e-22    -2.4542495411512918625982696540964906785104158029192006192709965118962859066726655   -2.07589e-44
1e-23    -2.454249541151291862598269654096490678510415782367896073206963008659677059280514    -2.07589e-46
1e-24    -2.4542495411512918625982696540964906785104157821623830277463226492253489547190167   -2.07589e-48
1e-25    -2.4542495411512918625982696540964906785104157821603278972917162456329120769558945   -2.07589e-50
1e-26    -2.4542495411512918625982696540964906785104157821603073459871701815950932199162853   -2.07589e-52
1e-27    -2.4542495411512918625982696540964906785104157821603071404741247209547032652631545   -2.07589e-54
1e-28    -2.454249541151291862598269654096490678510415782160307138418994266348299215848014    -2.07589e-56
1e-29    -2.4542495411512918625982696540964906785104157821603071383984429618022351762963574   -2.07589e-58
1e-30    -2.4542495411512918625982696540964906785104157821603071383982374487567745358893485   -2.07589e-60
1e-31    -2.4542495411512918625982696540964906785104157821603071383982353936263199294850767   -2.07589e-62
1e-32    -2.4542495411512918625982696540964906785104157821603071383982353730750153834210827   -2.07589e-64
1e-33    -2.454249541151291862598269654096490678510415782160307138398235372869502337960469    -2.07589e-66
1e-34    -2.4542495411512918625982696540964906785104157821603071383982353728674472075058117   -2.07589e-68
1e-35    -2.454249541151291862598269654096490678510415782160307138398235372867426656201317    -2.07589e-70
1e-36    -2.4542495411512918625982696540964906785104157821603071383982353728674264506882447   -2.07589e-72
1e-37    -2.4542495411512918625982696540964906785104157821603071383982353728674264486331475   -2.07959e-74
1e-38    -2.4542495411512918625982696540964906785104157821603071383982353728674264486125589   -2.07268e-76
1e-39    -2.4542495411512918625982696540964906785104157821603071383982353728674264486123516   0         
1e-40    -2.4542495411512918625982696540964906785104157821603071383982353728674264486123516   0         
1e-41    -2.4542495411512918625982696540964906785104157821603071383982353728674264486123516   0         
1e-42    -2.4542495411512918625982696540964906785104157821603071383982353728674264486123516   0

間隔が小さくなるにしたがって、誤差がBigFloat型の限界まで着実に小さくなっていることがわかる。誤差は隠れていて見えないかもしれないので、その場合はスクロールして確認しよう。

BigFloat型は、精度を自由に設定できるので、1000桁までの計算もできるが、この辺で終わりにする。