Newton法による非線形方程式(f(x)=0)の解（多次元）

############################################
# 多次元ニュートン法
# （３点(0.5,1.0),(0.0,1.5),(0.5,2.0)
#   を通る円の中心と半径）
#      coded by Y.Suganuma
############################################

############################################
# Newton法による非線形方程式(f(x)=0)の解
#      x0 : 初期値
#      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
#      eps2 : 終了条件２（｜f(x(k))｜＜eps2）
#      max : 最大試行回数
#      w : f(x)の微分の逆行列を入れる (n x 2n) 
#      n : xの次数
#      x : 解
#      fn : f(x)とその微分を計算する関数名
#      return : 実際の試行回数
#            （負の時は解を得ることができなかった）
#      coded by Y.Suganuma
############################################

def newton(x0, eps1, eps2, max, w, n, x, &fn) 

	ind = 0
	sw  = 0
	g   = Array.new()
	x1  = Array.new()
	for i1 in 0 ... n 
		x1[i1] = x0[i1]
		x[i1]  = x1[i1]
	end

	while sw == 0 && ind >= 0 

		sw   = 1
		ind += 1
		fn.call(0, x1, g, w, eps2)

		sw1 = 0
		for i1 in 0 ... n 
			if g[i1].abs() > eps2 
				sw1 = 1
				break
			end
		end
		if sw1 > 0 
			if ind <= max 
				sw1 = fn.call(1, x1, g, w, eps2)
				if sw1 == 0 
					for i1 in 0 ... n 
						x[i1] = x1[i1]
						for i2 in 0 ... n 
							x[i1] -= w[i1][i2+n] * g[i2]
						end
					end
					sw1 = 0
					for i1 in 0 ... n 
						if (x[i1]-x1[i1]).abs() > eps1 && (x[i1]-x1[i1]).abs() > eps1*x[i1].abs() 
							sw1 = 1
							break
						end
					end
					if sw1 > 0 
						for i1 in 0 ... n 
							x1[i1] = x[i1]
						end
						sw = 0
					end
				else 
					ind = -1
				end
			else 
				ind = -1
			end
		end
	end

	return ind
end

############################################
# 線形連立方程式を解く（逆行列を求める）
#      w : 方程式の左辺及び右辺
#      n : 方程式の数
#      m : 方程式の右辺の列の数
#      eps : 逆行列の存在を判定する規準
#      return : =0 : 正常
#               =1 : 逆行列が存在しない
#      coded by Y.Suganuma
############################################

def gauss(w, n, m, eps) 

	nm  = n + m
	ind = 0

	for i1 in 0 ... n 

		y1 = 0.0
		m1 = i1 + 1
		m2 = 0
					# ピボット要素の選択
		for i2 in i1 ... n 
			y2 = w[i2][i1].abs()
			if y1 < y2 
				y1 = y2
				m2 = i2
			end
		end
					# 逆行列が存在しない
		if y1 < eps 
			ind = 1
			break
					# 逆行列が存在する
		else    # 行の入れ替え
			for i2 in i1 ... nm 
				y1        = w[i1][i2]
				w[i1][i2] = w[m2][i2]
				w[m2][i2] = y1
			end
						# 掃き出し操作
			y1 = 1.0 / w[i1][i1]

			for i2 in m1 ... nm 
				w[i1][i2] *= y1
			end

			for i2 in 0 ... n 
				if i2 != i1 
					for i3 in m1 ... nm 
						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
					end
				end
			end
		end
	end

	return ind
end

########################################
# 関数値（f(x)）とその微分の計算
#      x : 変数値
#      y : 関数値
#      w : 微分の逆行列（後半）
#      eps : 逆行列の存在を判定する規準
#      return : =0 : 正常
#               =1 : 逆行列が存在しない
########################################
snx = Proc.new { |sw, x, y, w, eps| 
	if sw == 0   # 関数値
		y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
		y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
		y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
	else   # 微分
		for i1 in 0 ... 3 
			for i2 in 0 ... 3 
				w[i1][i2+3] = 0.0
			end
			w[i1][i1+3] = 1.0
		end

		w[0][0] = -2.0 * (0.5 - x[0])
		w[0][1] = -2.0 * (1.0 - x[1])
		w[0][2] = -2.0 * x[2]
		w[1][0] = -2.0 * (0.0 - x[0])
		w[1][1] = -2.0 * (1.5 - x[1])
		w[1][2] = -2.0 * x[2]
		w[2][0] = -2.0 * (0.5 - x[0])
		w[2][1] = -2.0 * (2.0 - x[1])
		w[2][2] = -2.0 * x[2]
		sw = gauss(w, 3, 3, eps)
	end
}

			# データの設定
eps1 = 1.0e-7
eps2 = 1.0e-10
max  = 20
n    = 3
x    = [0, 0, 0]
x0   = [0, 0, 2]
w    = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]]
			# 実行と結果
ind  = newton(x0, eps1, eps2, max, w, n, x, &snx)

print("ind = ", ind, "\n")
print("x ", x[0], " ", x[1], " ", x[2], "\n")
snx.call(0, x, x0, w, eps1)
print("f ", x0[0], " ", x0[1], " ", x0[2], "\n")