Shade3D 公式

38 Newton 法( 2/2 ) [ Shade Labo ]


#1

関連記事:
  • 37 Newton 法 ( 1/2 )

**修正**     **38 - 3** [ script 38 - 5 ] の関数 newton_s( ) の収束判定部分を修正 ( 16 / 10 / 30 )

## 38 - 1  非線形連立方程式の解法
非線形連立方程式を数値微分を用いた **Newton 法**で解くには、次のようにします。
**< Newton 法による方程式の解法 >**

          [ 式 38 - 1 ]


**< Newton 法による連立方程式の解法 >**
ここでは簡単のために 二元連立方程式 のケースで記述しますが、N元連立方程式も同様に拡張できます。

          [ 式 38 - 2 ]


     

     [ 式 38 - 3 ]


          [ 式 38 - 4 ]



[ script 38 - 1 ] は次の連立程式の解を、初期解 ( x = 3, y = 1.5 ) , ( x = 1.5, y = -1 ) を与えて数値微分の Newton 法で解を求める script です。

          [ 式 38 - 5 ]


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

     


[ scrip 38 - 1 ]
#  newton_2(fc as function list, coeff as list, x0 as float list, d0 as float list, threshold as float, maxN as int, report as bool) as float list
#
#	Newton 法により 関数渡しによって与える 二元連立方程式の実根解を求める
#
#	fc :		関数 list
#	coeff :		各関数 に与える係数を格納した list の list
#				係数がなければ None list を渡す
#
#				2つの関数値の有効桁が同程度になるよう、必要なら scaling を施しておく必要がある 
#
#	x0 :		初期解 list
#	d0 :		導関数値をもとめるための差分 list
#
#	threshold :	収束判定しきい値
#	maxN :		打ち切り計算回数
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton_2(fc, coeff, x0, d0, threshold, maxN, report) :
	
	#  check
	if (len(fc) != 2) or (len(coeff) != 2) or (len(x0) != 2) or  (len(d0) != 2):
		return None
		
		
	x = [x0[0], x0[1]]
	dx = [d0[0], d0[1]]
	f1 = fc[0](x, coeff[0])
	f2 = fc[1](x, coeff[1])
	

	if (f1 == 0) and (f2 == 0) :			#  初期解が正しい解
		return x0
		
	else :
		result = False
		for kk in range(maxN) :	
		
			for i in range(2) :
				if abs(dx[i]) > d0[i] :		
					dx[i] = d0[i]			#  この方が収束が早い
			
			#  F[]
			F = [-f1, -f2]
											
			#  J[[,], [,]]		jacobi 行列	
			J_00 = (f1 - fc[0]([x[0] - dx[0], x[1]], coeff[0]))/dx[0]
			J_01 = (f1 - fc[0]([x[0], x[1] - dx[1]], coeff[0]))/dx[1]
			J_10 = (f2 - fc[1]([x[0] - dx[0], x[1]], coeff[1]))/dx[0]
			J_11 = (f2 - fc[1]([x[0], x[1] - dx[1]], coeff[1]))/dx[1]
			J = [[J_00, J_01], [J_10, J_11]]
			
			#  [J][dx] = [F] なる二元連立方程式を解く。
			det_J = J[0][0]*J[1][1] - J[0][1]*J[1][0]
			if det_J == 0 :
				if kk == 0 :	
					return None									#  初回にここに入るのは、関数におかしな所がある
				else :
					if (abs(f1) <= 10*threshold) and (abs(f2) <= 10*threshold) :
						dx = [0, 0]								#  収束したと見なす	ただし精度が落ちる
						result = True							#  関数 fc[0], fc[1] の 変数 x[0], x[1] 近傍における勾配に著しい差がある場合に起こりうる
					else :
						return None								#  収束計算の途中で計算不能となった
			else :
				det_J = 1/det_J		
				dx = [det_J*(J[1][1]*F[0] -  J[0][1]*F[1]), det_J*(- J[1][0]*F[0] +  J[0][0]*F[1])]
				
			#  x[ ] の更新
			for i in range(2) :
				x[i] += dx[i]
		
		
			#  関数値の更新
			f1 = fc[0](x, coeff[0])
			f2 = fc[1](x, coeff[1])
	
			if report :
				print 'f1 = ', f1, '     x[0] = ', x[0], '     dx[0] = ', dx[0]
				print 'f2 = ', f2, '     x[1] = ', x[1], '     dx[1] = ', dx[1]
				
			#  収束判定
			if (abs(f1) <= threshold) and (abs(f2) <= threshold) :
				result = True
				break
			
			
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None
		
				
		

