三重県松阪市 意非多神社 文化11年(1814)
和算の館
http://wasan.jp/kyoto/isida2.html
キーワード：3次元，球7個，正四面体
#Julia #SymPy #算額 #和算 #数学
正四面体の中に，甲球 4 個，乙球 3 個を容れる。乙球の直径が 155 寸のとき，甲球の直径はいかほどか。
「算額（その1886）」も参照のこと。
甲球と乙球は，正四面体の頂点 D と辺 AD の中点 E を結ぶ直線に接する。
正四面体の底面を \(x-y\) 平面に置き，底面の重心を原点とする。
正四面体が内接する球の半径を \(R\)
正四面体の高さを \(\sqrt{2}R\)
甲球の半径と中心座標を \(r_1,\ (x_{11},\ 0,\ r_1),\ (x_{12},\ r_1,\ r_1),\ (0,\ 0,\ z_1)\)
乙球の半径と中心座標を \(r_2,\ (x_{21},\ 0,\ z_2),\ (x_{22},\ r_2,\ z_2)\)
とおき，以下の連立方程式の解を求める。
解析解を求める
まず，下にある 4 個の甲球について，算法助術の公式80 より，\(z_1 = r_1 (3 + 2\sqrt{6})/3\) がわかる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r2::positive, R::positive, r1::positive, x11::positive, x12::negative, z1::positive, x21::positive, x22::negative, z2::positive
z1 = r1*(3 + 2√Sym(6))/3 # 算法助術公式80
eq1 = (x12 - x11)^2 + r1^2 - 4r1^2
eq2 = x11^2 + (z1 - r1)^2 - 4r1^2
eq3 = x12^2 + r1^2 + (z1 - r1)^2 - 4r1^2
eq4 = (x22 - x21)^2 + r2^2 - 4r2^2
eq5 = x21^2 + (z2 - z1)^2 - (r1 + r2)^2
eq6 = x22^2 + r2^2 + (z2 - z1)^2 - (r1 + r2)^2
eq7 = dist3(0, √Sym(2)*R, R/2, 0, x11, r1, r1) |> simplify
eq8 = dist3(0, √Sym(2)*R, R/2, 0, x21, z2, r2) |> simplify;
次に，eq2, eq3 を解いて \(x_{11},\ x_{12}\) が求まる。
(ans_x11, ans_x12) = solve([eq2, eq3], (x11, x12))[1]
(2*sqrt(3)*r1/3, -sqrt(3)*r1/3)
更に，eq4, eq5, eq6 を解いて \(x_{21},\ x_{22},\ z_2\) が求まる。
(ans_x21, ans_x22, ans_z2) = solve([eq4, eq5, eq6], (x21, x22, z2))[2] # 2 of 2
(2*sqrt(3)*r2/3, -sqrt(3)*r2/3, r1*(3 + 2*sqrt(6))/3 + sqrt(3)*sqrt(3*r1^2 + 6*r1*r2 - r2^2)/3)
eq7 の \(x_{11}\) に，前に求めた解を代入し，新たな方程式を解いて \(R\) を求める。
eq7_2 = eq7(x11 => ans_x11) |> factor
\(\displaystyle \frac{2 \left(3 R^{2} - 8 \sqrt{3} R r_{1} - 3 \sqrt{2} R r_{1} + 4 r_{1}^{2} + 4 \sqrt{6} r_{1}^{2}\right)}{27}\)
ans_R = solve(eq7_2, R)[2] # 2 of 2
\(\displaystyle \frac{2 r_{1} \left(2 \sqrt{3} + 3 \sqrt{2}\right)}{3}\)
eq8
\(\displaystyle \frac{2 R^{2}}{9} - \frac{8 R x_{21}}{9} - \frac{2 \sqrt{2} R z_{2}}{9} - r_{2}^{2} + \frac{8 x_{21}^{2}}{9} + \frac{4 \sqrt{2} x_{21} z_{2}}{9} + \frac{z_{2}^{2}}{9}\)
eq8_2 = eq8(x21 => ans_x21, z2 => ans_z2, R => ans_R) |> simplify |> x -> x*27
\(12 \sqrt{6} r_{1}^{2} + 38 r_{1}^{2} - 24 \sqrt{6} r_{1} r_{2} - 26 r_{1} r_{2} - 6 \sqrt{3} r_{1} \sqrt{3 r_{1}^{2} + 6 r_{1} r_{2} - r_{2}^{2}} - 4 \sqrt{2} r_{1} \sqrt{3 r_{1}^{2} + 6 r_{1} r_{2} - r_{2}^{2}} + 4 r_{2}^{2} + 8 \sqrt{2} r_{2} \sqrt{3 r_{1}^{2} + 6 r_{1} r_{2} - r_{2}^{2}}\)
# k = r1/r2 として式を定義
Sqrt(x) = sympy.sqrt(x)
@syms k
eq8_3 = 12Sqrt(6)k^2 + 38k^2 - 24Sqrt(6)k - 26k -
(6Sqrt(3)k + 4Sqrt(2)k - 8Sqrt(2))*sqrt(3k^2 + 6k - 1) + 4
\(12 \sqrt{6} k^{2} + 38 k^{2} - 24 \sqrt{6} k - 26 k - \sqrt{3 k^{2} + 6 k - 1} \left(4 \sqrt{2} k + 6 \sqrt{3} k - 8 \sqrt{2}\right) + 4\)
# 根号を含まない部分 A と、根号の係数 B を抽出
sqrt_term = sqrt(3*k^2 + 6*k - 1)
\(\sqrt{3 k^{2} + 6 k - 1}\)
B = eq8_3.coeff(sqrt_term)
\(- 6 \sqrt{3} k - 4 \sqrt{2} k + 8 \sqrt{2}\)
A = simplify(eq8_3 - B * sqrt_term)
\(12 \sqrt{6} k^{2} + 38 k^{2} - 24 \sqrt{6} k - 26 k + 4\)
# A^2 - B^2 * D = 0 を解く
eq8_4 = A^2 - B^2 * (3*k^2 + 6*k - 1) |> factor
\(16 \left(k - 1\right) \left(48 \sqrt{6} k^{3} + 118 k^{3} - 105 \sqrt{6} k^{2} - 250 k^{2} + 18 \sqrt{6} k + 60 k - 9\right)\)
ans_k = solve(eq8_4, k)[4].evalf()
\(1.93484692283495\)
\(k=r_1/r_2\) としたので，もとに戻すと \(r_1 = 1.93484692283495 \cdot r_2\) である。
ここが，今日のハイライト！！
1.93484692283495 という数値が，どういうものか SymPy の nsimplify() で探索する。
sympy.nsimplify(ans_k, [Sqrt(6)]) |> factor
\(\displaystyle \frac{3 \left(\sqrt{6} + 4\right)}{10}\)
\(\displaystyle 1.93484692283495 = \frac{3\left(\sqrt{6} + 4\right)}{10}\) である。
術は \( (\sqrt{0.54} + 1.2)\cdot 155 = \displaystyle \frac{3\left(\sqrt{6} + 4\right)}{10}\cdot 155\) で一致している。
(sqrt(Sym(54)/100) + Sym(12)/10) |> factor
\(\displaystyle \frac{3 \left(\sqrt{6} + 4\right)}{10}\)
数値解を求める
最後に nsimplify が使えるなら，数値解を求めるのが簡単である。
include("julia-source.txt"); # julia-source.txt ソース
function driver(r2)
function H(u)
function parameters()
eq1 = (x12 - x11)^2 + r1^2 - 4r1^2
eq2 = x11^2 + (z1 - r1)^2 - 4r1^2
eq3 = x12^2 + r1^2 + (z1 - r1)^2 - 4r1^2
eq4 = (x22 - x21)^2 + r2^2 - 4r2^2
eq5 = x21^2 + (z2 - z1)^2 - (r1 + r2)^2
eq6 = x22^2 + r2^2 + (z2 - z1)^2 - (r1 + r2)^2
eq7 = dist3(0, √2R, R/2, 0, x11, r1, r1)
eq8 = dist3(0, √2R, R/2, 0, x21, z2, r2)
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8]
end;
(R, r1, x11, x12, z1, x21, x22, z2) = u
return parameters()
end;
iniv = BigFloat[10, 2, 2.2, -1.1, 5, 1.2, -0.6, 7.6]*r2
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
r2 = 1/2
(R, r1, x11, x12, z1, x21, x22, z2) = driver(r2)
@printf("乙円の直径が %g のとき，甲円の直径は %.15g である。\n", 2r2, 2r1)
乙円の直径が 1 のとき，甲円の直径は 1.93484692283495 である。
sympy.nsimplify(1.93484692283495, [sqrt(Sym(6))]) |> simplify |> factor
\(\displaystyle \frac{3 \left(\sqrt{6} + 4\right)}{10}\)
描画関数プログラムのソースを見る
function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAMincho")
(R, r1, x11, x12, z1, x21, x22, z2) = driver(r2)
p1 = plot(xlabel="x-axis", ylabel="y-axis")
plot!([0.5, 0.5, -1, 0.5].*R, [-√3, √3, 0, -√3].*R/2, color=:green, lw=0.5)
rotate(x11, 0, r1)
circle(0, 0, r1)
rotate(x21, 0, r2, :blue)
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(xlabel="x-axis", ylabel="z-axis")
plot!([0.5, 0, -1, 0.5].*R, [0, √2, 0, 0].*R, color=:green, lw=0.5)
circle(x11, r1, r1)
circle(x12, r1, r1)
circle(0, z1, r1)
circle(x21, z2, r2, :blue)
circle(x22, z2, r2, :blue)
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)
point(0, √2R, "√2R", :green, :left, :bottom, delta=delta/2)
point(R/2, 0, "R/2", :green, :left, :bottom, delta=delta/2)
xlims!(-R - delta, R/2 + 8delta)
end
plot(p1, p2)
end;
draw(155/2, true)
「算額あれこれ」の全ページの索引
以下のアイコンをクリックして応援してください
