算額（その880）
六十四 加須市不動岡 総願寺 慶応二丙寅(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
二等辺三角形内に，甲円 1 個，乙円 4 個，丙円 2 個が入っている。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。
二等辺三角形の底辺と高さを 2a，b
甲円の半径と中心座標を r1, (0, y1)
乙円の半径と中心座標を r2, (r2, r2), (3r2, r2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms d, a::positive, b::positive,
r1::positive, y1::positive,
r2::positive, r3::positive, y3::positive
eq1 = numerator(apart(dist(a, 0, 0, b, 0, y1) - r1^2, d))
eq2 = numerator(apart(dist(a, 0, 0, b, 3r2, r2) - r2^2, d))
eq3 = numerator(apart(dist(a, 0, 0, b, 2r2, y3) - r3^2, d))
eq4 = r2^2 + (y1 - r2)^2 - (r1 + r2)^2 |> expand
eq5 = 4r2^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq6 = r2^2 + (y3 - r2)^2 - (r2 + r3)^2 |> expand;
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(a, b, r1, y1, r3, y3) = u
return [
a^2*b^2 - 2*a^2*b*y1 - a^2*r1^2 + a^2*y1^2 - b^2*r1^2, # eq1
a^2*b^2 - 2*a^2*b*r2 - 6*a*b^2*r2 + 6*a*b*r2^2 + 8*b^2*r2^2, # eq2
a^2*b^2 - 2*a^2*b*y3 - a^2*r3^2 + a^2*y3^2 - 4*a*b^2*r2 + 4*a*b*r2*y3 + 4*b^2*r2^2 - b^2*r3^2, # eq3
-r1^2 - 2*r1*r2 + r2^2 - 2*r2*y1 + y1^2, # eq4
-r1^2 - 2*r1*r3 + 4*r2^2 - r3^2 + y1^2 - 2*y1*y3 + y3^2, # eq5
r2^2 - 2*r2*r3 - 2*r2*y3 - r3^2 + y3^2, # eq6
]
end;
r2 = 1/2
iniv = BigFloat[3, 3, 1, 1.5, 1, 2]
res = nls(H, ini=iniv)
([2.726702475495804, 2.665648262572091, 0.73841681234051, 1.6329943517456873, 0.3542486889354094, 1.1926332525571277], true)
乙円の直径が 1 のとき，甲円の直径は 1.47683362468102 である。
算額では「術曰置二十八個開平方加八個除以球個得甲円径合問」で，数式で表すと (sqrt(28) + 8)/9 = 1.47683362468102 となり，数値解が正しいことがわかる。
その他のパラメータは以下のとおりである。
a = 2.7267; b = 2.66565; r1 = 0.738417; y1 = 1.63299; r3 = 0.354249; y3 = 1.19263
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 1/2
(a, b, r1, y1, r3, y3) = res[1]
@printf("乙円の直径が %g のとき，甲円の直径は %.15g である\n", 2r2, 2r1)
@printf("a = %g; b = %g; r1 = %g; y1 = %g; r3 = %g; y3 = %g\n", a, b, r1, y1, r3, y3)
plot([a, 0, -a, a], [0, b, 0, 0], color=:green, lw=0.5)
circle(0, y1, r1, :magenta)
circle2(r2, r2, r2)
circle2(3r2, r2, r2)
circle2(2r2, y3, r3, :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, y1, "甲円:r1,(0, y1)", :magenta, :center, delta=-delta)
point(r2, r2, "乙円:r2\n(r2,r2)", :red, :center, delta=-delta)
point(3r2, r2, "乙円:r2\n(3r2,r2)", :red, :center, delta=-delta)
point(2r2, y3, " 丙円:r3,(2r2,y3)", :blue, :left, :vcenter)
point(a, 0, "a", :green, :left, :bottom, delta=delta)
point(0, b, " b", :green, :left, :vcenter)
end
end;
