Newton法

　　元来，関数 Gauss は，終了状態を表す変数 ind を return していました．しかし，この方法では，ブロック内で関数 Gauss を呼んだとき，ブロックに引数として渡した配列の値の変化を，ブロックを呼んだ側で知ることができません．そのため，関数 Gauss によって目的とする配列を return するように修正し，その配列の 1 行目の最後の要素に ind を記憶しています．ブロック付き関数呼び出しと C++ などにおける関数名を引数とする場合との違いが大きく出てしまいました．しかし，このままの仕様では，ブロック付き関数呼び出しにおいて，複数の種類のデータを返したいような場合は不便だと思います．

#*****************************/
# Newton法による最小値の計算 */
#      coded by Y.Suganuma   */
#*****************************/

#******************************************************/
# 線形連立方程式を解く（逆行列を求める）              */
#      w : 方程式の左辺及び右辺                       */
#      n : 方程式の数                                 */
#      m : 方程式の右辺の列の数                       */
#      eps : 正則性を判定する規準                     */
#      return : =0 : 正常                             */
#               =1 : 逆行列が存在しない               */
#******************************************************/
def Gauss(w, n, m, eps)

	ind = 0
	nm  = n + m

	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

	w[0][nm] = ind
	return w
end

#******************************************************/
# 与えられた点xから，dx方向へk*dxだけ進んだ点における */
# 関数値，微係数，及び，Hesse行列の逆行列を計算する   */
#******************************************************/
					# 関数１
snx1 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = w[0] - 1.0
		y1   = w[1] - 2.0
		x1 * x1 + y1 * y1
	elsif sw == 1
		dx[0] = -2.0 * (x[0] - 1.0)
		dx[1] = -2.0 * (x[1] - 2.0)
	else
		h = [[2.0, 0.0, 1.0, 0.0, 0.0], [0.0, 2.0, 0.0, 1.0, 0.0]]
		Gauss(h, 2, 2, eps)
	end
}
					# 関数２
snx2 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = w[1] - w[0] * w[0]
		y1   = 1.0 - w[0]
		100.0 * x1 * x1 + y1 * y1
	elsif sw == 1
		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
		dx[1] = -200.0 * (x[1] - x[0] * x[0])
	else
		h[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
		h[0][1] = -400.0 * x[0]
		h[1][0] = h[0][1]
		h[1][1] = 200.0
		h[0][2] = 1.0
		h[0][3] = 0.0
		h[1][2] = 0.0
		h[1][3] = 1.0
		Gauss(h, 2, 2, eps)
	end
}
					# 関数３
