Shade3D 公式

37 Newton 法( 1/2 ) [ Shade Labo ]


#1

37 - 1  非線形方程式の解法


非線形方程式の数値解法として既に **08 Bezier Line 等分割** の **8 - 3 Newton 法** で Newton 法について紹介しました。
**< Newton 法 >**
関数 f(x) = 0 なる x の値を収束計算により求める
          [ 式 37 - 1 ]
導関数 f’(x) を差分を用いた数値微分で置き換えることもできます。

          [ 式 37 - 2 ]


計算手順は次のようになります。
  1. 最初の近似解 x と、最初の数値微分のための差分 dx の2つを与え、 f(x) と f’(x) を求める

  2. f(x) が十分に 0 に近ければ x が求める解であるとみなす

  3. そうでなければ、次の近似解を x + dx として f(x) を求め 2)の判定に戻る


[ script 37- 1 ] は次の 3次方程式の解を、初期解 x = 0, 初期差分 dx = 1e-5 を与えて数値微分の Newton 法で解を求める script です。

          [ 式 37 - 3 ]


結果は次のようになり、7回の trial で解が求まります。

     


[ scrip 37 - 1 ]
#  newton(fc as function, coeff as list, x0 as float, dx0 as float, threshold as float, maxN as int, report as bool) as float
#
#	Newton 法により 関数渡しによって与える関数の実根解を求める
#
#	fc :		関数
#	coeff :		関数 func( ) に与える係数を格納した list
#				係数がなければ None を渡す
#
#	x0 :		初期解
#	dx0 :		初期導関数値をもとめるための差分
#
#	threshold :	収束判定しきい値
#	maxN :		打ち切り計算回数	
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton(fc, coeff, x0, dx0, threshold, maxN, report) :

	#  check
	if threshold == None and epsilon == None :	
		return None
	
	x = x0
	dx = dx0
	f = fc(x, coeff)	
	ff = fc(x - dx, coeff)
	
	if f == 0 :								#  初期解が正しい解
		result = x0
		
	else :
		result = False
		for kk in range(maxN) :
			dx = - f*dx/(f - ff)
			ff = f
			x += dx
			f = fc(x, coeff)
			
			if report :
				print 'f = ', f, '     x = ', x, '     dx = ', dx
				
			#  収束判定
			if abs(f) <= threshold :
				result = True
				break
				
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None
						
					
	
	
def func(x, coeff) :
	[a, b, c, d] = coeff
	return a*x**3 + b*x**2 + c*x + d
	
	
	
coeff = [6, 30, 2, -50]	
x = 2
dx = 1e-5

print ''
x = newton(func, coeff, x, dx, 1e-12, 20, True)

if x != None :
	print 'x = ', x, '     f = ', func(x, coeff)	
	
			
			
					
	
	
def func(x, coeff) :
	[a, b, c, d] = coeff
	return a*x**3 + b*x**2 + c*x + d
	
	
	
coeff = [6, 30, 2, -50]	
x = 2
dx = 1e-5

print ''
x = newton(func, coeff, x, dx,1e-12, 20, True)

if x != None :
	print 'x = ', x, '     f = ', func(x, coeff)

#2

37 - 2  Newton 法の初期解


Newton 法で収束解を得るには **初期解が重要** で、これが予想される解にある程度近くなければ、収束しないか、あるいは別の解を拾ってしまいます。

また Newton 法では 解の存在判定はできません。( 解が求まらない ≠ 解が存在しない )

初期解の当たりをつけるには、関数 f(x) をグラフ出力して、f(x) = 0 を与える x のおおよその値の見つけます。


[ scrip 37 - 2 ] では、初期解の当たりをつける newton_check( ) で [ 式 37 - 3 ] をグラフ出力します。 ( x の範囲は -5 ~ 5 )


          [ fig 37 - 1 ]


[ scrip 37 - 2 ]
#  newton_check(xshade as xshade, f as function, coeff as float list, x1 as float, x2 as float, nx as int, size as float)
#
#	newton() より解を求める際に与える初期解を見つけるため、関数 func(x) の x = x1 ~ x2 の区間をグラフにして出力する
#	グラフサイズは X, Y 方向共に 1 なる大きさに規格化されたものを size 倍する
#	グラフの X 方向は x の範囲 x1 ~ x2 に関係なく、x >= 0 として出力される
#
#	func :		関数
#	coeff :		関数 func( ) に与える係数を格納した list
#				係数がなければ None を渡す
#	x1, x2		グラフ出力する区間
#	nx			グラフのプロット数
#	size		グラフサイズ

