算額（その1295）

算額（その1295）

二十七 群馬県太田市細谷 冠稲荷神社 文化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)

#和算