数値微分(Juliaで)
Juliaで数値微分を書いてみます。
数値微分の基本
数値微分は、直感的には微小幅 について関数の値を計算して差分を取れば、それが勾配になります。しかし、数学的にもう少しちゃんと根拠が示されています。
まず、微分を求めたい1変数関数を 、点をとします。そして、微小幅 について、関数 を の周りでテイラー展開します。
ただし、は、 以上の項です。同様に、関数 をの周りでテイラー展開します。
ここで上2つの式を引き算すると、都合の良いことに の項が消えます。
以上の項を無視することによって、
という結果が得られました。これは中心差分の公式と呼ばれます。上の式変形からわかるように、の項まで考慮しているので(計算の途中で消えたけどw)、誤差のオーダーが (二次精度)となります。
Juliaで実装
上式をJuliaで実装してみます。Float32
をベースとして使って、 としてJuliaの組み込み定数であるeps(Float16)
を使ってみます。
コード自体は非常に単純です。
function diff(f) h = eps(Float16) return (x -> (f(x+h) - f(x-h))/(2h)) end
受け取った関数 の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
なんとなく良さそうです。