算額（その1189）

算額（その1189）

(12) 京都市伏見区御香宮門前町 御香宮神社（ごこうのみやじんじゃ）文久3年(1863)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円6個，直線上

直線上に互いに隣同士外接する甲，乙，丙，丁，戊の 5 個の円が載っており，それらの全てに外接する己円が載っている。甲，乙，丙円の直径がそれぞれ 4 寸，1.2 寸，1 寸のとき，戊円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, y6)
己円の半径と中心座標を r6, (0, y6)
とおき，以下の連立方程式を解く。
なお，上の図は甲，乙，丙円の直径がそれぞれ 4 寸，1.29 寸，1 寸のときものである。乙円の直径が違うだけで，図は大きく変わる。与えられた条件のときの図は，後で掲示する。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive, r4::positive, x4::positive,
     r5::positive, r6::positive, x6::positive, y6::positive
#@syms r1, x1, r2, x2, r3, x3, r4, x4, r6, x6, y6, r5
#(r1, r2, r3) = (4//2, 12//20, 1//2)
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (x1 - x6)^2 + (r1 - y6)^2 - (r1 + r6)^2 |> expand
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x2 - x6)^2 + (r2 - y6)^2 - (r2 + r6)^2 |> expand
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand
eq6 = (x3 - x6)^2 + (r3 - y6)^2 - (r3 + r6)^2 |> expand
eq7 = (x4 - x6)^2 + (r4 - y6)^2 - (r4 + r6)^2 |> expand
eq8 = x4^2 + (r4 - r5)^2 - (r4 + r5)^2 |> expand
eq9 = x6^2 + (y6 - r5)^2 - (r6 + r5)^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)
   (x1, x2, x3, r4, x4, r5, r6, x6, y6) = u
   return [
       -4*r1*r2 + x1^2 - 2*x1*x2 + x2^2,  # eq1
       -2*r1*r6 - 2*r1*y6 - r6^2 + x1^2 - 2*x1*x6 + x6^2 + y6^2,  # eq2
       -4*r2*r3 + x2^2 - 2*x2*x3 + x3^2,  # eq3
       -2*r2*r6 - 2*r2*y6 - r6^2 + x2^2 - 2*x2*x6 + x6^2 + y6^2,  # eq4
       -4*r3*r4 + x3^2 - 2*x3*x4 + x4^2,  # eq5
       -2*r3*r6 - 2*r3*y6 - r6^2 + x3^2 - 2*x3*x6 + x6^2 + y6^2,  # eq6
       -2*r4*r6 - 2*r4*y6 - r6^2 + x4^2 - 2*x4*x6 + x6^2 + y6^2,  # eq7
       -4*r4*r5 + x4^2,  # eq8
       -2*r5*r6 - 2*r5*y6 - r6^2 + x6^2 + y6^2,  # eq9
   ]
end;
(r1, r2, r3) = (4, 1.2, 1)./2
#(r1, r2, r3) = (4.0, 1.29, 1.0)./2
iniv = BigFloat[11.5, 9.3, 8.2, 0.94, 6.8, 12.5, 1, 8.5, 2]
res = nls(H, ini=iniv)

   ([11.502173707608492, 9.311283477587828, 8.215838362577495, 0.9375000000000001, 6.84653196881458, 12.500000000000012, 1.0125, 8.489699641330079, 1.9874999999999998], true)

甲，乙，丙円の直径がそれぞれ 4 寸，1.2 寸，1 寸のとき，戊円の半径は 12.500000000000012 である（直径は 25 寸）。
このときの図は以下のようになる。

「術」でも，戊円の直径は 25 寸としており，間違いではない。

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

   r1 = 2;  x1 = 11.5022
   r2 = 0.6;  x2 = 9.31128
   r3 = 0.5;  x3 = 8.21584
   r4 = 0.9375;  x4 = 6.84653
   r5 = 12.5
   r6 = 1.0125;  x6 = 8.4897;  y6 = 1.9875

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, x2, x3, r4, x4, r5, r6, x6, y6) = res[1]
   @printf("r1 = %g;  x1 = %g\n", r1, x1)
   @printf("r2 = %g;  x2 = %g\n", r2, x2)
   @printf("r3 = %g;  x3 = %g\n", r3, x3)
   @printf("r4 = %g;  x4 = %g\n", r4, x4)
   @printf("r5 = %g\n", r5)
   @printf("r6 = %g;  x6 = %g;  y6 = %g\n", r6, x6, y6)
   plot()
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :magenta)
   circle(x4, r4, r4, :purple)
   circle(0, r5, r5, :green)
   circle(x6, y6, r6, :orange)
   segment(-r5, 0, x1+r1, 0, :black)
   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(x1, r1, "甲:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(x2, r2, "乙:r2,(x2,r2)", :blue, :left, delta=-delta/2)
       point(x3, r3, "丙:r3,(x3,r3)", :magenta, :center, delta=-6delta)
       point(x4, r4, "丁:r4,(x4,r4)", :purple, :right, delta=-delta/2)
       point(0, r5, "戊:r5,(0,r5)", :green, :center, delta=-delta/2)
       point(x6, y6, "己:r6,(x6,y6)", :orange, :center, delta=-delta/2)
   end
end;