def fc1(w, coeff) :
	[x, y] = w
	return (x - 3)**2 + y**2 - 3
	
	
def fc2(w, coeff) :
	[x, y] = w
	return (x - 1)**2 - y - 1.5
	
	
	
	
x0 = 3
y0 = 1.5
dx0 = 1e-5
dy0 = 1e-5
coeff = [None, None]

print ''
[x, y] = newton_2([fc1, fc2], coeff, [x0, y0], [dx0, dy0], 1e-12, 30, True)
print 'x = ', x, '   y = ', y
print 'fc1 = ',fc1([x, y], coeff), '   fc2 = ', fc2([x, y], coeff)



x0 = 1.5
y0 = -1
dx0 = 1e-5
dy0 = 1e-5
coeff = [None, None]

print ''
[x, y] = newton_2([fc1, fc2], coeff, [x0, y0], [dx0, dy0], 1e-12, 30, True)
print 'x = ', x, '   y = ', y
print 'fc1 = ',fc1([x, y], coeff), '   fc2 = ', fc2([x, y], coeff)

#2

38 - 2  例題 - デルタ12面体を作る ( その1 )


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

デルタ12面体は 上下2つのパーツの組合せで構成されています。
          [ fig 38 - 2 ]

下半分のパーツについて、2つの回転角度 **θ**, **φ** を変数として、newton_2( ) を用いて解きます。

     

     [ fig 38 - 3 ]


**θ**, **φ** を求めるために 関数 **newton_2( )** に引数として渡す関数 **func_1**, **func_2** と newton_2( ) を呼び出す main 文は次のようになります。
[ scrip 38 - 2 ]
def func(w, r) :			#  r :	辺長
	[theta, phi] = w
	
	p1 = vec3(-r/2, 0, 0)
	p2 = vec3(r/2, 0, 0)
	
	p3 = vec3(0, 0, -r*3**0.5/2)
	Mr1 = matrix().rotate_x(theta)
	p3 = p3*Mr1
	
	axis = p3 - p2
	Mt1 = matrix().translate(-p2[0], -p2[1], -p2[2])
	Mt2 = matrix().translate(p2[0], p2[1], p2[2])
	Q = quaternion().rotation(axis, phi)
	Mq = Q.matrix()
	p4 = p1*(Mt1*Mq*Mt2)
	
	return [p3, p4]
	
		

def func_1(w, r) :
	[p3, p4] = func(w, r)
	return p4[2]
	
	
	
def func_2(w, r) :
	[p3, p4] = func(w, r)
	return p4[0] + p3[2]
	
	
	

from math import pi	


r = 1000.			#  辺長

theta = pi/4
phi = pi*2./3
dt = 1e-5
dp = 1e-5
coeff = [r, r]

print ''
data = newton_2([func_1, func_2], coeff, [theta, phi], [dt, dp], 1e-12, 30, False)

if data == None :
	print 'error'
else :
	[theta, phi] = data
	print 'theta = ', theta*180/pi, '   phi = ', phi*180/pi


以上から、デルタ12面体を作成する script を [ script 38 - 3 ] に記します。
[ script 38 - 3 ]
from vec3 import *
from matrix import *
from quaternion import *



