本当は怖いHPC

HPC屋の趣味&実益ブログ

数値微分(Juliaで)

Juliaで数値微分を書いてみます。

数値微分の基本

数値微分は、直感的には微小幅  h について関数の値を計算して差分を取れば、それが勾配になります。しかし、数学的にもう少しちゃんと根拠が示されています。

まず、微分を求めたい1変数関数を  f、点を x_0とします。そして、微小幅  h > 0 について、関数  f(x + h) x の周りでテイラー展開します。

 { \displaystyle
f(x + h) = f(x) + h f^{(1)}(x) + \frac{h^{2}}{2} f^{(2)}(x)  + \frac{h^{3}}{6} f^{(3)}(x) + \dots
}

ただし、 \dots は、 h^{4} 以上の項です。同様に、関数  f(x - h)xの周りでテイラー展開します。

 { \displaystyle
f(x - h) = f(x) - h f^{(1)}(x) + \frac{h^{2}}{2} f^{(2)}(x)  - \frac{h^{3}}{6} f^{(3)}(x) + \dots
}

ここで上2つの式を引き算すると、都合の良いことに  h^{2} の項が消えます。

 {\displaystyle
f(x + h) - f(x - h) = 2 h f^{(1)}(x) + \frac{h^{3}}{6} f^{(3)}(x) + \dots
}

 h^{3} 以上の項を無視することによって、

{\displaystyle
f^{(1)}(x)  = \frac{ f(x + h) - f(x - h)  }{2h}
}

という結果が得られました。これは中心差分の公式と呼ばれます。上の式変形からわかるように、 h^{1}の項まで考慮しているので(計算の途中で消えたけどw)、誤差のオーダーが  O(h^{2})(二次精度)となります。

Juliaで実装

上式をJuliaで実装してみます。Float32をベースとして使って、 h としてJuliaの組み込み定数であるeps(Float16)を使ってみます。

コード自体は非常に単純です。

function diff(f)
  h = eps(Float16)
  return (x -> (f(x+h) - f(x-h))/(2h))
end

受け取った関数  f の1次の導関数を返します。試しに、いくつか解析的に自明な関数の計算をしてみます。

julia> f = x -> x  # 恒等関数:微分係数は常に1
(anonymous function)

julia> df = diff(f)
(anonymous function)

julia> df(0)
1.0f0

julia> df(1000)
1.0f0

julia> g = x -> x^2 + 3x + 1   # g(x) = x^2 + 3x + 1
(anonymous function)

julia> dg = diff(g)  # 導関数は g'(x) = 2x + 3
(anonymous function)

julia> dg(0)
3.0f0

julia> dg(5)
13.0f0

julia> h = x -> sin(x)  # h(x) = sin(x)
(anonymous function)

julia> dh = diff(h)  # 導関数は h'(x) = cos(x)
(anonymous function)

julia> dh(0)   # cos(0) = 1
0.9999998f0

julia> dh(pi/2)    # cos(π / 2) = 0
0.0

なんとなく良さそうです。

【広告】