def newton_check(xshade, f, coeff, x1, x2, nx, size) :

	v = []
	min_y = f(x1, coeff)
	max_y = min_y
	dx = float(x2 - x1)/nx
	n = nx + 1
	
	for i in range(n) :
		x = x1 + i*dx
		y = f(x, coeff)
		min_y = min(min_y, y)
		max_y = max(max_y, y)
		v.append(vec3(x, y, 0))
		
		
	scene = xshade.scene()
	s = scene.user_to_native_unit
	ob = scene.active_shape()
	mx1 = matrix(ob.world_to_local_matrix)
	
	x = abs(x2 - x1)
	y = max_y - min_y
	Mt = matrix().translate(-x1, 0, 0)
	Ms = matrix().scale(s*size/x, s*size/y, 1)
	M = Mt*Ms*mx1
	M.transform(v)
	
	u = [vec3(0, 0, 0), vec3(size*s, 0, 0)]
	mx1.transform(u)
	

	scene.create_part('newton check')
	ob = scene.active_shape()
	scene.create_line(None, v, False)
	scene.create_line(None, u, False)

	ob.activate()
	
	


	
def func(x, coeff) :
	[a, b, c, d] = coeff
	return a*x**3 + b*x**2 + c*x + d
	

	
	
	

coeff = [6, 30, 2, -50]	

newton_check(xshade, func, coeff, -5., 5., 200, 2000)


[ fig 37 - 1 ] より、 [ 式 37 - 3 ] の3つの解のおおよその値が解ります。

それぞれの初期解を x = -5, -1, 2 として [ script 37 - 3 ] で解を求めると、次の 3つの解が得られます。

  • x = -4.51796380931
  • x = -1.62035795676
  • x = 1.13832176607

[ scrip 37 - 3 ]
#  newton(fc as function, coeff as list, x0 as float, dx0 as float, threshold as float, maxN as int, report as bool) as float
#
#	Newton 法により 関数渡しによって与える関数の実根解を求める
#
#	fc :		関数
#	coeff :		関数 func( ) に与える係数を格納した list
#				係数がなければ None を渡す
#
#	x0 :		初期解
#	dx0 :		初期導関数値をもとめるための差分
#
#	threshold :	収束判定しきい値 ( 関数値 f に対する判定値 )
#	maxN :		打ち切り計算回数	
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton(fc, coeff, x0, dx0, threshold, maxN, report) :

	#  check
	if threshold == None and epsilon == None :	
		return None
	
	x = x0
	dx = dx0
	f = fc(x, coeff)	
	ff = fc(x - dx, coeff)
	
	if f == 0 :								#  初期解が正しい解
		result = x0
		
	else :
		result = False
		for kk in range(maxN) :
			dx = - f*dx/(f - ff)
			ff = f
			x += dx
			f = fc(x, coeff)
			
			if report :
				print 'f = ', f, '     x = ', x, '     dx = ', dx
				
			#  収束判定
			if abs(f) <= threshold :
				result = True
				break
				
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None
			
					
	
	
def func(x, coeff) :
	[a, b, c, d] = coeff
	return a*x**3 + b*x**2 + c*x + d

	
	

coeff = [6, 30, 2, -50]	

for xx in [-5, -1, 2] :
	print ''
	x = newton(func, coeff, xx, 1e-5, 1e-12, 20, True)

	if x != None :
		print 'x = ', x, '     f = ', func(x, coeff)

#3

37 - 3  例題 - デルタ16面体を作る


**関数 newton( )** を用いて 16 個の正三角形からなるデルタ16面体を作ってみます。
          [ fig 37 - 2]

デルタ16面体は 中央部分の 8つの面からなる **正反四角柱** と上下2つの **正四角錐** の組合せで構成されます。
          [ fig 37 - 3 ]

これら2つのパーツを三角関数を用いないで Newton 法を用いて作ります。

正反四角柱 は2つの正方形を水平に上下 2H 離して置き、一方を 45度回転させたとき、点 P0, P1, P2 が正三角形になるような高さ H を求めます。


          [ fig 37 - 4 ]


**H** を求めるために 関数 **newton( )** に引数として渡す関数 **func_1** と newton( ) を呼び出す main 文は次のようになります。
[ scrip 36 - 3 ]
def func_1(H, coeff) :
	[r] = coeff					#  辺長
	p0 = vec3(-r/2, H, r/2)
	p2 = vec3(0, -H, r*2**0.5/2)
	v = p2 - p0
	
	return v.abs() - r


r = 1000.		#  辺長
coeff = [r]
print ''

#  H :	正反四角柱の高さ / 2
x, dx, threshold = r/2, r/1000, r*1e-12
H = newton(func_1, coeff, x, dx, threshold, 20, True)
print H


**正四角錐** では、点 **P0**, **P1**, **P2** が正三角形になるような高さ h を求めます。
          [ fig 37 - 5 ]
**h** を求めるために 関数 **newton( )** に引数として渡す関数 **func_2** と newton( ) を呼び出す main 文は次のようになります。
[ scrip 36 - 4 ]
def func_2(h, coeff) :
	[r] = coeff					#  辺長
	p0 = vec3(-r/2, 0, r/2)
	p2 = vec3(0, h, 0)
	v = p2 - p0
	
	return v.abs() - r


r = 1000.		#  辺長
coeff = [r]
print ''

#  h :	正四角錐の高さ
x, dx, threshold = r, r/1000, r*1e-12
h = newton(func_2, coeff, x, dx, threshold, 20, True)
print h


