算額（その1405）
八七 群馬県碓氷郡松井田町峠 熊野神社 安政4年(1857)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：3次元，球2個，四角台
#Julia, #SymPy, #算額, #和算
上面，下面の一辺が 125 寸，2000 寸の正方形の四角台の中に小球と大球を容れる。小球と大球は互いに接するほかに，小球は側面と上面に 5 点で，大球は側面と下面に 5 点で接している。小球の直径はいかほどか。
上面の正方形，下面の正方形の一辺の長さをそれぞれ 2j, 2k
四角台の高さを i
四角台を含む四角錐の高さを h
大球の半径と中心座標を r1, (0, 0, r1)
小球の半径と中心座標を r2, (0, 0, 2r1 + r2)
とおく。
3次元の図形をx軸の正の方向から原点方向を見た y-z 平面へ投影された図を考える。
上面と下面の正方形の一辺の中点 i, k を結ぶ直線は小球と大球の共通接線になる。o, i, j, k は y-z 平面に平行な平面上にある。これを y-z 平面に投影した図は以下のようになる。
点と線の位置関係は以下の 3 本の方程式で表される。図形を描くために必要な変数は r1, r2, h の 3 変数なので，連立方程式を解けばよい。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms r1::positive, r2::positive, h::positive, j::positive, k::positive
eq1 = r2/(h - 2r1 - r2) - r1/(h - r1)
eq2 = r1/(h - r1) - k/sqrt(h^2 + k^2)
eq3 = j/(h - 2r1 - 2r2) - k/h
res = solve([eq1, eq2, eq3], (r1, r2, h))[2]; # 2 of 4
r1, r2, h ともに，ルートの中は j と k の高次式になっている。
小球の半径は，以下のようである。
res[2] |> display
res[2] |> println
sqrt((2*j^(21/2)*k^4 - 12*j^(19/2)*k^5 + 30*j^(17/2)*k^6 - 40*j^(15/2)*k^7 + 30*j^(13/2)*k^8 - 12*j^(11/2)*k^9 + 2*j^(9/2)*k^10 + j^11*k^(7/2) - 5*j^10*k^(9/2) + 9*j^9*k^(11/2) - 5*j^8*k^(13/2) - 5*j^7*k^(15/2) + 9*j^6*k^(17/2) - 5*j^5*k^(19/2) + j^4*k^(21/2))/(j^(19/2)*k^3 - 5*j^(17/2)*k^4 + 9*j^(15/2)*k^5 - 5*j^(13/2)*k^6 - 5*j^(11/2)*k^7 + 9*j^(9/2)*k^8 - 5*j^(7/2)*k^9 + j^(5/2)*k^10 + 2*j^9*k^(7/2) - 12*j^8*k^(9/2) + 30*j^7*k^(11/2) - 40*j^6*k^(13/2) + 30*j^5*k^(15/2) - 12*j^4*k^(17/2) + 2*j^3*k^(19/2)))
SymPy では簡約化できないので，手動で簡約化する。まず根号の中身の分子と分母をそれぞれ変数変換して因子分解し，再度分数式に戻して平方根を取る。
# 分子を変数変換
@syms t, u
num = res[2].args[1](j^(1//2) => t, k^(1//2) => u) |> numerator |> factor
num |> println
t^8*u^7*(t - u)^6*(t + u)^8
# 分母を変数変換
den = res[2].args[1](j^(1//2) => t, k^(1//2) => u) |> denominator |> factor
den |> println
t^5*u^6*(t - u)^6*(t + u)^8
# 分数式に戻す
r2 = num/den
r2 |> println
t^3*u
# 変数の逆変換
r2 = sqrt(r2(t => j^(1//2), u => k^(1//2))) |> factor
r2 |> println
j^(3/4)*k^(1/4)
小球の直径は上面の正方形の一辺の長さの 3 乗と下面の正方形の一辺の長さを掛けて，4 乗根を求めることで得られる。
# 実際の数値を代入して答えを求める
r2 = fourthroot(125^3 * 2000)
250.0
大球の直径は小球の直径の 4 倍，四角台の高さは大球と小球の直径の和（小球の直径の 5 倍）である。
# 大球の直径は小球の直径の何倍か？
res[1]/res[2] |> println
(res[1]/res[2])(j => 125//2, k => 2000//2) |> println
sqrt(k)/sqrt(j)
4
# 四角錐の高さは小球の直径の何倍か？
res[3]/res[2] |> println
(res[3]/res[2])(j => 125//2, k => 2000//2) |> println
2*(sqrt(j)*k + k^(3/2))/(sqrt(j)*(-j + k))
32/3
1. 連立方程式の解き方（別法）
割り算を避けるなど，なるべく式を簡潔にして，逐次的に解を求めていく。
@syms r1::positive, r2::positive, h::positive, j::positive, k::positive
eq1 = r2*(h - r1) - r1*(h - 2r1 - r2) |> expand
eq2 = r1*sqrt(h^2 + k^2) - k*(h - r1) |> expand
eq3 = j*h - k*(h - 2r1 - 2r2) |> expand;
eq2 は r2 を含まないので，そのぶん簡潔なので，まず eq1, eq2 から r1, h を求める（式に r2 が含まれる）。
res = solve([eq1, eq2], (r1, h))[1]
(k^(2/3)*r2^(1/3), 2*k^(4/3)*r2^(2/3)/(k^(2/3)*r2^(1/3) - r2))
eq3 の r1, h に代入して r2 を求める。
eq13 = eq3(r1 => res[1], h => res[2]) |> simplify |> numerator;
ans_r2 = solve(eq13, r2)[1]
ans_r2 |> println
j^(3/4)*k^(1/4)
求められた r2 を，先に求めた r1, h の r2 に代入する。
@syms j, k
j = 125/2
k = 2000/2
r2 = j^(3/4)*k^(1/4)
r1 = k^(2/3)*r2^(1/3)
h = 2*k^(4/3)*r2^(2/3)/(k^(2/3)*r2^(1/3) - r2)
println("r2 = $r2\nr1 = $r1\nh = $h")
r2 = 125.0
r1 = 499.9999999999999
h = 1333.3333333333328
function draw(j, k, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = fourthroot(j^3 * k)
r1 = 4r2
i = 2r1 + 2r2
h = r2*32/3
@printf("上面と下面の正方形の一辺の長さが %g, %g のとき，小球，大球の直径は %g, %g，四角台の高さは %g である。\n",2j, 2k, 2r2, 2r1, 2r1 + 2r2)
plot([k, 0, -k, k], [0, h, 0, 0], color=:blue, lw=0.5)
circle(0, r1, r1)
circle(0, 2r1 + r2, r2, :green)
segment(j, i, -j, i, :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, r1, "大球:r1,(0,0,r1)", :red, :center, delta=-delta)
point(0, 2r1 + r2, "小球:r2,(0,0,2r1+r2)", :green, :left, delta=-delta, deltax=-5delta)
point(j, i, "(j,0,i)", :blue, :left, :bottom, delta=delta)
point(k, 0, "(k,0,k)", :blue, :left, :bottom, delta=delta)
point(0, h, "(0,0,h)", :blue, :left, :bottom, delta=delta)
xlims!(-k - 3delta, k + 15delta)
end
end;
draw(125/2, 2000/2, true)