#  newton_2(fc as function list, coeff as list, x0 as float list, d0 as float list, threshold as float, maxN as int, report as bool) as float list
#
#	Newton 法により 関数渡しによって与える 二元連立方程式の実根解を求める
#
#	fc :		関数 list
#	coeff :		各関数 に与える係数を格納した list の list
#				係数がなければ None list を渡す
#
#				2つの関数値の有効桁が同程度になるよう、必要なら scaling を施しておく必要がある 
#
#	x0 :		初期解 list
#	d0 :		初期導関数値をもとめるための差分 list
#
#	threshold :	収束判定しきい値
#	maxN :		打ち切り計算回数
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton_2(fc, coeff, x0, d0, threshold, maxN, report) :
	
	#  check
	if (len(fc) != 2) or (len(coeff) != 2) or (len(x0) != 2) or  (len(d0) != 2):
		return None

			
	x = [x0[0], x0[1]]
	dx = [d0[0], d0[1]]
	f1 = fc[0](x, coeff[0])
	f2 = fc[1](x, coeff[1])
	

	if (f1 == 0) and (f2 == 0) :			#  初期解が正しい解
		return x0
		
	else :
		result = False
		for kk in range(maxN) :	
		
			for i in range(2) :
				if abs(dx[i]) > d0[i] :		
					dx[i] = d0[i]			#  この方が収束が早い
					
			#  F[]
			F = [-f1, -f2]
											
			#  J[[,], [,]]		jacobi 行列
			J_00 = (f1 - fc[0]([x[0] - dx[0], x[1]], coeff[0]))/dx[0]
			J_01 = (f1 - fc[0]([x[0], x[1] - dx[1]], coeff[0]))/dx[1]
			J_10 = (f2 - fc[1]([x[0] - dx[0], x[1]], coeff[1]))/dx[0]
			J_11 = (f2 - fc[1]([x[0], x[1] - dx[1]], coeff[1]))/dx[1]
			J = [[J_00, J_01], [J_10, J_11]]
			
			#  [J][dx] = [F] なる二元連立方程式を解く。
			det_J = J[0][0]*J[1][1] - J[0][1]*J[1][0]
			if det_J == 0 :
				if kk == 0 :	
					return None									#  初回にここに入るのは、関数におかしな所がある
				else :
					if (abs(f1) <= 10*threshold) and (abs(f2) <= 10*threshold) :
						dx = [0, 0]								#  収束したと見なす	ただし精度が落ちる
						result = True							#  関数 fc[0], fc[1] の 変数 x[0], x[1] 近傍における勾配に著しい差がある場合に起こりうる
					else :
						return None								#  収束計算の途中で計算不能となった
			else :
				det_J = 1/det_J		
				dx = [det_J*(J[1][1]*F[0] -  J[0][1]*F[1]), det_J*(- J[1][0]*F[0] +  J[0][0]*F[1])]
				
			#  x[ ] の更新
			for i in range(2) :
				x[i] += dx[i]
		
		
			#  関数値の更新
			f1 = fc[0](x, coeff[0])
			f2 = fc[1](x, coeff[1])
	
			if report :
				print 'f1 = ', f1, '     x[0] = ', x[0], '     dx[0] = ', dx[0]
				print 'f2 = ', f2, '     x[1] = ', x[1], '     dx[1] = ', dx[1]
				
			#  収束判定
			if (abs(f1) <= threshold) and (abs(f2) <= threshold) :
				result = True

			if result :
				break
			
			
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None
		
				

				
def func(w, r) :			#  r :	辺長
	[theta, phi] = w
	
	p1 = vec3(-r/2, 0, 0)
	p2 = vec3(r/2, 0, 0)
	
	p3 = vec3(0, 0, -r*3**0.5/2)
	Mr1 = matrix().rotate_x(theta)
	p3 = p3*Mr1
	
	axis = p3 - p2
	Mt1 = matrix().translate(-p2[0], -p2[1], -p2[2])
	Mt2 = matrix().translate(p2[0], p2[1], p2[2])
	Q = quaternion().rotation(axis, phi)
	Mq = Q.matrix()
	p4 = p1*(Mt1*Mq*Mt2)
	
	return [p3, p4]
	
		

def func_1(w, r) :
	[p3, p4] = func(w, r)
	return p4[2]
	
	
	
def func_2(w, r) :
	[p3, p4] = func(w, r)
	return p4[0] + p3[2]
	
	
	

from math import pi	
scene = xshade.scene()

r = 1000*scene.user_to_native_unit		#  辺長

theta = pi/4
phi = pi*2./3
dt = 1e-5
dp = 1e-5
coeff = [r, r]

data = newton_2([func_1, func_2], coeff, [theta, phi], [dt, dp], 1e-12, 30, False)
[theta, phi] = data



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

#  pLs :	3 x 3角面
p1 = vec3(-r/2, 0, 0)
p2 = vec3(r/2, 0, 0)
[p3, p4] = func([theta, phi], r)
p5 = vec3(p4)
p5[0] *= -1

pLs = []
pLs.append([p3, p2, p1])
pLs.append([p2, p3, p4])
pLs.append([p3, p1, p5])


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

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

