算額（その2095）
二十七 群馬県太田市細谷 冠稲荷神社 文化11年(1814)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，長方形
長方形の中に，甲乙丙丁戊己の 6 個の円を容れる。甲円，丁円の直径がそれぞれ 169 寸，36 寸のとき，己円の直径はいかほどか。
長方形の長辺と短辺を a, b; b = 2r1
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, b - r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, r5)
己円の半径と中心座標を r6, (r6, r6)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
x2::positive, r3::positive, x3::positive, r4::positive, x4::positive,
r5::positive, x5::positive, r6::positive;
b = 2r1
eq1 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (b - r3 - r2)^2 - (r2 + r3)^2
eq4 = (x4 - x3)^2 + (b - r3 - r4)^2 - (r4 + r3)^2
eq5 = (x5 - x3)^2 + (b - r3 - r5)^2 - (r5 + r3)^2
eq6 = (r6 - x3)^2 + (b - r3 - r6)^2 - (r6 + r3)^2
eq7 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq8 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
eq9 = (x5 - r6)^2 + (r5 - r6)^2 - (r5 + r6)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (a, r2, x2, r3, x3, x4, r5, x5, r6))
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, r2, x2, r3, x3, x4, r5, x5, r6) = u
return [
(r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - x2)^2, # eq1
(r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - x3)^2, # eq2
-(r2 + r3)^2 + (x2 - x3)^2 + (2*r1 - r2 - r3)^2, # eq3
-(r3 + r4)^2 + (-x3 + x4)^2 + (2*r1 - r3 - r4)^2, # eq4
-(r3 + r5)^2 + (-x3 + x5)^2 + (2*r1 - r3 - r5)^2, # eq5
-(r3 + r6)^2 + (r6 - x3)^2 + (2*r1 - r3 - r6)^2, # eq6
(r2 - r4)^2 - (r2 + r4)^2 + (x2 - x4)^2, # eq7
(r4 - r5)^2 - (r4 + r5)^2 + (x4 - x5)^2, # eq8
(r5 - r6)^2 - (r5 + r6)^2 + (-r6 + x5)^2, # eq9
]
end;
(r1, r4) = (169/2, 36/2)
iniv = BigFloat[350, 27, 171, 67, 115, 127, 20, 89, 36]
res = nls(H, ini=iniv)
([350.3565608160057, 26.68421052631579, 170.88675952162194, 66.89583333333333, 115.48770876656474, 127.05454353959865, 19.6047261009667, 89.48407269786442, 36.202314049586775], true)
甲円の直径が 169 寸，丁円の直径が 36 寸のとき，己円の直径は 72.4046 寸である。
すべてのパラメータは以下のとおりである。
a = 350.357; b = 169; r1 = 84.5; r2 = 26.6842; x2 = 170.887; r3 = 66.8958; x3 = 115.488; r4 = 18; x4 = 127.055; r5 = 19.6047; x5 = 89.4841; r6 = 36.2023
function draw(r1, r4, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r2, x2, r3, x3, x4, r5, x5, r6) = res[1]
b = 2r1
@printf("甲円の直径が %g，丁円の直径が %g のとき，己円の直径は %g である。\n", 2r1, 2r4, 2r6)
@printf("a = %g; b = %g; r1 = %g; r2 = %g; x2 = %g; r3 = %g; x3 = %g; r4 = %g; x4 = %g; r5 = %g; x5 = %g; r6 = %g\n", a, b, r1, r2, x2, r3, x3, r4, x4, r5, x5, r6)
plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
circle(a - r1, r1, r1)
circle(x2, r2, r2, :magenta)
circle(x3, b - r3, r3, :blue)
circle(x4, r4, r4, :orange)
circle(x5, r5, r5, :purple)
circle(r6, r6, r6, :brown)
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(a, b, "(a,b)", :green, :right, :bottom, delta=2delta)
point(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-2delta)
point(x2, r2, "乙円:r2,(x2,r2)", :magenta, :left, :bottom, delta=2delta)
point(x3, b - r3, "丙円:r3,(x3,b-r3)", :blue, :center, delta=-2delta)
point(x4, r4, " 丁円:r4,(x4,r4)", :orange, :left, :vcenter)
point(x5, r5, "戊円:r5,(x5,r5) ", :purple, :right, :vcenter)
point(r6, r6, "己円:r6\n(r6,r6) ", :brown, :center, :bottom, delta=2delta)
end
end;
draw(169/2, 36/2, true)
#和算
