[ソースコード]↓「二次曲線・二次曲面の標準形の計算・グラフ描画プログラムソースコード(Python)_1」からの続き↓
[QuadraticSurfaceCalculation]
import numpy as np
import math
# 二次形式の行列のランク(階数)
def drank(dis):
return np.linalg.matrix_rank(dis)
# 二次曲面の判別式(行列式 ≠ 0：有心二次曲面・柱面(楕円面・双曲面・錐面)／行列式 = 0：無心二次曲面・柱面(放物面・柱面))
def ddet(dis):
return np.linalg.det(dis)
# 二次曲面の分類の行列のランク(階数)
def classrank(cla):
return np.linalg.matrix_rank(cla)
# 二次曲面の分類の行列の行列式
def classdet(cla):
return np.linalg.det(cla)
# 二次形式の行列の固有値・固有ベクトル
def quadraticeigen(dis):
return np.linalg.eig(dis)
# 二次式部分の座標変換(二次形式の行列の固有値)
def eigenvalue2(dis):
e = quadraticeigen(dis)
return [e[0][0], e[0][1], e[0][2]]
# 二次形式の行列の固有ベクトル
def eigenvector(dis):
e = quadraticeigen(dis)
return e[1]
# 二次行列の回転角度
def rotateangle(dis):
# 回転行列を表す行列：
# [[cosYcosZ, -cosYsinZ, sinY],
# [sinXsinYcosZ+cosXsinZ, -sinXsinYsinZ+cosXcosZ, -sinXcosY],
# [-cosXsinYcosZ+sinXsinZ, cosXsinYsinZ+sinXcosZ, cosXcosY] ]
e = eigenvector(dis)
siny = e[0][2]
ay = math.asin(siny)
cosx = e[2][2] / np.cos(ay)
cosz = e[0][0] / np.cos(ay)
ax = math.acos(cosx)
az = math.acos(cosz)
return [ax / np.pi, ay / np.pi, az / np.pi]
# 一次式部分の座標変換(一次式部分の行列と二次形式の行列の固有ベクトルとの積)
def eigenvalue1(cla, dis):
con1 = np.dot([cla[3][0]*2, cla[3][1]*2, cla[3][2]*2], eigenvector(dis))
return [con1[0], con1[1], con1[2]]
# 座標変換式の平方完成・標準形・平行移動量
def standardtype(cla, dis):
# 判別行列のランク・行列値
dr = drank(dis)
dd = ddet(dis)
cr = classrank(cla)
cd = classdet(cla)
# 座標変換
e2 = eigenvalue2(dis)
e1 = eigenvalue1(cla, dis)
cx2 = e2[0]
cy2 = e2[1]
cz2 = e2[2]
cx1 = e1[0]
cy1 = e1[1]
cz1 = e1[2]
d = cla[3][3]
# 二次曲面の分類(行列値・行列のランクにより分類)
# classrank - drank = 0：錐面
# classrank - drank = 1：有心二次曲面
# classrank - drank = 2：無心二次曲面
if dr == 1: # ddet = 0
# classrank = 1：平面(y^2 = 0)(classdet = 0)
# classrank = 2：平行2平面・虚の平行2平面(y^2 = 1)(classdet = 0)
# classrank = 3：放物柱面(y^2 = x)(classdet = 0)
if round(cx2, 8) != 0:
# 座標変換式の平方完成・標準形
sy2 = 0
sz2 = 0
qx1 = cx1/(2*cx2)
if cr == 1: # cx2*(x+qx1)^2 = 0
sx2 = cx2
qy1 = 0
qz1 = 0
elif cr == 2: # cx2*(x+qx1)^2 + qd = 0
qd = d - (cx2*((cx1/(2*cx2))**2))
sx2 = cx2 / (-qd)
qy1 = 0
qz1 = 0
elif cr == 3:
if round(cy1, 8) != 0: # y = sx2*(x+qx1)^2 + qy1
sx2 = -(cx2/cy1)
qy1 = -( (cx2/cy1) * ( ( cx1/(2*cx2) )**2) ) + (d/cy1)
qz1 = 0
elif round(cz1, 8) != 0: # z = sx2*(x+qx1)^2 + qz1
sx2 = -(cx2/cz1)
qy1 = 0
qz1 = -( (cx2/cz1) * ( ( cx1/(2*cx2) )**2) ) + (d/cz1)
if round(cy2, 8) != 0:
# 座標変換式の平方完成・標準形
sx2 = 0
sz2 = 0
qy1 = cy1/(2*cy2)
if cr == 1: # cy2*(y+qy1)^2 = 0
sy2 = cy2
qx1 = 0
qz1 = 0
elif cr == 2: # cy2*(y+qy1)^2 + qd = 0
qd = d - (cy2*((cy1/(2*cy2))**2))
sy2 = cy2 / (-qd)
qx1 = 0
qz1 = 0
elif cr == 3:
if round(cx1, 8) != 0: # x = sy2*(y+qy1)^2 + qx1
sy2 = -(cy2/cx1)
qx1 = -( (cy2/cx1) * ( ( cy1/(2*cy2) )**2) ) + (d/cx1)
qz1 = 0
elif round(cz1, 8) != 0: # z = sy2*(y+qy1)^2 + qz1
sy2 = -(cy2/cz1)
qx1 = 0
qz1 = -( (cy2/cz1) * ( ( cy1/(2*cy2) )**2) ) + (d/cz1)
if round(cz2, 8) != 0:
# 座標変換式の平方完成・標準形
sx2 = 0
sy2 = 0
qz1 = cz1/(2*cz2)
if cr == 1: # cz2*(z+qz1)^2 = 0
sz2 = cz2
qx1 = 0
qy1 = 0
elif cr == 2: # cz2*(z+qz1)^2 + qd = 0
qd = d - (cz2*((cz1/(2*cz2))**2))
sz2 = cz2 / (-qd)
qx1 = 0
qy1 = 0
elif cr == 3:
if round(cx1, 8) != 0: # x = sz2*(z+qz1)^2 + qx1
sy2 = -(cy2/cx1)
qx1 = -( (cz2/cx1) * ( ( cz1/(2*cz2) )**2) ) + (d/cx1)
qz1 = 0
elif round(cy1, 8) != 0: # y = sz2*(z+qz1)^2 + qy1
sy2 = -(cz2/cy1)
qx1 = 0
qz1 = -( (cz2/cy1) * ( ( cz1/(2*cz2) )**2) ) + (d/cy1)
elif dr == 2: # ddet = 0
# classrank = 2：交わる2平面・1直線(x^2 + y^2 = 0)(classdet = 0)
# classrank = 3：楕円柱面・虚楕円柱面・双曲柱面(x^2 + y^2 = 1)(classdet = 0)
# classrank = 4：(x^2 + y^2 = z)／classdet < 0：楕円放物面／classdet > 0：双曲放物面
if round(cx2, 8) == 0:
# 座標変換式の平方完成・標準形
sx2 = 0
qy1 = cy1/(2*cy2)
qz1 = cz1/(2*cz2)
if cr == 2: # cy2*(y+qy1)^2 + cz2*(z+qz1)^2 = 0
sy2 = cy2
sz2 = cz2
qx1 = 0
elif cr == 3: # cy2*(y+qy1)^2 + cz2*(z+qz1)^2 + qd = 0
qd = d - (cy2*((cy1/(2*cy2))**2)) - (cz2*((cz1/(2*cz2))**2))
sy2 = cy2 / (-qd)
sz2 = cz2 / (-qd)
qx1 = 0
elif cr == 4: # x = sy2*(y+qy1)^2 + sz2*(z+qz1)^2 + qx1
sy2 = -(cy2/cx1)
sz2 = -(cz2/cx1)
qx1 = -( (cy2/cx1)*( (cy1/(2*cy2) )**2) ) - ( (cz2/cx1)*( (cz1/(2*cz2) )**2) ) + (d/cx1)
elif round(cy2, 8) == 0:
# 座標変換式の平方完成・標準形
sy2 = 0
qx1 = cx1/(2*cx2)
qz1 = cz1/(2*cz2)
if cr == 2: # cx2*(x+qx1)^2 + cz2*(z+qz1)^2 = 0
sx2 = cx2
sz2 = cz2
qy1 = 0
elif cr == 3: # cx2*(x+qx1)^2 + cz2*(z+qz1)^2 + qd = 0
qd = d - (cx2*((cx1/(2*cx2))**2)) - (cz2*((cz1/(2*cz2))**2))
sx2 = cx2 / (-qd)
sz2 = cz2 / (-qd)
qy1 = 0
elif cr == 4: # y = sx2*(x+qx1)^2 + sz2*(z+qz1)^2 + qy1
sx2 = -(cx2/cy1)
sz2 = -(cz2/cy1)
qy1 = -( (cx2/cy1)*( (cx1/(2*cx2) )**2) ) - ( (cz2/cy1)*( (cz1/(2*cz2) )**2) ) + (d/cy1)
elif round(cz2, 8) == 0:
# 座標変換式の平方完成・標準形
sz2 = 0
qx1 = cx1/(2*cx2)
qy1 = cy1/(2*cy2)
if cr == 2: # cx2*(x+qx1)^2 + cy2*(y+qy1)^2 = 0
sx2 = cx2
sy2 = cy2
qz1 = 0
elif cr == 3: # cx2*(x+qx1)^2 + cy2*(y+qy1)^2 + qd = 0
qd = d - (cx2*((cx1/(2*cx2))**2)) - (cy2*((cy1/(2*cy2))**2))
sx2 = cx2 / (-qd)
sy2 = cy2 / (-qd)
qz1 = 0
elif cr == 4: # z = sx2*(x+qx1)^2 + sy2*(y+qy1)^2 + qz1
sx2 = -(cx2/cz1)
sy2 = -(cy2/cz1)
qz1 = -( (cx2/cz1)*( (cx1/(2*cx2) )**2) ) - ( (cy2/cz1)*( (cy1/(2*cy2) )**2) ) + (d/cz1)
elif dr == 3: # ddet ≠ 0
# 座標変換式の平方完成：cx2*(x+qx1)^2 + cy2*(y+qy1)^2 + cz2*(z+qz1)^2 + qd = 0
qx1 = cx1/(2*cx2)
qy1 = cy1/(2*cy2)
qz1 = cz1/(2*cz2)
qd = cd / dd
#qd = d - (cx2*((cx1/(2*cx2))**2)) - (cy2*((cy1/(2*cy2))**2)) - (cz2*((cz1/(2*cz2))**2))
if cr == 3: # 錐面・1点・虚の錐面(classdet = 0)
# 標準形：sx2*x^2 + sy2*y^2 + sz2*z^2 = 0
sx2 = cx2
sy2 = cy2
sz2 = cz2
elif cr == 4: # 楕円面・虚楕円面・二葉双曲面・一葉双曲面(classdet ≠ 0)
# 標準形：sx2*x^2 + sy2*y^2 + sz2*z^2 = 1
sx2 = cx2 / (-qd)
sy2 = cy2 / (-qd)
sz2 = cz2 / (-qd)
return [[sx2, sy2, sz2], [qx1, qy1, qz1]]
[QuadraticSurfaceDraw]
# IPython Notebookモジュールのインポート
import importnb
import numpy as np
import sympy
# グラフ描画モジュールのインポート
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 別ファイルをインポート(jupyterモジュールインポート用)
with __import__('importnb').Notebook():
import QuadraticSurfaceCalculation
# グラフ描画(a2*x^2 + b2*y^2 + c2*z^2 + ab*xy + bc*yz + ca*zx + a1*x + b1*y + c1*z + d = 0)
def graphdraw(varlist):
# 変数名を定義
varname = ['a2', 'b2', 'c2', 'ab', 'bc', 'ca', 'a1', 'b1', 'c1', 'd']
# 動的に変数を定義
for i in range(len(varlist)):
globals()[varname[i]] = varlist[i]
# 二次形式の行列
dismatrix = [[a2, ab/2, ca/2], [ab/2, b2, bc/2], [ca/2, bc/2, c2]]
# 二次曲面の分類の行列
classmatrix = [[a2, ab/2, ca/2, a1/2], [ab/2, b2, bc/2, b1/2], [ca/2, bc/2, c2, c1/2], [a1/2, b1/2, c1/2, d]]
# 判別行列のランク・行列値
dr = QuadraticSurfaceCalculation.drank(dismatrix)
dd = QuadraticSurfaceCalculation.ddet(dismatrix)
cr = QuadraticSurfaceCalculation.classrank(classmatrix)
cd = QuadraticSurfaceCalculation.classdet(classmatrix)
# 座標変換
e2 = QuadraticSurfaceCalculation.eigenvalue2(dismatrix)
e1 = QuadraticSurfaceCalculation.eigenvalue1(classmatrix, dismatrix)
cx2 = e2[0]
cy2 = e2[1]
cz2 = e2[2]
cx1 = e1[0]
cy1 = e1[1]
cz1 = e1[2]
# 座標変換式の平方完成・標準形・平行移動量
st = QuadraticSurfaceCalculation.standardtype(classmatrix, dismatrix)
sx2 = st[0][0]
sy2 = st[0][1]
sz2 = st[0][2]
qx1 = st[1][0]
qy1 = st[1][1]
qz1 = st[1][2]
# 二次行列の回転角度
ra = QuadraticSurfaceCalculation.rotateangle(dismatrix)
# 文字を定義
x = sympy.Symbol('x')
y = sympy.Symbol('y')
z = sympy.Symbol('z')
# 数式設定
tos = float(a2)*x**2 + float(b2)*y**2 + float(c2)*z**2 + float(ab)*x*y + float(bc)*y*z + float(ca)*z*x + float(a1)*x + float(b1)*y + float(c1)*z
tor = tos + float(d)
tsq = round(cx2,2)*x**2 + round(cy2,2)*y**2 + round(cz2,2)*z**2 + round(cx1,2)*x + round(cy1,2)*y + round(cz1,2)*z + round(d,2)
if cr == 1:
if dr == 1: # 平面
# a2*x^2 + b2*y^2 + c2*z^2 + ab*xy + bc*yz + ca*zx + a1*x + b1*y + c1*z + d を因数分解
st = sympy.factor(tor)
elif cr == 2:
if dr == 1: # 平行2平面・虚の平行2平面
# a2*x^2 + b2*y^2 + c2*z^2 + ab*xy + bc*yz + ca*zx + a1*x + b1*y + c1*z を(-d)で割って因数分解
st = sympy.factor(tos/(-d))
elif dr == 2: # 交わる2平面・1直線
# a2*x^2 + b2*y^2 + c2*z^2 + ab*xy + bc*yz + ca*zx + a1*x + b1*y + c1*z + d を因数分解
st = sympy.factor(tor)
else: # 柱面・放物面・錐面・楕円面・一葉双曲面・二葉双曲面
# 数式設定
tsx = round(sx2,2)*x**2
tsy = round(sy2,2)*y**2
tsz = round(sz2,2)*z**2
# 描画領域の分割行列指定
fig, ax = plt.subplots(1, 1, figsize=(15, 15), subplot_kw=dict(projection='3d'))
# データを作成
x = np.linspace(-np.pi, np.pi, 200)
y = np.linspace(-2*np.pi, 2*np.pi, 100)
t = np.linspace(-2, 2, 100)
if round(dd, 8) == 0:
if round(cd, 8) == 0: # 柱面
pass
else: # classdet < 0：楕円放物面／classdet > 0：双曲放物面
x2,y2 = np.meshgrid(x,y)
if round(cx2, 8) == 0:
X = sy2 * (x2**2) + sz2 * (y2**2)
Y = x2
Z = y2
elif round(cy2, 8) == 0:
X = x2
Y = sx2 * (x2**2) + sz2 * (y2**2)
Z = y2
elif round(cz2, 8) == 0:
X = x2
Y = y2
Z = sx2 * (x2**2) + sy2 * (y2**2)
else:
if cd < 0:
if cx2 > 0 and cy2 > 0 and cz2 > 0: # 楕円面
x2,y2 = np.meshgrid(x,y)
X = 0 + sx2 * np.cos(y2) * np.sin(x2)
Y = 0 + sy2 * np.sin(y2) * np.sin(x2)
Z = 0 + sz2 * np.cos(x2)
else: # 二葉双曲面
y2,t2 = np.meshgrid(y,t)
r = (((t2**2)-1))**(1/2)
if cx2 < 0 and cy2 < 0 and cz2 > 0:
X = 0 + sx2 * r * np.cos(y2)
Y = 0 + sy2 * r * np.sin(y2)
Z = 0 + sz2 * t2
elif cx2 < 0 and cy2 > 0 and cz2 < 0:
X = 0 + sx2 * r * np.sin(y2)
Y = 0 + sy2 * t2
Z = 0 + sz2 * r * np.cos(y2)
elif cx2 > 0 and cy2 < 0 and cz2 < 0:
X = 0 + sx2 * t2
Y = 0 + sy2 * r * np.cos(y2)
Z = 0 + sz2 * r * np.sin(y2)
elif round(cd, 8) == 0:
if (cx2 > 0 and cy2 > 0 and cz2 >0) or (cx2 < 0 and cy2 < 0 and cz2 < 0): # 点楕円面・虚錐面
X = 0
Y = 0
Z = 0
else: # 錐面
y2,t2 = np.meshgrid(y,t)
if (cx2 < 0 and cy2 < 0 and cz2 > 0) or (cx2 > 0 and cy2 > 0 and cz2 < 0):
X = 0 + sx2 * t2 * np.cos(y2)
Y = 0 + sy2 * t2 * np.sin(y2)
Z = 0 + sz2 * t2
elif (cx2 < 0 and cy2 > 0 and cz2 < 0) or (cx2 > 0 and cy2 < 0 and cz2 > 0):
X = 0 + sx2 * t2 * np.sin(y2)
Y = 0 + sy2 * t2
Z = 0 + sz2 * t2 * np.cos(y2)
elif (cx2 > 0 and cy2 < 0 and cz2 < 0) or (cx2 < 0 and cy2 > 0 and cz2 > 0):
X = 0 + sx2 * t2
Y = 0 + sy2 * t2 * np.cos(y2)
Z = 0 + sz2 * t2 * np.sin(y2)
elif cd > 0:
if cx2 > 0 and cy2 > 0 and cz2 > 0: # 虚楕円面
pass
else: # 一葉双曲面
y2,t2 = np.meshgrid(y,t)
r = (((t2**2)+1))**(1/2)
if cx2 > 0 and cy2 > 0 and cz2 < 0:
X = 0 + sx2 * r * np.cos(y2)
Y = 0 + sy2 * r * np.sin(y2)
Z = 0 + sz2 * t2
elif cx2 > 0 and cy2 < 0 and cz2 > 0:
X = 0 + sx2 * r * np.sin(y2)
Y = 0 + sy2 * t2
Z = 0 + sz2 * r * np.cos(y2)
elif cx2 < 0 and cy2 > 0 and cz2 > 0:
X = 0 + sx2 * t2
Y = 0 + sy2 * r * np.cos(y2)
Z = 0 + sz2 * r * np.sin(y2)
# 縦と横の比を同じにする
ax.set_aspect('equal')
# 軸のラベルの設定
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# 視点の設定
ax.view_init(elev=20, azim=20)
# タイトル
titleorigin = rf'${sympy.latex(tor)} = 0$'
if cr <= 2:
titlerotate = ''
titlesquare = ''
titlemove = ''
sl = sympy.latex(st)
if cr == 1:
if dr == 1: # 平面
titlestandard = rf'$→ {sl} = 0$'
elif cr == 2:
if dr == 1: # 平行2平面
titlestandard = rf'$→ {sl} = 1$'
elif dr == 2: # 1直線・交わる2平面
titlestandard = rf'$→ {sl} = 0$'
else:
titlerotate = rf'$→ rotateX = {round(ra[0],2)} π / rotateY = {round(ra[1],2)} π / rotateZ = {round(ra[2],2)} π$'
titlesquare = rf'$→ {sympy.latex(tsq)} = 0$'
titlemove = rf'$→ moveX = {round(qx1,2)} / moveY = {round(qy1,2)} / moveZ = {round(qz1,2)}$'
sl = sympy.latex(tsx + tsy + tsz)
if cr == 3:
if dr == 1: # 放物柱面
if round(cx1, 8) != 0:
titlestandard = rf'$→ x = {sl}$'
elif round(cy1, 8) != 0:
titlestandard = rf'$→ y = {sl}$'
elif round(cz1, 8) != 0:
titlestandard = rf'$→ z = {sl}$'
elif dr == 2: # 楕円柱面・虚楕円柱面・双曲柱面
titlestandard = rf'$→ {sl} = 1$'
elif dr == 3: # 錐面
titlestandard = rf'$→ {sl} = 0$'
elif cr == 4:
if dr == 2: # 楕円放物面・双曲放物面
if round(cx2, 8) == 0:
titlestandard = rf'$→ x = {sl}$'
elif round(cy2, 8) == 0:
titlestandard = rf'$→ y = {sl}$'
elif round(cz2, 8) == 0:
titlestandard = rf'$→ z = {sl}$'
elif dr == 3: # 楕円面・一葉双曲面・二葉双曲面
titlestandard = rf'$→ {sl} = 1$'
ax.set_title(titleorigin + '\n' + titlerotate + '\n' + titlesquare + '\n' + titlemove + '\n' + titlestandard,
color='black', loc='center', size=16)
if (round(dd, 8) == 0) and (round(cd, 8) == 0):
pass
else:
# 回転・平行移動した曲面を描画
ax.plot_surface(X, Y, Z, alpha=0.5, antialiased=True, shade=True, color='#cccccc')
# 軸を自動調整
ax.autoscale()
fig.show()