はじめに
円錐を様々な角度で切断することを考える。この場合、切断面は円、楕円、放物線、双曲線になることが知られている。
これらの曲線を総称して2次曲線という。
そこで、今回はPythonを用いて円錐を描写し、様々な角度で切断するシミュレーションを実施する。
問題設定
今回用いる円錐の方程式を以下のように定義する。
z=r-\sqrt{x^2+y^2}
今回は、$r=1$とする。
この場合、円錐は以下のように図示する事ができる。
この円錐を3パターンの方法で切断することを考え、その際の断面をプロットするプログラムを作成した。
双曲線
縦に切断することを考える。
楕円
円錐に対して斜めに切断する。
放物線
円錐に対して母線と平行に切断する。
プログラム
python cone.py
"""
円錐(コーン)を色々な角度の平面で切断すると、切り口が
円・楕円・放物線・双曲線のいずれか(=2次曲線/円錐曲線)になる、
という古典的な事実を数値的に(点群として)確認・可視化するプログラム。
【円錐の作り方】
z = r - sqrt(x^2 + y^2)
は、頂点 (0, 0, r) から下に向かって開いていく「円錐面」の式で、
両辺を整理すると
x^2 + y^2 = (r - z)^2
という陰関数(円錐の標準形)になる。この円錐は
半径方向に 1 進むと z が 1 下がる
という形なので、軸(z軸)と母線(斜面)のなす半頂角は 45° になる。
この「45°」という値が、後述する各切断面が
楕円/放物線/双曲線のどれになるかを決める基準値として効いてくる。
【切断面が何になるかの判定(一般論)】
平面を z = A*x + B*y + C の形で表したとき(A, B は平面の傾き)、
円錐 x^2+y^2=(1-z)^2 との交線の式を x, y について整理すると
(1-A^2) x^2 - 2AB xy + (1-B^2) y^2 + (線形項) = 0
という2次曲線の一般形になり、判別式 B^2-4AC に相当する量が
4*(A^2 + B^2 - 1)
となる。つまり平面の傾きの大きさ A^2+B^2 を
円錐の傾き(45°、勾配の2乗で言うと 1)と比較すればよく、
A^2 + B^2 < 1 → 判別式 < 0 → 楕円(母線より緩い傾きで切る)
A^2 + B^2 = 1 → 判別式 = 0 → 放物線(母線と平行に切る)
A^2 + B^2 > 1 → 判別式 > 0 → 双曲線(母線より急な傾きで切る)
平面が軸(z軸)と平行/軸を含むような「垂直な平面」(x=定数 など)は
A→∞ に相当する極端なケースで、これも双曲線になる。
"""
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib # matplotlibで日本語(軸ラベル等)を文字化けさせずに表示するための設定
import math
from mpl_toolkits.mplot3d.axes3d import Axes3D # 3D描画(Axes3D)を使うために必要
r = 1 # 円錐の「高さ方向のスケール」。頂点は (0, 0, r) = (0, 0, 1)
n = 1000 # x, y方向の格子点数(解像度)。大きいほど滑らかだが計算が重くなる
# x, y をそれぞれ -1〜1 の範囲で n 分割し、円錐を描くための格子(メッシュ)を作る
x_ary = np.linspace(-1, 1, n)
y_ary = x_ary
X, Y = np.meshgrid(x_ary, y_ary)
#Z=r-(X**2+Y**2)**0.5
Z = np.zeros((n, n)) # 円錐の高さ z を格子点ごとに計算して入れていく配列(後のループで埋める)
# 各切断面ごとに、切り口(交線)付近の点を集めておくためのリスト
# 1: 平面 x = 一定 で切った断面 → 双曲線になるはず
x1_ary = []
y1_ary = []
z1_ary = []
# 2: 平面 x + y + 10z = 一定(ほぼ水平に近い平面)で切った断面 → 楕円になるはず
x2_ary = []
y2_ary = []
z2_ary = []
# 3: 平面 z = (y-x)/√2 - 定数 で切った断面 → 放物線になる
x3_ary = []
y3_ary = []
z3_ary = []
for i in range(n):
for k in range(n):
x = X[i][k]
y = Y[i][k]
z = r - (x**2 + y**2) ** 0.5 # 円錐面の高さ z = r - (原点からの距離)
Z[i][k] = z
# --- 切断面1: x ≈ 0.555 の「垂直な平面」(z軸と平行)で円錐を切る ---
# ちょうどの等式では格子点がほぼヒットしないので、薄い帯(0.55<x<0.56)で
# 「その平面付近を通る点」を代わりに拾い、切り口の形を点群として近似する
if 0.55 < x < 0.56:
y1_ary.append(y)
z1_ary.append(z)
# --- 切断面2: x + y + 10z ≈ 一定 という、ほとんど水平に近い平面で切る ---
# z の係数(10)が x, y の係数(1)よりずっと大きい = 傾き A,B が小さい
# → A^2+B^2 < 1 なので楕円になる
if 0.1 < y + x - 0.5 + 10 * z < 0.2:
x2_ary.append(x)
y2_ary.append(y)
z2_ary.append(z)
# --- 切断面3: (y-x)/√2 - z ≈ 一定、つまり z = (y-x)/√2 - 定数 という平面で切る ---
# 傾きは A=-1/√2, B=1/√2 なので A^2+B^2 = 1/2+1/2 = 1 となり、
# ちょうど円錐の母線(45°)と平行な平面 → 放物線になる
if 0.1 < (y - x) / math.sqrt(2) - z < 0.2:
x3_ary.append(x)
y3_ary.append(y)
z3_ary.append(z)
# --- 円錐そのものを3Dで描画 ---
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap='bwr', linewidth=0) # 高さzに応じて色が変わる曲面プロット
fig.colorbar(surf) # 高さと色の対応を示すカラーバー
plt.savefig("円錐.png")
plt.show()
# --- 切断面1(x≈一定)による切り口を (y, z) 平面上の散布図として描画 ---
# 前述の通り、この切り口は理論上「双曲線」になる
# ※直前が3D軸(ax)なので、新しくFigureを作ってから2D散布図を描く必要がある
plt.figure()
plt.scatter(y1_ary, z1_ary)
plt.savefig("双曲線.png")
plt.show()
plt.close()
# 図の作成
fig = plt.figure()
# 3次元の座標軸
ax2 = fig.add_subplot(projection="3d")
# --- 切断面2(x+y+10z≈一定)による切り口を3D散布図として描画 ---
# 傾きが緩やか(A^2+B^2<1)なので、理論上「楕円」になる
ax2.scatter(x2_ary, y2_ary, z2_ary)
plt.savefig("楕円.png")
plt.show()
# 図の作成
fig = plt.figure()
# 3次元の座標軸
ax3 = fig.add_subplot(projection="3d")
# --- 切断面3((y-x)/√2-z≈一定)による切り口を3D散布図として描画 ---
# A^2+B^2=1(母線と平行)なので、理論上「放物線」になる
ax3.scatter(x3_ary, y3_ary, z3_ary)
plt.savefig("放物線.png")
plt.show()
まとめ
今回は、円錐の切断面が2次曲線になることを視覚的に確かめた。
ただし、証明についてはかなり難しいので省略した。
数値計算ではこのように、証明が難しい命題の真偽を確認する際の最初のステップ（実験）としては、
有効である。
ただし、あくまで扱えるのは有限のパターンのみであることに注意したい。
参考文献