はじめに

筆者は数学もプログラミングも得意ではないため、間違っているところやもっといい感じにできるところがあれば助言お願いいたしますm(_ _)m

四次元立方体とは

四次元立方体（正八胞体）とは、我々のよく知る立方体を四次元に拡張したものです。
今回は一辺の長さが$1$であるものを考えます。

頂点は座標が$(x, y, z, u)$の各要素が$0$か$1$であることを考えれば、$2^4=16$個あることがわかります。
辺の数も数えておきましょう。$(x, y, z, u)$のうち、どれか3つを固定すると辺が一つ決まります。

例：$x=0, y=1, z=1$を固定すると、$(0, 1, 1, 0)$と$(0, 1, 1, 1)$を結ぶ辺が決まる。

例のように考えると、辺の数は${}_4C_3 \times 2^3=32$個あることがわかります。

今回やること

この記事では、

四次元の方向ベクトル$\mathbf{e}$を設定する。
$\mathbf{e}$に対し、三次元の直交補空間を定め、原点を揃えた三次元座標系$(X, Y, Z)$を設定する。
前項で設定した四次元立方体を、$-\mathbf{e}$の方向に平行移動させて、三次元座標系$(X, Y, Z)$への射影を考える。
matplotlibで描画し、アニメーションにしてみる

みたいなことをやります。

三次元の場合

話を簡単にするために、まずは三次元立方体について考えていきましょう。
例えば$-\mathbf{e}$が下図の黒矢印の方向だとすると、直交補空間は$\mathbf{e}$と垂直な平面なので、黄緑色の平面となる。

黄緑の平面（直交補空間）に適当な$(X, Y)$座標系をとると、下図のように三角形だったり六角形だったりというような断面になることがわかります。

では、それぞれの頂点を計算で求めることを考察してみましょう。

$|\mathbf{e}|=1$であることを仮定します。座標系$(x, y, z)$の正規直交基底を$\mathbf{e}_x=(1, 0, 0), \mathbf{e}_y=(0, 1, 0), \mathbf{e}_z=(0, 0, 1)$とします。
$\mathbf{e}$の$z$成分が非ゼロだとすると、$\mathbf{e}, \mathbf{e}_x, \mathbf{e}_y$は線型独立と言えます。
Gram-Schmidtの直交化などによって、$\mathbf{e}, \mathbf{e}_x, \mathbf{e}_y$から正規直交基底$\mathbf{e}, \mathbf{e}_X, \mathbf{e}_Y$を作ると、$\mathbf{e}_X, \mathbf{e}_Y$が$\mathbf{e}$に対する直交補空間の正規直交基底となっています。

立方体の辺と直交補空間の交点が頂点となるので、立方体が$-t\mathbf{e}$だけ平行移動したときに、辺上の点$\mathbf{r}$が直交補空間に含まれる条件を考えると、

(\mathbf{r}-t\mathbf{e})\cdot\mathbf{e}=0 \iff
\mathbf{r}\cdot\mathbf{e}=t

となります。
また、辺上の点$\mathbf{r}$は、うまく立方体の頂点座標$\mathbf{node}$(辺の端点のどちらか)を選べば、

\mathbf{r}=\mathbf{node}+a\mathbf{e}_i\ (0\leq a\leq 1,\ i:x, y, z のどれか)

と書けます。このとき、$\mathbf{node}$の$i$成分は$0$となる。$a$について解くと、

a = \frac{t - \mathbf{node}\cdot\mathbf{e}}{\mathbf{e}_i\cdot\mathbf{e}}

となり、それぞれの辺に対して$\mathbf{node}$をうまく選ぶ($\mathbf{r}(0,1,a)$に対しては$\mathbf{node}(0,1,0)$)ことによって、与えられた$t$から$a$を求めることができ、$0\leq a\leq 1$を満たせば、$\mathbf{r}$は立方体の辺と直交補空間の交点となります。

交点の$(X, Y)$座標は$X =(\mathbf{r}-t\mathbf{e})\cdot\mathbf{e}_X, Y =(\mathbf{r}-t\mathbf{e})\cdot\mathbf{e}_Y$によってそれぞれ求められます。

交点を全て求めた後は、それぞれの点が同一の面上にある場合に線分で結べば、知りたい図形の形が分かります。

四次元の場合

四次元の場合も三次元の場合と同じように考えられるはずです。

$\mathbf{e}_x=(1, 0, 0, 0), \mathbf{e}_y=(0, 1, 0, 0), \mathbf{e}_z=(0, 0, 1, 0), \mathbf{e}_u=(0, 0, 0, 1)$、$\mathbf{e}$の$u$成分が非ゼロ、$|\mathbf{e}|=1$であることを仮定します。