Mt = matrix().translate(0, -(p3[1] + p4[1])/2, 0)
ob1.move_object(Mt)


#  ob1 を回転複製
Mr = matrix().rotate_y(pi)
ob1.copy_object(Mr)
ob2 = ob1.bro

Mr1 = matrix().rotate_z(pi)
Mr2 = matrix().rotate_y(pi/2)
M = Mr1*Mr2

ob1.copy_object(M)
ob2.copy_object(M)


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

#3

38 - 3  例題 - デルタ12面体を作る ( その2 )


**修正**     [ script 38 - 5 ] の関数 newton_s( ) の収束判定部分を修正 ( 16 / 10 / 30 )

12 個の正三角形からなるデルタ12面体を N元連立方程式 を解く関数 **newton_s( )** を用いて作ってみます。

ここでは 4変数をとり、4元連立方程式として解きます。


newton_2( )newton_s( ) との実質的な違いは、後者は連立方程式 [ J ] [ dx ] = [ F ] を解くのに Gauss の消去法( 関数 gausian_elimination( ) ) を用いている部分のみです。

     

     [ fig 38 - 4 ]


**p3(y)**, **P3(z)**, **P4(x)**, **P4(y)** を求めるために 関数 **newton_s( )** に引数として渡す関数 **func_1**, **func_2**, **func_3**, **func_4** と newton_s( ) を呼び出す main 文は次のようになります。
[ script 38 - 4 ]
def func(w, r) :			#  r :	辺長
	[y3, z3, x4, y4] = w
	
	p1 = vec3(-r/2, 0, 0)
	p2 = vec3(r/2, 0, 0)
	p3 = vec3(0, y3, z3)
	p4 = vec3(x4, y4, 0)
	
	return [p1, p2, p3, p4]
	
		