snx3 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = 1.5 - w[0] * (1.0 - w[1])
		y1   = 2.25 - w[0] * (1.0 - w[1] * w[1])
		z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
		x1 * x1 + y1 * y1 + z1 * z1
	elsif sw == 1
		dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
		dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
	else
		x1      = 1.0 - x[1]
		x2      = 1.0 - x[1] * x[1]
		x3      = 1.0 - x[1] * x[1] * x[1]
		h[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
		h[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                  4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                  6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
		h[1][0] = h[0][1]
		h[1][1] = 2.0 * x[0] * x[0] +
                  4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                  12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                  18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
		h[0][2] = 1.0
		h[0][3] = 0.0
		h[1][2] = 0.0
		h[1][3] = 1.0
		Gauss(h, 2, 2, eps)
	end
}

#***************************************************************/
# 黄金分割法(与えられた方向での最小値)                         */
#      a,b : 初期区間 a < b                                    */
#      eps : 許容誤差                                          */
#      val : 関数値                                            */
#      ind : 計算状況                                          */
#              >= 0 : 正常終了(収束回数)                       */
#              = -1 : 収束せず                                 */
#      max : 最大試行回数                                      */
#      w : 位置                                                */
#      dw : 傾きの成分                                         */
#      h : Hesse行列の逆行列                                   */
#      snx : 関数値を計算する関数の名前                        */
#            （微分は，その符号を変える)                       */
#      return : 結果(w＋y＊dwのy)                              */
#***************************************************************/
def gold(a, b, eps, val, ind, max, w, dw, h, &snx)
					# 初期設定
	tau    = (Math.sqrt(5.0) - 1.0) / 2.0
	x      = 0.0
	count  = 0
	ind[0] = -1
	x1     = b - tau * (b - a)
	x2     = a + tau * (b - a)
	f1     = snx.call(0, x1, w, dw, h, eps)
	f2     = snx.call(0, x2, w, dw, h, eps)
					# 計算
	while count < max && ind[0] < 0
		count += 1
		if f2 > f1
			if (b-a).abs() < eps
				ind[0] = 0
				x      = x1
				val[0] = f1
			else
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = snx.call(0, x1, w, dw, h, eps)
			end
		else
			if (b-a).abs() < eps
				ind[0] = 0
				x      = x2
				val[0] = f2
				f1     = f2
			else
				a  = x1
				x1 = x2
				x2 = b - (1.0 - tau) * (b - a)
				f1 = f2
				f2 = snx.call(0, x2, w, dw, h, eps)
			end
		end
	end
					# 収束した場合の処理
	if ind[0] == 0
		ind[0] = count
		fa     = snx.call(0, a, w, dw, h, eps)
		fb     = snx.call(0, b, w, dw, h, eps)
		if fb < fa
			a  = b
			fa = fb
		end
		if fa < f1
			x      = a
			val[0] = fa
		end
	end

	return x
end

#*******************************************************/
# Newton法                                             */
#      opt_1 : =0 : １次元最適化を行わない             */
#              =1 : １次元最適化を行う                 */
#      max : 最大繰り返し回数                          */
#      n : 次元                                        */
#      eps : 収束判定条件                              */
#      step : きざみ幅                                 */
#      y : 最小値                                      */
#      x : x(初期値と答え)                             */
#      dx : 関数の微分値                               */
#      h : Hesse行列の逆行列                           */
#      snx : 関数値，その微分，及び，Hesse行列の逆行列 */
#      return : >=0 : 正常終了(収束回数)               */
#               =-1 : 収束せず                         */
#               =-2 : １次元最適化の区間が求まらない   */
#               =-3 : 黄金分割法が失敗                 */
#*******************************************************/
def Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx)

	sw1 = Array.new(n)
	y2  = Array.new(n)
	wk  = Array.new(n)
	for i1 in 0 ... n
		wk[i1] = 0.0
	end
	count = 0
	sw    = 0

	y1 = snx.call(0, 0.0, x, dx, h, eps)

	while count < max && sw == 0
					# 傾きの計算
		snx.call(1, 0.0, x, wk, h, eps)
					# Hesse行列の逆行列の計算
		h      = snx.call(2, 0.0, x, wk, h, eps)
		sw1[0] = h[0][2*n]   # もっと良い方法はあるか？
					# 収束していない
		if sw1[0] == 0
						# 方向の計算
			count += 1
			for i1 in 0 ... n
				dx[i1] = 0.0
				for i2 in 0 ... n
					dx[i1] += h[i1][n+i2] * wk[i2]
				end
			end
						# １次元最適化を行わない
			if opt_1 == 0
							# 新しい点
				for i1 in 0 ... n
					x[i1] += dx[i1]
				end
							# 新しい関数値
				y2[0] = snx.call(0, 0.0, x, dx, h, eps)
								# 関数値の変化が大きい
				if (y2[0]-y1).abs() > eps
					y1 = y2[0]
								# 収束（関数値の変化＜ｅｐｓ）
				else
					sw   = count
					y[0] = y2[0]
				end
						# １次元最適化を行う
			else
							# 区間を決める
				sw1[0] = 0
				f1     = y1
				sp     = step
				f2     = snx.call(0, sp, x, dx, h, eps)
				if f2 > f1
					sp = -step
				end
				for i1 in 0 ... max
					f2 = snx.call(0, sp, x, dx, h, eps)
					if f2 > f1
						sw1[0] = 1
						break
					else
						sp *= 2.0
						f1  = f2
					end
				end
							# 区間が求まらない
				if sw1[0] == 0
					sw = -2
							# 区間が求まった
				else
								# 黄金分割法
					k = gold(0.0, sp, eps, y2, sw1, max, x, dx, h, &snx)
								# 黄金分割法が失敗
					if sw1[0] < 0
						sw = -3
								# 黄金分割法が成功
					else
									# 新しい点
						for i1 in 0 ... n
							x[i1] += k * dx[i1]
						end
										# 関数値の変化が大きい
						if (y1-y2[0]).abs() > eps
							y1 = y2[0]
										# 収束（関数値の変化＜ｅｐｓ）
						else
							sw   = count
							y[0] = y2[0]
						end
					end
				end
			end
					# 収束（傾き＜ｅｐｓ）
		else
			sw = count
			y[0] = y1
		end
	end

	if sw == 0
		sw = -1
	end

	return sw
end
				# データの入力
str   = gets()
a     = str.split(" ")
fun   = Integer(a[1])
n     = Integer(a[3])
max   = Integer(a[5])
opt_1 = Integer(a[7])
str   = gets()
a     = str.split(" ")
eps   = Float(a[1])
step  = Float(a[3])
x     = Array.new(n)
dx    = Array.new(n)
y     = Array.new(1)
sw    = 0
str   = gets()
a     = str.split(" ")
h     = Array.new(n)
for i1 in 0 ... n
	x[i1]  = Float(a[i1+1])
	dx[i1] = 0.0
	h[i1]  = Array.new(2*n+1)
end
				# 実行
case fun
	when 1
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx1)
	when 2
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx2)
	when 3
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx3)
end
				# 結果の出力
if sw < 0
	printf("   収束しませんでした！")
	case sw
		when -1
			printf("（収束回数）\n")
		when -2
			printf("（１次元最適化の区間）\n")
		when -3
			printf("（黄金分割法）\n")
	end
else
	printf("   結果＝")
	for i1 in 0 ... n
		printf("%f ", x[i1])
	end
	printf(" 最小値＝%f  回数＝%d\n", y[0], sw)
end