$\mathbf{e}$と$\mathbf{e}_{x, y, z}$からGram-Schmidtの直交化によって、正規直交基底$\mathbf{e}, \mathbf{e}_X, \mathbf{e}_Y, \mathbf{e}_Z$を作る。$\mathbf{e}_X, \mathbf{e}_Y, \mathbf{e}_Z$が$\mathbf{e}$に対する直交補空間を張る。
$a = (t - \mathbf{node}\cdot\mathbf{e})/(\mathbf{e}_i\cdot\mathbf{e})$に従って交点座標を求める。
交点同士が同じ面上にあれば、線分で結ぶ。$\mathbf{p}_1$と$\mathbf{p}_2$同じ面上にあるかどうかは$(\mathbf{p}_1-\mathbf{p}_2)\cdot\mathbf{e}_i=0\ (i=x, y, z, u)$となる$i$が二つ以上あるかどうかを確認すれば良い。

といった手順で考えていきます。

サンプルコード

無駄な定義とか、センスのない部分もありますが、このようなコードを書きました。

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.art3d as art3d
import matplotlib.animation as animation

print('方向ベクトル(x, y, z, u)を入力してください。(u>0。他の要素は非負。半角スペース区切り。)')
e = np.array(list(map(float, input().split())))
ex = np.array([1, 0, 0, 0])
ey = np.array([0, 1, 0, 0])
ez = np.array([0, 0, 1, 0])
eu = np.array([0, 0, 0, 1])
basis = [ex, ey, ez, eu]

# 頂点
class Node:
    def __init__(self, x, y, z, u):
        self.r = np.array([x, y, z, u])

# 辺
# node1→node2のような方向を持っていると考える。
class Edge:
    def __init__(self, node1, node2):
        self.node1 = node1
        self.node2 = node2
        self.dir = node2.r - node1.r

# 点p1とp2が同一面上にあるかチェックする関数
def is_on_same_face(p1, p2):
    c = 0
    for b in basis:
        if np.dot(b, p1-p2) == 0:
            c += 1
    if c >= 2:
        return True
    else:
        return False

def draw_line(ax, p1, p2, im):
    line = art3d.Line3D([p1[0], p2[0]],[p1[1], p2[1]],[p1[2], p2[2]], color='c')
    im.append(ax.add_line(line))

# 正規化
def normalize(v):
    return v / np.linalg.norm(v, ord=2)

# Gram-Schmidtの直交化
def gram_schmidt(e):
    eX = ex - np.dot(e, ex) * e
    eX = normalize(eX)
    eY = ey - np.dot(e, ey) * e - np.dot(eX, ey) * eX
    eY = normalize(eY)
    eZ = ez - np.dot(e, ez) * e - np.dot(eX, ez) * eX - np.dot(eY, ez) * eY
    eZ = normalize(eZ)
    return (eX, eY, eZ)

node_list = []
edge_list = []

for i in range(16):
    # 頂点を登録
    x = i//8
    y = (i - 8*x) // 4
    z = (i - 8*x - 4*y) // 2
    u = i - 8*x - 4*y - 2*z
    node = Node(x, y, z, u)
    node_list.append(node)
    # 辺を登録
    # node1が(*, *, *, 0)、node2が(*, *, *, 1)のようになるようにする。
    if x == 1:
        edge = Edge(node_list[i-8], node)
        edge_list.append(edge)
    if y == 1:
        edge = Edge(node_list[i-4], node)
        edge_list.append(edge)
    if z == 1:
        edge = Edge(node_list[i-2], node)
        edge_list.append(edge)
    if u == 1:
        edge = Edge(node_list[i-1], node)
        edge_list.append(edge)

e = normalize(e)
eX, eY, eZ = gram_schmidt(e)

fig = plt.figure(figsize=(10, 7),dpi=120)
ax = fig.add_subplot(111, projection='3d')
ax.set_box_aspect((1,1,1))

ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

ims = []

for t in np.linspace(0, 2, 200):
    intersection_list = []
    projected_intersection_list = []

    im = []

    for edge in edge_list:
        # 交点を計算し、登録
        if np.dot(e, edge.dir) == 0:
            continue
        a = (t - np.dot(e, edge.node1.r)) / np.dot(e, edge.dir)
        if 0 <= a <= 1:
            intersection_list.append(a * edge.dir + edge.node1.r)

    for p in intersection_list:
        # 交点を直交補空間上の座標に変換、登録
        x = np.dot(p-t*e, eX)
        y = np.dot(p-t*e, eY)
        z = np.dot(p-t*e, eZ)
        projected_intersection_list.append(np.array([x, y, z]))

    for i in range(len(intersection_list)):
        for j in range(i + 1, len(intersection_list)):
            # 2点が同一面上にある場合は線分で結ぶ
            if is_on_same_face(intersection_list[i], intersection_list[j]):
                draw_line(ax, projected_intersection_list[i], projected_intersection_list[j], im)
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims, interval=10)
plt.show()

結果

入力を$(1, 1, 1, 1)$とした場合のアニメーションはこのようになりました。

正四面体→正八面体→正四面体のように変化しています。三次元の場合は正三角形→正六角形→正三角形のように変化することと似ていますね。

皆さんも色々入力を変えて試してみてください。