def func_1(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p2 - p3
	return v.abs() - r
	
	
def func_2(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p3 - p4
	return v.abs() - r
	
	
def func_3(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p4 - p2
	return v.abs() - r
	
			
def func_4(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	return p4[0] + p3[2]
	
	
	

from math import pi
scene = xshade.scene()

r = 1000*scene.user_to_native_unit		#  辺長

fc = [func_1, func_2, func_3, func_4]
coeff = [r, r, r, r]
y3 = 0.75*r
z3 = -0.75*r
x4 = 0.6*r
y4 = 0.9*r
x0 = [y3, z3, x4, y4]
dx0 = [1e-5, 1e-5, 1e-5, 1e-5]

data = newton_s(fc, coeff, x0, dx0, 1e-12, 30, False)

if data == None :
	print 'error'
else :
	print data


以上から、デルタ12面体を作成する script を [ script 38 - 5 ] に記します。
[ script 38 - 5 ]
**修正**     関数 newton_s( ) の収束判定部分を修正 ( 16 / 10 / 30 )
from vec3 import *
from matrix import *
from quaternion import *



#  gaussian_elimination(A as float list, B_ as float list, pivot_check as bool = True) as float list
#
#	連立一次方程式 [A][x] = [B_] の解を ガウスの消去法 により求める
#
#	解が求められないなら、None が返される
#
#	関数 newton_s( ) から gaussian_elimination() が呼ばれる時、2回目以降の収束計算内で pivot サイズが 0 になる場合がある
#	これは渡された連立方程式の関数の導関数の勾配に、当該の変数値近傍で著しい差がある場合に生じる
#	この場合、pivot_check = False としておけば True が返され、newton_s( ) では状況に応じて収束したと見なすこともできるが、精度はやや落ちる

def gaussian_elimination(A, B_, pivot_check = True) :

	n = len(B_)
	B = [ b for b in B_]
	w = [0. for i in range(n)]
	nk = n
	
	for k in range(nk) :
		#   ピボット選択
		p = k
		pMax = abs(A[k][k])
		for i in range(k + 1, n) :
			if abs(A[i][k]) > pMax :
				p = i
				pMax = abs(A[i][k])
					
		if p != k :				#  ピボット処理
			for i in range(k, n) :
				A[k][i], A[p][i] = A[p][i], A[k][i]
			B[k], B[p] = B[p], B[k]
				
		#  前進消去
		if A[k][k] == 0 :
			if pivot_check :
				return None
			else :
				return True		
		else :
			for i in range(k + 1, n) :
				w[i] = A[i][k]/A[k][k]
				for j in range(k + 1, n) :
					A[i][j] += - A[k][j]*w[i]
				B[i] += - B[k]*w[i]
					
	#  後退代入
	for i in range(n - 1, -1, -1) :
		for j in range(i + 1, n) :
			B[i] += - A[i][j]*B[j]
		B[i] /= A[i][i]
		
		
	return B
	
	
	

#  newton_s(fc as function list, coeff as list, x0 as float list, d0 as float list, threshold as float, maxN as int, report as bool) as float list
#
#	Newton 法により 関数渡しによって与える N元連立方程式の実根解を求める
#
#	fc :		関数 list
#	coeff :		各関数 に与える係数を格納した list の list
#				係数がなければ None list を渡す
#
#				2つの関数値の有効桁が同程度になるよう、必要なら scaling を施しておく必要がある 
#
#	x0 :		初期解 list
#	d0 :		初期導関数値をもとめるための差分 list
#
#	threshold :	収束判定しきい値 ( 関数値 f に対する判定値 )
#	maxN :		打ち切り計算回数
#
#	report :	収束計算の途中経過の report flag
#
#	与えられた関数の導関数は差分により求める
#	解が得られなかったなら、None を返す

def newton_s(fc, coeff, x0, d0, threshold, maxN, report) :
	
	#  check
	n = len(fc)
	if (len(coeff) != n) or (len(x0) != n) or (len(d0) != n) :
		return None


	x = [xx for xx in x0]
	dx = [ddx for ddx in d0]
#	x = [x0[0], x0[1]]					#  これは不可
#	dx = [d0[0], d0[1]]					#  これは不可

	f = []
	result = True
	for i in range(n) :
		ff = fc[i](x, coeff[i])
		f.append(ff)
		if ff != 0 :
			result = False
			
			
	if result :							#  初期解が正しい解
		return x0
		
	else :
		for kk in range(maxN) :	
		
			for i in range(2) :
				if abs(dx[i]) > d0[i] :		
					dx[i] = d0[i]		#  この方が収束が早い
					
			#  F[]
			F = [- ff for ff in f]
											
			#  J[[,], [,]]		jacobi 行列
			J = []
			for i in range(n) :
				JJ = []
				for j in range(n) :
					x1 = []
					for k in range(n) :
						if k == j :
							x1.append(x[k] - dx[k])
						else :
							x1.append(x[k])

					JJ.append((f[i] - fc[i](x1, coeff[i]))/dx[j])	
	
				J.append(JJ)
			
			
			#  Gaussの消去法により [J][dx] = [F] なる連立方程式を解く。
			dx = gaussian_elimination(J, F, kk == 0)
			if dx == None :
				return None				#  収束計算の初回にここに入るのは、関数におかしな所がある
			
			elif dx == True :			#  収束計算の2回目以降に pivot = 0 となった ( これは渡された連立方程式の関数の導関数の勾配に、当該の変数値近傍で著しい差がある場合に生じる )		
				if (abs(f1) <= 10*threshold) and (abs(f2) <= 10*threshold) :
					result = True		#  収束したと見なすが、精度が落ちる
				else :
					return None			#  収束計算の途中で計算不能となった
										
			else :	
				#  x[ ] の更新
				for i in range(n) :
					x[i] += dx[i]
		
		
				#  関数値の更新
				f = []
				for i in range(n) :
					ff = fc[i](x, coeff[i])
					f.append(ff)

				if report :
					for i in range(n) :
						print 'f[' + str(i) + '] = ', f[i], '     x[' + str(i) + '] = ', x[i], '     dx[' + str(i) + '] = ', dx[i]
				
			#  収束判定
			result = True					#  16/10/30 修正
			for ff in f :					#  16/10/30 修正
				if abs(ff) > threshold :	#  16/10/30 修正
					result = False			#  16/10/30 修正
					break					#  16/10/30 修正

			if result :
				break
			
			
		if report :
			if result :
				print '繰り返し計算回数 : ' + str(kk + 1)
			else :
				print '収束せず'
				
		if result :
			return x
		else :
			return None

				

				
def func(w, r) :			#  r :	辺長
	[y3, z3, x4, y4] = w
	
	p1 = vec3(-r/2, 0, 0)
	p2 = vec3(r/2, 0, 0)
	p3 = vec3(0, y3, z3)
	p4 = vec3(x4, y4, 0)
	
	return [p1, p2, p3, p4]
	
		

def func_1(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p2 - p3
	return v.abs() - r
	
	
def func_2(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p3 - p4
	return v.abs() - r
	
	
def func_3(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	v = p4 - p2
	return v.abs() - r
	
			
def func_4(w, r) :
	[p1, p2, p3, p4] = func(w, r)
	return p4[0] + p3[2]
	
	
	

from math import pi	
scene = xshade.scene()

r = 1000*scene.user_to_native_unit		#  辺長

fc = [func_1, func_2, func_3, func_4]
coeff = [r, r, r, r]
y3 = 0.75*r
z3 = -0.75*r
x4 = 0.6*r
y4 = 0.9*r
x0 = [y3, z3, x4, y4]
dx0 = [1e-5, 1e-5, 1e-5, 1e-5]

data = newton_s(fc, coeff, x0, dx0, 1e-12, 30, False)



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

#  pLs :	3 x 3角面
[p1, p2, p3, p4] = func(data, r)
p5 = vec3(p4)
p5[0] *= -1

pLs = []
pLs.append([p3, p2, p1])
pLs.append([p2, p3, p4])
pLs.append([p3, p1, p5])


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

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

Mt = matrix().translate(0, -(p3[1] + p4[1])/2, 0)
ob1.move_object(Mt)


#  ob1 を回転複製
Mr = matrix().rotate_y(pi)
ob1.copy_object(Mr)
ob2 = ob1.bro

Mr1 = matrix().rotate_z(pi)
Mr2 = matrix().rotate_y(pi/2)
M = Mr1*Mr2

ob1.copy_object(M)
ob2.copy_object(M)


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

#4

38 - 4  応用例


**newton( )** , **newton_2( )** , **newton_s( )** を利用すれば、5つの正多面体をベースにして **半正多面体** のようなものを **三角関数を一切使用せずに作れます**。
          [ fig 38 - 5 ]

#5

38 - 5  変数の scaling ( 重要 )


Newton 法を用いて多元連立方程式を解く場合、次のことに留意しなければなりません。
     **各方程式の変数のオーダーが大きく違わないこと**
例えば [ script 38 - 3 ] では変数 **θ**, **φ** について解を求め、解析中に扱われる θ, φ の値はおおよそ **-π ~ π** であり、この大きさが **変数のオーダー** です。

[ script 38 - 5 ] では4つの座標値を変数とし、そのオーダーはおおよそ 辺長 になります。


[ script 38 - 3 ], [ script 38 - 5 ] はともに各変数のオーダーが同程度になっていますから問題ありません。

しかし、例えば 角度 θ距離 L を変数にして解く場合、モデルによっては両者のオーダーが大きく異なるケースが考えられます。

この場合、Jacobi 行列の各要素のオーダーに著しい差が生じる可能性が大きくなり、連立方程式 [J] [dx] = [F] で得られる差分 [dx] の精度が低下し、収束速度の低下や、解が求まらないケースも考えられます。


          [ fig 38 - 5 ]


内部的には倍精度で計算が行われますから、ある程度オーダーに差があっても実質的な問題は生じません。

しかし収束判定条件は変数のオーダーの違いに関わりなく一律の値を基準にして行われますから、得られる解の精度にバラツキが生じます。

やはりお行儀よくしておくに越したことはなく、次のように handling する習慣をつけておくことを推奨します。


     予め変数のオーダーを揃えた形で方程式を組み立てる




radian : 無次元変数

  • 角度( radian )= 弧長 / 半径 = 長さ / 長さ = 無次元 ( - π ~ π )

**長さや距離** : 基準長で除して無次元化 ( 基準長としては bounding box size などを利用 )
  • 長さ = 長さ / 基準長 = 長さ / 長さ = 無次元 ( 大きくても数十のオーダーに収まるような基準長を用いる )

**bezier parameter** : 無次元変数 ( 0 ~ 1 )
例えば 38 - 4 で紹介した半正多面体を作る場合に、**角度 θ** と **距離 L** を変数とした連立方程式を組み立てるなら、基準長として **辺長** をとれば簡単です。

つまり 辺長 r の多面体 についての方程式ではなく、単位辺長の多面体 に対する方程式を組み立て、得られた距離の解に後から辺長を乗じるようにすればいいわけです。