算額（その2120）

算額（その2120）

九 武州崎玉郡騎西町 久伊豆神社 文化4年(1807)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，接線4本，直線上
#Julia, #SymPy, #算額, #和算

直線上に接線 4 本で仕切られている区画に 8 個の円を容れる。乙円と丙円の直径がそれぞれ 44 寸，20 寸のとき，丁円の直径はいかほどか。

甲円の半径と中心座標を r1,(0, y1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4), (x42, r2)
接線が通る 2 点の座標対を [(-a, 0), (x01, y01)], [(b, 0), (x02, y02)]
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, y1::positive, r2::positive, x2::positive,
     r3::positive, r4::positive, x4::positive, x42::positive,
     a::positive, b::positive, x01::positive, y01::positive,
     x02::positive, y02::positive;
eq1 = dist2(-a, 0, x01, y01, 0, y1, r1)
eq2 = dist2(-a, 0, x01, y01, x2, r2, r2)
# eq3 = dist2(-a, 0, x01, y01, 0, r3, r3)
# eq4 = dist2(-a, 0, x01, y01, x42, r2, r4)
eq3 = r3/a - r2/(a + x2)
eq4 = r4/(a - x4) - r3/a
eq5 = dist2(b, 0, x02, y02, 0, y1, r1)
eq6 = dist2(b, 0, x02, y02, x2, r2, r2)
eq7 = dist2(b, 0, x02, y02, 0, r3, r3)
eq8 = dist2(b, 0, x02, y02, x4, r4, r4)
eq9 = dist2(b, 0, x02, y02, x42, r2, r4)
eq10 = dist2(a, 0, -x01, y01, x2, r2, r2)
eq11 = dist2(a, 0, -x01, y01, x4, r4, r4)
eq12 = dist2(a, 0, -x01, y01, x42, r2, r4);

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)
   (r1, y1, x2, r4, x4, x42, a, b, x01, y01, x02, y02) = u
   return [
-a^2*r1^2 + a^2*y01^2 - 2*a^2*y01*y1 + a^2*y1^2 - 2*a*r1^2*x01 - 2*a*x01*y01*y1 + 2*a*x01*y1^2 - r1^2*x01^2 - r1^2*y01^2 + x01^2*y1^2,  # eq1
y01*(-2*a^2*r2 + a^2*y01 - 2*a*r2*x01 - 2*a*r2*x2 + 2*a*x2*y01 - r2^2*y01 - 2*r2*x01*x2 + x2^2*y01),  # eq2
-r2/(a + x2) + r3/a,  # eq3
r4/(a - x4) - r3/a,  # eq4
-b^2*r1^2 + b^2*y02^2 - 2*b^2*y02*y1 + b^2*y1^2 + 2*b*r1^2*x02 + 2*b*x02*y02*y1 - 2*b*x02*y1^2 - r1^2*x02^2 - r1^2*y02^2 + x02^2*y1^2,  # eq5
y02*(-2*b^2*r2 + b^2*y02 + 2*b*r2*x02 + 2*b*r2*x2 - 2*b*x2*y02 - r2^2*y02 - 2*r2*x02*x2 + x2^2*y02),  # eq6
y02*(-2*b^2*r3 + b^2*y02 + 2*b*r3*x02 - r3^2*y02),  # eq7
y02*(-2*b^2*r4 + b^2*y02 + 2*b*r4*x02 + 2*b*r4*x4 - 2*b*x4*y02 - r4^2*y02 - 2*r4*x02*x4 + x4^2*y02),  # eq8
b^2*r2^2 - 2*b^2*r2*y02 - b^2*r4^2 + b^2*y02^2 - 2*b*r2^2*x02 + 2*b*r2*x02*y02 + 2*b*r2*x42*y02 + 2*b*r4^2*x02 - 2*b*x42*y02^2 + r2^2*x02^2 - 2*r2*x02*x42*y02 - r4^2*x02^2 - r4^2*y02^2 + x42^2*y02^2,  # eq9
y01*(-2*a^2*r2 + a^2*y01 - 2*a*r2*x01 + 2*a*r2*x2 - 2*a*x2*y01 - r2^2*y01 + 2*r2*x01*x2 + x2^2*y01),  # eq10
y01*(-2*a^2*r4 + a^2*y01 - 2*a*r4*x01 + 2*a*r4*x4 - 2*a*x4*y01 - r4^2*y01 + 2*r4*x01*x4 + x4^2*y01),  # eq11
a^2*r2^2 - 2*a^2*r2*y01 - a^2*r4^2 + a^2*y01^2 + 2*a*r2^2*x01 - 2*a*r2*x01*y01 + 2*a*r2*x42*y01 - 2*a*r4^2*x01 - 2*a*x42*y01^2 + r2^2*x01^2 + 2*r2*x01*x42*y01 - r4^2*x01^2 - r4^2*y01^2 + x42^2*y01^2,  # eq12
   ]
end;

(r2, r3) = (44/2, 20/2)
iniv = BigFloat[28, 55, 40, 5.5, 15, 10, 33, 6.6, 56, 59, 39, 77]
res = nls(H, ini=iniv)

   ([27.5, 55.0, 39.7994974842648, 5.5, 14.9248115565993, 9.9498743710662, 33.166247903554, 6.6332495807108, 54.652326389566944, 58.25225211084647, 39.43354494756409, 77.70448053191785], true)

乙円と丙円の直径がそれぞれ 44 寸，20 寸のとき，丁円の直径は 2*5.5 = 11 寸である。

すべてのパラメータは以下のとおりである。

    r2 = 22;  r3 = 10;  r1 = 27.5;  y1 = 55;  x2 = 39.7995;  r4 = 5.5;  x4 = 14.9248;  x42 = 9.94987;  a = 33.1662;  b = 6.63325;  x01 = 54.6523;  y01 = 58.2523;  x02 = 39.4335;  y02 = 77.7045

function draw(r2, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y1, x2, r4, x4, x42, a, b, x01, y01, x02, y02) = [27.5, 55.0, 39.7994974842648, 5.5, 14.9248115565993, 9.9498743710662, 33.166247903554, 6.6332495807108, 54.652326389566944, 58.25225211084647, 39.43354494756409, 77.70448053191785]
   plot()
   circle(0, y1, r1)
   circle2(x2, r2, r2, :blue)
   circle(0, r3, r3, :green)
   circle2(x4, r4, r4, :magenta)
   circle2(x42, r2, r4, :magenta)
   segment(-a, 0, x01, y01)
   segment(b, 0, x02, y02)
   segment(a, 0, -x01, y01)
   segment(-b, 0, -x02, y02)
   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(x01, y01, " (x01,y01)", :black, :left, :vcenter)
       point(x02, y02, " (x02,y02)", :black, :left, :vcenter)
       point(-a, 0, "-a", :black, :center, delta=-1.5delta)
       point(b, 0, "b", :black, :center, delta=-1.5delta)
       point(0, y1, "甲円:r1,(0,y1)", :red, :center, delta=-delta/2)
       point(x2, r2, "乙円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(0, r3, "丙円:r3\n(0,r3)", :green, :center, delta=-delta/2)
       point(x4, r4, "丁円:r4,(x4,r4)", :black, :left, delta=-delta/2)
       point(x42, r2, "丁円:r4,(x42,r2)", :black, :left, :bottom, delta=delta)
       plot!(xlims=(-(x2 + r2 + 3delta), x2 + r2 + 13delta), ylims=(-5delta, y1 + r1 + 2delta))
   end
end;

draw(44/2, 20/2, true)