四面体の体積

筑波大の三谷さんのツィッターに四面体の体積を計算する話があった.

角封筒の上の辺の中央と下の辺の両端を結ぶ折目をつけ, それを山折りにすると, 合同な 二等辺三角形4枚による四面体が出来る. 元の封筒の横と縦が1と√2の時の 体積を計算せよという問題である. ちょっとやってみたところ, 存外簡単であった.

下の図で, 赤い線は元の封筒を展開したもの. 下の辺と上の辺は繋っている. この斜めの 線を谷折りにすると, ABCを底面, ACD'を手前の面, ABDを向うの面, 辺DEとD'E'は 繋っていて斜面だが上の面になる.

錐体の体積を計算するには, 底面積と高さが必要だが, 底面積は簡単だ. 辺BCが1, 辺OAは辺ABが√2だから, √7/2である.

高さについては, 三角形ABDは, 辺ABを折目として頂点Dが起き上がるから, Dに対応する 平面上の点はABと直交する破線で示す線に沿って下がり, AO上のGまで来ると下から上って来たD’と 出会って頂点になる. GAの距離はADの距離にtan αを掛けたものだが, tan αはOB/OAだしADは1だから, GA=1/√7である. よって 高さは√6/√7.

求める体積は (√7/2)*(√6)/((√7)*2*3)=1/2√6

三谷さんのツィッターにはヒントと称して, 四面体の6辺の長さからその体積を 求める式があった. 平面三角形のHeronの公式の三次元版である. そこにある式の 添字の付け方を, 私流に修正したのが次の図である.

このように添字を付けると, 四面体の体積Vの2乗は

の一番上の式で得られる. しかしこれも目がくらくらするので, 始めの3項と終りの4項を それぞれ6変数のtfと3変数のtgに分け, 2行目以下のように定義する. 最初のV²を計算する式を書いたLatexは次の様だ. 計算式と同様に分解してある.

\newcommand{\sq}[1]{l_#1^2}
\newcommand{\tf}[6]{\sq#1\sq#2(\sq#3+
\sq#4+\sq#5+\sq#6-\sq#1-\sq#2)}
\newcommand{\tg}[3]{-\sq#1\sq#2\sq#3}
$V^2=\frac{1}{144}[\tf012345+\tf234501\\
+\tf450123\tg025\tg034\tg124\tg135]$

これを使ってpythonで計算したのが次だ.

import math
def sq(x):
  return x*x
def tf(a,b,c,d,e,f):
  return sq(a)*sq(b)
    *(sq(c)+sq(d)+sq(e)+sq(f)-sq(a)-sq(b))
def tg(a,b,c):
  return -sq(a)*sq(b)*sq(c)
def v2(l0,l1,l2,l3,l4,l5):
  return (tf(l0,l1,l2,l3,l4,l5)+tf(l2,l3,l4,l5,l0,l1)
    +tf(l4,l5,l0,l1,l2,l3)
    +tg(l0,l2,l5)+tg(l0,l3,l4)+tg(l1,l2,l4)
    +tg(l1,l3,l5))/144
l0=1;l1=1;l2=math.sqrt(2);l3=math.sqrt(2)
l4=math.sqrt(2);l5=math.sqrt(2)
print(math.sqrt(v2(l0,l1,l2,l3,l4,l5)))

=> 0.20412414523193143

この値は確かに1/2√6 である.

ところで, 最初の封筒の四面体の体積は, 形状が特殊であったせいか, 簡単であった. 底面の3点の座標と, その各点から頂点までの距離が与えられている時にも体積は 得られるだろうか. 以下はその計算である. 式が沢山あるので, Latexでpdfに したもので示す. 底辺の3点, A, B, Cの座標をそれぞれ(x₀,y₀,0), (x₁,y₁,0),(x₂,y₂,0); 頂点Dの座標を(x,y,z)とする. A,B,Cから頂点までの距離をそれぞれ l₀,l₁,l₂とする.

今回もpythonで計算すると

x0=0;y0=-1/2;x1=math.sqrt(7)/2;y1=0;x2=0;y2=1/2
l0=math.sqrt(2);l1=1;l2=math.sqrt(2)
a0=2*(x0-x1); b0=2*(y0-y1)
c0=sq(l1)-sq(l0)+sq(x0)-sq(x1)+sq(y0)-sq(y1)
a1=2*(x1-x2); b1=2*(y1-y2)
c1=sq(l2)-sq(l1)+sq(x1)-sq(x2)+sq(y1)-sq(y2)

x=(c0*b1-c1*b0)/(a0*b1-a1*b0)
y=(a0*c1-a1*c0)/(a0*b1-a1*b0)

z=math.sqrt(sq(l0)-sq(x0-x)-sq(y0-y))
print(x,y,z)

#=> 0.9449111825230684 0.0 0.9258200997725514

期待通りの値が得られた. でもこんな実数の値より, 1/2√6 の方が嬉しい.