74 岩手県川崎村門崎字宮畑 伊吹神社 明治17年(1884)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：3次元，球
#Julia #SymPy #算額 #和算 #数学
大球の中に小球 5 個，中球 2 個を容れる。小球は互いに外接しあい，すべてが中球とも外接し，大球に内接している。中球は大球に内接している（中球同士は外接していない）。小球の直径が与えられたとき，大球の直径はいかほどか。
大球の半径と中心座標を \(r_1, (0, 0, 0)\)
中球の半径と中心座標を \(r_2, (0, 0, z)\)
小球の半径と中心座標を \(r_3, (l, 0, 0); l = r_3/\sin(\pi/10)\)
とおき，以下の方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1, r2, r3, z
l = r3/sind(Sym(36))
eq1 = z^2 + l^2 - (r2 + r3)^2
eq2 = z + r2 - r1
eq3 = l + r3 - r1;
res = solve([eq1, eq2, eq3], (r1, r2, z))[1]
(r3 + 2*sqrt(2)*r3/sqrt(5 - sqrt(5)), r3*(-15*sqrt(2) - 10*sqrt(5 - sqrt(5)) + 2*sqrt(25 - 5*sqrt(5)) + 5*sqrt(10))/(5*(-3*sqrt(5 - sqrt(5)) - 3*sqrt(2) + sqrt(10) + sqrt(5)*sqrt(5 - sqrt(5)))), 2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30))
# r1
res[1] |> factor |> println
r3*(sqrt(5 - sqrt(5)) + 2*sqrt(2))/sqrt(5 - sqrt(5))
大球の半径 \(r_1\) は，小球の半径 \(r_3\) の \( \frac{\sqrt{5 - \sqrt{5}} + 2\sqrt{2}}{\sqrt{5 - \sqrt{5}}} = 2.70130161670408\) 倍である。
術では単に「3 倍」とあるが，理解不能である。
# r2
@syms d
apart(res[2]/r3, d) |> factor |> println
-(-30 - 10*sqrt(5) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 3*sqrt(10)*sqrt(5 - sqrt(5)))/20
# z
res[3] |> simplify |> println
2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30)
描画関数プログラムのソースを見る
function draw(r3, more=false)
pyplot(size=(500, 250), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, z) = (r3 + 2*sqrt(2)*r3/sqrt(5 - sqrt(5)), r3*(-15*sqrt(2) - 10*sqrt(5 - sqrt(5)) + 2*sqrt(25 - 5*sqrt(5)) + 5*sqrt(10))/(5*(-3*sqrt(5 - sqrt(5)) - 3*sqrt(2) + sqrt(10) + sqrt(5)*sqrt(5 - sqrt(5)))), 2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30))
l = r3/sind(Sym(36))
p1 = plot()
circle(0, 0, r1, :green)
circle(0, 0, r2, :blue)
rotate(l, 0, r3, angle=72)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
end
p2 = plot()
circle(0, 0, r1, :green)
circle22(0, z, r2, :blue)
circle(l, 0, r3)
circle(l*cosd(72), 0, r3)
circle(l*cosd(144), 0, r3)
plot(p1, p2)
end;
draw(1, true)
以下のアイコンをクリックして応援してください
