数値微分の検討、微小hの値について
$${xの関数f(x)}$$の微分は
$$
f'(x) = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}
$$
と定義されていて、数値計算ではごく小さな$${h}$$を使って
$$
f'(x) \approx \frac{f(x+h) - f(x)}{h}
$$
で求まるとされています。
しかし、この$${h}$$はいくらならよいのかがよくわかりません。あるところには$${h}$$は計算機イプシロンの2倍が良いなどと書かれていましたが、$${x}$$の値によっては丸め誤差の影響を受けてしまうでしょう。そこで、簡単な例で$${h}$$の値と数値微分の結果を調べました。
例題
関数$${f(x)を f(x) =x^2}$$ とします。微分の真値は$${f'(x) = 2x}$$です(Fig.1)。
$${x = 1}$$ のとき、$${h}$$ を 0.1から始めて順に半分にしていき、
$${f'(x) = \frac{f(x+h) - f(x)}{h}}$$を求めました。結果をFig.2に示します。
$${f'(1)の値は2でなければなりませんが、hが10^{-13}}$$より小さくなるあたりから2からずれています。丸め誤差の影響が出ているようです。
また、$${hが10^{-1}近辺でも2からずれています。}$$拡大したものを、Fig.3に示します。
当然ですが$${h}$$の値が大きいと微分の値は正確には求まりません。Fig.3上では$${hが10^{-4}}$$程度より小さくなれば良さそうに見えます。
一般的な関数の数値微分を求める場合、関数に与える値によって$${h}$$の適正な値は変わります。そこで、$${h}$$を小さくしながら繰り返し$${f'}$$を求め、直前に求めた$${f'}$$との差異が必要な精度に達したときにその値を微分とすればよいと思います。簡単なコードを示します(Excel VBA)。ここでは精度として Accuracy = 0.000001 を与えています。
関数funcは$${x^2}$$のままですが、これは任意の関数に変更できます。条件分岐をして値が決まるものでもよいし、ひとつの式として表現できなくても何らかの計算で$${x}$$に対して値が一意に決まればそれでよいです(実際、そういう計算で使っています)。
ただし、関数に連続性がない点$${x}$$や、微分が不可能な点$${x}$$では、奇妙な結果が得られたり計算できない可能性があります。
'微分の計算
Sub TestDifferential()
Dim df As Double
df = differential(1)
Debug.Print df
End Sub
Function differential(x As Double)
Const Accuracy = 0.000001
Const iMax = 1000
Dim dif As Double
Dim difpost As Double
Dim h As Double
Dim i As Long
i = 0
h = 1
difpost = (func(x + h) - func(x)) / h
h = h / 2
Do
dif = (func(x + h) - func(x)) / h
'Debug.Print dif
If Abs(dif - difpost) < Accuracy Then
Exit Do
End If
h = h / 2
difpost = dif
i = i + 1
Loop While i < iMax
differential = dif
End Function
Function func(x As Double) As Double
func = x ^ 2
End Function
一点の微分を求めるために、繰り返し微分を計算しなくてはなりませんが、計算結果の精度が保証される(完全ではない)ので有効な数値微分の方法だと思います。
追記
不連続点、微分可能性については、
$${f'(x) = \frac{f(x+h) - f(x)}{h}}$$と$${f'(x) = \frac{f(x) - f(x-h)}{h}}$$を両方求めて比較して、一致するなら微分が可能といえますし、一致しないなら微分不可能と判定すればよいと思います。
応援してやろうということで、お気持ちをいただければ嬉しいです。もっと勉強したり、調べたりする糧にしたいと思います。