読者です 読者をやめる 読者になる 読者になる

猫尾製作所

あまりアテにしないでね

数値微分(プログラミング編)

Ruby プログラミング 数学 数値計算

iruby の環境があれば、Rubyインタプリタで実行することが出来ます。iruby は UNIX 環境上のみならず、Windows 環境上でも実行が可能です。また、スマートフォンでも実行可能である場合もあります。

Ruby で数値微分を行うプログラムを書いてみました。
 

# 関数 func : f(x) の x = a における微分係数 f'(a) を求める

# f(x) を定義する
def func(x)
	# 次の行の "y =" 以降を書き換えて関数の形を定義
	y = Math.log(x*x + 1) - Math.exp(Math.sin(x))
	
	return y
end

# 数値微分する(f'(x) の近似値を返す)
def diff(x)
	# 微小量を定義(あまり小さすぎてもダメ)
	dx = 1.0e-8
	
	# 右からの変化量 と 左からの変化量 の平均値
	df = func(x + dx) - func(x - dx)
	dfx = df / (dx*2)
	return (dfx)
end

#  f(x) = Math.log(x*x + 1) - Math.exp(Math.sin(x)) の厳密な導関数
# f'(x) = (2*x) / (x*x + 1) - Math.exp(Math.sin(x))*Math.cos(x)
def diff2(x)
	yp = (2*x) / (x*x + 1) - Math.exp(Math.sin(x))*Math.cos(x)
	return yp
end

 
上のコードの『関数の形を定義』と『微小量の定義』の部分を書き換えることで、任意の関数について任意のところでの微分係数を数値計算することができます。

iruby を読み込みます。

% irb
irb(main):001:0> require "diff.rb"
=> true

 
いくつかの点における  f(x) の近似値を出力してみます。

irb(main):002:0> func(0)
=> -1.0
irb(main):003:0> func(0.5)
=> -1.39200274512787
irb(main):004:0> func(1)
=> -1.62662964415591
irb(main):005:0> func(2)
=> -0.8731398155809
irb(main):006:0> func(-0.25)
=> -0.720200586428596
irb(main):007:0> func(-1.5)
=> 0.809852857038888

 
いくつかの点における数値微分の値を出力してみます。

irb(main):012:0> diff(0)
=> -0.999999993922529
irb(main):013:0> diff(0.5)
=> -0.617424211757367
irb(main):014:0> diff(-2)
=> -0.632373087228189
irb(main):015:0> diff(3.14)
=> 1.57988281346277
irb(main):016:0> diff(5.0)
=> 0.275886247180779

 
x を -5.0 から 5.0 まで 1.0 ごとに変化させて、数値微分の値を出力します。

irb(main):022:0> for i in (-5.0..5.0).step(1.0) do p diff(i) end
-1.1246583886404
0.922620446708322
0.259694743487415
-0.632373087228189
-1.23291132303471
-0.999999993922529
-0.253380783021839
1.83311685830745
1.74003854658267
0.777254394179749
0.275886247180779

 
導関数の厳密な値(の近似値)を求めます(今回だけ比較のために、わざわざ定義しておいたのです)

irb(main):023:0> for i in (-5.0..5.0).step(1.0) do p diff2(i) end
-1.12465840267138
0.922620454849071
0.259694725471414
-0.632373088725048
-1.23291133013811
-1.0
-0.253380767493447
1.83311686799583
1.74003856751332
0.777254412528328
0.275886251985847

 
導関数の近似値と厳密値(といえども計算機上である限り厳密ではないですが)の差を -5.0 から 5.0 まで 1.0 ごとに出力します。

1.40309812657335e-08
8.14074863075831e-09
1.80160016904907e-08
1.4968595252185e-09
7.10340075649185e-09
6.07747097092215e-09
1.55283923497507e-08
9.68838609338718e-09
2.09306478815563e-08
1.83485787674798e-08
4.80506789912738e-09

 
微小量を 1.0e-10 に(1.0e-8 から)変えてやって同じことを行います。
iruby を再起動して再び require します。

irb(main):001:0> require "diff.rb"
=> true
irb(main):002:0> for i in (-5.0..5.0).step(1.0) do p (diff(i)-diff2(i)).abs end
2.47872609593358e-06
6.96479023898355e-07
2.37335293973295e-07
4.59057805102248e-08
5.70634368646239e-08
8.27403709990904e-08
9.92524654019888e-07
1.23376140148679e-06
6.56305397139789e-07
5.95664551572561e-07
2.71258593809165e-07

 
かえって誤差が大きくなってしまいました。試行錯誤の結果ですと 1.0e-5 ぐらいがベストのようです。h = 1.0e-5 と定義して誤差を計算した結果を以下に示します。

irb(main):003:0> for i in (-5.0..5.0).step(1.0) do p (diff(i)-diff2(i)).abs end
8.65800764415781e-11
8.05295830019759e-11
8.18423107062927e-12
1.94155802546447e-12
9.0756291371008e-12
1.00008890058234e-12
3.69344554940199e-11
7.38549221779294e-11
1.40463196629526e-11
7.69451169446711e-12
8.90443274670361e-12