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

猫尾製作所

あまりアテにしないでね

数値微分の応用(Newton 法)

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

今日は数値微分の応用として、関数を定義するだけで(微分した関数を式の形で得られなくても、自動的に数値微分を行う)その根(解)を求めるプログラムを書きました。

ニュートン法とは次のような漸化式で、値域も定義域も実数(実数の部分集合でも可)である関数  f'(x) = a の根(解)を近似的に求めるメソッドのことです。

 x_{n+1} = x_n - \frac{f(x_n) - a}{f'(x_n)}
 x_0 は適当な解の予測値とします。そして、 \epsilon = \frac{f(x_n) - a}{f'(x_n)} の絶対値が適当な微小量より小さくなったら根の近似を打ち切るという点も重要です。あるいはある回数以上の近似計算を行っても、収束しない場合も近似計算を打ち切ります。

#! ruby -Ku
require "kconv"

# Newton 法 で方程式の近似解を求める

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

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

# x0 を初期値として Newton 法で f(x) = a の根を求める(a は任意の実数の定数)
def newton(x0, a)
	# 打ち切りの条件
	eps0 = 1.0e-12  # 打ち切り誤差の絶対値
	n0 = 30     # 計算回数

	x = x0
	n = 0
	while n < n0 do
		eps = (func(x) - a) / diff(x)
		x -= eps
		if eps.abs < eps0 then
			break
		end
		
		n += 1
	end
	
	if n < n0 then
		puts "#{n} 回の計算で x = #{x} に収束しました。"
	else
		puts "収束しませんでした。"
	end
end

# x0 を初期値として Newton 法で f(x) = 0 の根を求める
def newton0(x0)
	newton(x0, 0)
end

上のプログラムでは  f(x) = x^3 - x^2 + \cos x としています。 f(x) = a の右辺の定数  a は変更可能です。

実行例は次のようになります。

irb(main):001:0> require "newton.rb"
=> true
irb(main):002:0> newton0(0)
収束しませんでした。
=> nil
irb(main):003:0> newton0(-1)
5 回の計算で x = -0.680219358792778 に収束しました。
=> nil
irb(main):004:0> newton0(1)
8 回の計算で x = -0.680219358792778 に収束しました。
=> nil
irb(main):005:0> newton(-1, 1)
28 回の計算で x = -6.10574164002843e-09 に収束しました。
=> nil
irb(main):006:0> newton(-1, 2)
14 回の計算で x = 1.72416729384022 に収束しました。
=> nil
irb(main):007:0> newton(-1, 3)
6 回の計算で x = 1.91209402657632 に収束しました。
=> nil

関数を書き換えることで結構応用は広がります。