以上から、デルタ16面体を作成する script を [ script 37 - 5 ] に記します。
[ script 37 - 5 ]
#  newton(fc as function, coeff as list, x0 as float, dx0 as float, threshold as float, maxN as int, report as bool) as float
#
#	Newton 法により 関数渡しによって与える関数の実根解を求める
#
#	fc :		関数
#	coeff :		関数 func( ) に与える係数を格納した list
#				係数がなければ None を渡す
#
#	x0 :		初期解
#	dx0 :		初期導関数値をもとめるための差分
#
#	threshold :	収束判定しきい値 ( 関数値 f に対する判定値 )
#	maxN :		打ち切り計算回数	
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton(fc, coeff, x0, dx0, threshold, maxN, report) :

	#  check
	if threshold == None and epsilon == None :	
		return None
	
	x = x0
	dx = dx0
	f = fc(x, coeff)	
	ff = fc(x - dx, coeff)
	
	if f == 0 :								#  初期解が正しい解
		result = x0
		
	else :
		result = False
		for kk in range(maxN) :
			dx = - f*dx/(f - ff)
			ff = f
			x += dx
			f = fc(x, coeff)
			
			if report :
				print 'f = ', f, '     x = ', x, '     dx = ', dx
				
			#  収束判定
			if abs(f) <= threshold :
				result = True
				break
				
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None
			
			
					
	
	
def func_1(H, coeff) :
	[r] = coeff					#  辺長
	p0 = vec3(-r/2, H, r/2)
	p2 = vec3(0, -H, r*2**0.5/2)
	v = p2 - p0
	
	return v.abs() - r
	
	
	
def func_2(h, coeff) :
	[r] = coeff					#  辺長
	p0 = vec3(-r/2, 0, r/2)
	p2 = vec3(0, h, 0)
	v = p2 - p0
	
	return v.abs() - r
	

	

from matrix import *
from vec3 import *	
import labo

	
scene = xshade.scene()

r = 1000*scene.user_to_native_unit		#  辺長
coeff = [r]


#  H :	正反四角柱の高さ / 2
x, dx, threshold = r/2, r/1000, r*1e-12
H = newton(func_1, coeff, x, dx, threshold, 20, False)


#  h :	正四角錐の高さ
x, dx, threshold = r, r/1000, r*1e-12
h = newton(func_2, coeff, x, dx, threshold, 20, False)



#  面を構成する線形状を作成して polygon mesh に変換するが、
#  Shade の仕様により、線形状の面法線が一般ルールと逆向きの左回りで定義されている。
#
#  ポリゴン面としての線形状を左回りで作るのは、なんとも気持ち悪いので、
#  一般的なルール通り右回りで作成し、Shade に出力する直前に反転するものとする。

#  pLs1 :	正反四角柱 の面リスト
pL1 = [vec3(-r/2, H, r/2), vec3(r/2, H, r/2), vec3(r/2, H, -r/2), vec3(-r/2, H, -r/2)]
pL1.append(pL1[0])
a = 1/2**0.5
pL2 = [vec3(0, -H, r*a), vec3(r*a, -H, 0), vec3(0, -H, -r*a), vec3(-r*a, -H, 0)]
pL2.append(pL2[0])

pLs1 = []
for i in range(4) :
	pL = [pL1[i], pL2[i], pL1[i + 1]]
	pLs1.append(pL)
	pL = [pL2[i], pL2[i + 1], pL1[i + 1]]
	pLs1.append(pL)

	
#  pLs2 :	上側正四角錐 の面リスト
p0 = vec3(0, H + h, 0)
pLs2 = []
for i in range(4) :
	pL = [pL1[i], pL1[i + 1], p0]
	pLs2.append(pL)
	

#  pLs1, pL2 を線形状として出力し、ポリゴンメッシュ ob1, ob2 に変換
scene.create_part('デルタ16面体')
ob = scene.active_shape()
scene.cancel_transformation()

scene.create_part()
ob1 = scene.active_shape()
scene.begin_creating()
for pL in pLs1 :
	pL.reverse()						#  ここで反転しておく
	scene.create_line(None, pL, True)
scene.end_creating()
ob1.activate()
ob1.convert_to_polygon_mesh_with_subdivision_level(0)

scene.create_part()
ob2 = scene.active_shape()
scene.begin_creating()
for pL in pLs2 :
	pL.reverse()						#  ここで反転しておく
	scene.create_line(None, pL, True)
scene.end_creating()
ob2.activate()
ob2.convert_to_polygon_mesh_with_subdivision_level(0)
ob2 = scene.active_shape()


#  ob3 :		 下側正四角錐 polygon mesh
u1, v1 = pL1[0] - vec3(0, H, 0), pL1[1] - vec3(0, H, 0)
u2, v2 = pL2[1] - vec3(0, -H, 0), pL2[0] - vec3(0, -H, 0)
Q = labo.attitude_control_quaternion(u1, v1, u2, v2, 0)
Mq = Q.matrix()
ob3 = ob2.copy_object(Mq)


#  ob1, ob2, ob3 を一つの polygon mesh に統合
ob.activate()
ob.convert_to_polygon_mesh_with_subdivision_level(0)
ob = scene.active_shape()
ob.threshold = 0