算額（その1864） - 算額あれこれ

92 岩手県大船渡市盛町館下根城八幡宮年紀不明

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：3次元，球，半球
#Julia #SymPy #算額 #和算 #数学

半球の中に天球を設け，地球 2 個，人球 1 個を容れる。人球の直径が与えられたとき，半球の直径を得る術を述べよ。

半球の半径と中心座標を \(R, (0, 0)\)
天球の半径と中心座標を \(r_1, (0, 0, r_1)\)
地球の半径と中心座標を \(r_2, (x_2, r_2, r_2)\)
人球の半径と中心座標を \(r_3, (x_3, 0, z_3)\)
とおき以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive, r3::positive, x3::positive, z3::positive
r1 = R/2
# 天球と地球が外接する
eq1 = x2^2 + r2^2 + (r1 - r2)^2 - (r1 + r2)^2
# 天球と人球が外接する
eq2 = x3^2 + (z3 - r1)^2 - (r1 + r3)^2
# 地球と人球が外接する
eq3 = (x3 - x2)^2 + r2^2 + (z3 - r2)^2 - (r3 + r2)^2
# 人球が半球に内接する
eq4 = x3^2 + z3^2 - (R - r3)^2
# 地球が半球に内接する
eq5 = x2^2 + r2^2 + r2^2 - (R - r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (R, r2, x2, x3, z3))[1]

(r3*(1 + sqrt(14*sqrt(21) + 70)/2), r3*(2 + sqrt(14*sqrt(21) + 70))/8, r3*(2*sqrt(7) + 7*sqrt(2*sqrt(21) + 10))/8, r3*sqrt(2*sqrt(21) + 10), r3*(3 + sqrt(21))/2)

半球の半径 \(R\) は，人球の半径 \(r_3\) の \( (\sqrt{21} + 9)/2 = 6.7912878474779195\) 倍である。

ans_R = res[1] |> sympy.sqrtdenest |> factor
@show(ans_R)

ans_R = r3*(sqrt(21) + 9)/2

　\(\displaystyle \frac{r_{3} \left(\sqrt{21} + 9\right)}{2}\)

ans_R(r3 => 1/2).evalf() * 2

　\(\displaystyle 6.79128784747792\)

術は，「(sqrt(5.25) + 4.5)*人球径」である。5.25 = 21/4, 4.5=9/2 なので，上で得た解と同じである。

sqrt(Sym(21)/4) + Sym(9)/2

　\(\displaystyle \frac{\sqrt{21}}{2} + \frac{9}{2}\)

描画関数プログラムのソースを見る

function draw(r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r2, x2, x3, z3) = (r3*(1 + sqrt(14*sqrt(21) + 70)/2), r3*(2 + sqrt(14*sqrt(21) + 70))/8, r3*(2*sqrt(7) + 7*sqrt(2*sqrt(21) + 10))/8, r3*sqrt(2*sqrt(21) + 10), r3*(3 + sqrt(21))/2)
println( (r3, R, r2, x2, x3, z3))
r1 = R/2
p1 = plot(xlabl="x-axis", ylabel="z-axis")
circle(0, 0, R, beginangle=0, endangle=180)
circle(0, r1, r1, :green)
circle(x2, r2, r2, :blue)
circle(x3, z3, r3, :magenta)
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\n(0,0,r1)", :red, :center, delta=-delta)
point(x2, r2, "地円:r2,(x2,r2,r2)", :blue, :right, delta=-6delta, deltax=15delta)
point(x3, z3, "人円:r2,(x3,0,z3)", :magenta, :right, :bottom, delta=6delta, deltax=5delta)
end
p2 = plot(xlabl="x-axis", ylabel="y-axis")
circle(0, 0, R)
circle(0, 0, r1, :green)
circle22(x2, r2, r2, :blue)
circle(x3, 0, r3, :magenta)
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
plot(p1, p2)
end;