算額（その1045）

算額（その1045）

八十七 室根村矢越 矢越弥栄神社 昭和4年(1929)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

外円の中に甲円，乙円，丙円，丁円が入っている。甲円と丙円の直径がそれぞれ 2 寸，6 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (r3, y3)
丁円の半径と中心座標を r4, (0, y4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, y3::negative,
     r4::positive, y4::positive
eq1 = r1^2 + y1^2 - (R - r1)^2 |> expand
eq2 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq3 = r3^2 + y3^2 - (R - r3)^2 |> expand
eq4 = (x2 - r1)^2 + (y1 - y2)^2 - (r1 + r2)^2 |> expand
eq5 = r1^2 + (y1 - y4)^2 - (r1 + r4)^2 |> expand
eq6 = (x2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq7 = x2^2 + (y4 - y2)^2 - (r2 + r4)^2 |> expand
eq8 = r3^2 + (y4 - y3)^2 - (r3 + r4)^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)
   (R, y1, r2, x2, y2, y3, r4, y4) = u
   return [
       -R^2 + 2*R*r1 + y1^2,  # eq1
       -R^2 + 2*R*r2 - r2^2 + x2^2 + y2^2,  # eq2
       -R^2 + 2*R*r3 + y3^2,  # eq3
       -2*r1*r2 - 2*r1*x2 - r2^2 + x2^2 + y1^2 - 2*y1*y2 + y2^2,  # eq4
       -2*r1*r4 - r4^2 + y1^2 - 2*y1*y4 + y4^2,  # eq5
       -r2^2 - 2*r2*r3 - 2*r3*x2 + x2^2 + y2^2 - 2*y2*y3 + y3^2,  # eq6
       -r2^2 - 2*r2*r4 - r4^2 + x2^2 + y2^2 - 2*y2*y4 + y4^2,  # eq7
       -2*r3*r4 - r4^2 + y3^2 - 2*y3*y4 + y4^2,  # eq8
   ]
end;
(r1, r3) = (2, 6)./2
iniv = BigFloat[6, 5, 1, 3, 4, -1, 2, 3]
res = nls(H, ini=iniv)

   ([6.147114317029974, 5.049038105676658, 1.5, 3.0, 3.549038105676658, -0.950961894323342, 1.647114317029974, 2.598076211353316], true)

甲円，丙円の直径が 2, 6 のとき，乙円の直径は 3 である。

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

r1 = 1;  r3 = 3;  R = 6.14711;  y1 = 5.04904;  r2 = 1.5;  x2 = 3;  y2 = 3.54904;  y3 = -0.950962;  r4 = 1.64711;  y4 = 2.59808

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (2, 6)./2
   (R, y1, r2, x2, y2, y3, r4, y4) = res[1]
   @printf("甲円，丙円の直径が %g, %g のとき，乙円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r3 = %g;  R = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  y3 = %g;  r4 = %g;  y4 = %g\n", r1, r3, R, y1, r2, x2, y2, y3, r4, y4)
   plot()
   circle(0, 0, R)
   circle2(r1, y1, r1, :blue)
   circle2(x2, y2, r2, :green)
   circle2(r3, y3, r3, :orange)
   circle(0, y4, r4, :magenta)
   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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
   end
end;