# -*- coding: UTF-8 -*-
from math import *
import numpy as np
############################################
# Newton法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# dfn : f(x)の微分を計算する関数名
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# w : f(x)の微分の逆行列を入れる (n x 2n)
# n : xの次数
# x : 解
# return : 実際の試行回数
# （負の時は解を得ることができなかった）
# coded by Y.Suganuma
############################################
def newton(fn, dfn, x0, eps1, eps2, max, w, n, x) :
ind = 0
sw = 0
g = np.empty(n, np.float)
x1 = np.empty(n, np.float)
for i1 in range(0, n) :
x1[i1] = x0[i1]
x[i1] = x1[i1]
while sw == 0 and ind >= 0 :
sw = 1
ind += 1
fn(x1, g)
sw1 = 0
for i1 in range(0, n) :
if abs(g[i1]) > eps2 :
sw1 = 1
break
if sw1 > 0 :
if ind <= max :
sw1 = dfn(x1, w, eps2)
if sw1 == 0 :
for i1 in range(0, n) :
x[i1] = x1[i1]
for i2 in range(0, n) :
x[i1] -= w[i1][i2+n] * g[i2]
sw1 = 0
for i1 in range(0, n) :
if abs(x[i1]-x1[i1]) > eps1 and abs(x[i1]-x1[i1]) > eps1*abs(x[i1]) :
sw1 = 1
break
if sw1 > 0 :
for i1 in range(0, n) :
x1[i1] = x[i1]
sw = 0
else :
ind = -1
else :
ind = -1
return ind
############################################
# 線形連立方程式を解く（逆行列を求める）
# 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 range(0, n) :
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in range(i1, n) :
y2 = abs(w[i2][i1])
if y1 < y2 :
y1 = y2
m2 = i2
# 逆行列が存在しない
if y1 < eps :
ind = 1
break
# 逆行列が存在する
else : # 行の入れ替え
for i2 in range(i1, nm) :
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in range(m1, nm) :
w[i1][i2] *= y1
for i2 in range(0, n) :
if i2 != i1 :
for i3 in range(m1, nm) :
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
return ind
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import newton, gauss
############################################
# 多次元ニュートン法
# （３点(0.5,1.0),(0.0,1.5),(0.5,2.0)
# を通る円の中心と半径）
# coded by Y.Suganuma
############################################
# 関数値（f(x)）の計算
def snx(x, y) :
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];
# 関数の微分の計算
# x : 変数値
# w : 微分の逆行列（後半）
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
def dsnx(x, w, eps) :
for i1 in range(0, 3) :
for i2 in range(0, 3) :
w[i1][i2+3] = 0.0
w[i1][i1+3] = 1.0
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)
return sw
# データの設定
eps1 = 1.0e-7
eps2 = 1.0e-10
max = 20
n = 3
x = np.array([0, 0, 0], np.float)
x0 = np.array([0, 0, 2], np.float)
w = np.empty((n, 2*n), np.float)
# 実行と結果
ind = newton(snx, dsnx, x0, eps1, eps2, max, w, n, x)
print("ind = " + str(ind))
print("x " + str(x[0]) + " " + str(x[1]) + " " + str(x[2]))
snx(x, x0)
print("f " + str(x0[0]) + " " + str(x0[1]) + " " + str(x0[2]))