算額（その600）

算額（その600）

長崎市 鎮西大社諏訪神社 明治20年(1887)
米光丁: 長崎県の和算の概説
http://hyonemitsu.web.fc2.com/Nagasakiwasan.pdf

問題 13. 外円の中に大円，中円，小円を入れる。大円の直径が 36 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
中円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, x3, y3)
斜線と外円の交点座標を (x, -sqrt(R^2 - x^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, R::positive, r1::positive, r2::positive, r3::positive, x3::positive, y3::positive
eq1 = r1 + r2 - R
eq2 =x3^2 + (y3 - R + r2)^2 - (r2 - r3)^2
eq3 = (y3 - R + r2)/x3*(-R -sqrt(R^2 - x^2))/x + 1
eq4 = (r2 - 2r3)/(R - r2) - x/sqrt(x^2 +(sqrt(R^2 - x^2) + R)^2)
eq5 = distance(0, R, x, -sqrt(R^2 - x^2), x3, y3) - r3^2
eq6 = distance(0, R, x, -sqrt(R^2 - x^2), 0, r1 - R) - r1^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6], (x, R, r2, r3, x3, y3))

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 v, r.f_converged
end;

function H(u)
   (x, R, r2, r3, x3, y3) = u
   return [
       -R + r1 + r2,  # eq1
       x3^2 - (r2 - r3)^2 + (-R + r2 + y3)^2,  # eq2
       1 + (-R - sqrt(R^2 - x^2))*(-R + r2 + y3)/(x*x3),  # eq3
       -x/sqrt(x^2 + (R + sqrt(R^2 - x^2))^2) + (r2 - 2*r3)/(R - r2),  # eq4
       -r3^2 + (x3 - x*(R^2 - R*y3 + R*sqrt(R^2 - x^2) + x*x3 - y3*sqrt(R^2 - x^2))/(2*R*(R + sqrt(R^2 - x^2))))^2 + (y3 - (R^2 + R*y3 - R*sqrt(R^2 - x^2) - x*x3 + y3*sqrt(R^2 - x^2))/(2*R))^2,  # eq5
       -r1^2 + (-x + r1*x/(2*R))^2 + (-R + r1 - (R^2 - (R + sqrt(R^2 - x^2))*(2*R - r1)/2)/R)^2,  # eq6
   ]
end;

r1 = 36//2
iniv = BigFloat[22, 35, 17, 5, 10, 20]
res = nls(H, ini=iniv)

   (BigFloat[22.62741699796952078082701958735516925711475000603116912031953235176946800128799, 36.00000000000000000000000000000000000000000000000000023997217352842064306535952, 18.0000000000000000000000000000000000000000000000000002399721735284206430653598, 6.000000000000000000000000000000000000000000000000000298124281531485670596315888, 11.3137084989847603904135097936775846285573750030155858282654603937949442060392, 22.00000000000000000000000000000000000000000000000000175488066987747432400063655], true)

小円の直径は 12 寸である。

その他のパラメータは以下の通り。
x = 22.6274;  R = 36;  r2 = 18;  r3 = 6;  x3 = 11.3137;  y3 = 22

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 36//2
   (x, R, r2, r3, x3, y3) = res[1]
   @printf("小円の直径 = %g;  x = %g;  R = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", 2r3, x, R, r2, r3, x3, y3)
   plot()
   circle(0, 0, R, :blue)
   circle(0, r1 - R, r1, :green)
   circle(0, R - r2, r2, :orange)
   circle(x3, y3, r3)
   segment(0, R, x, -sqrt(R^2 - x^2))
   segment(0, R, -x, -sqrt(R^2 - x^2))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, R - r2, " R-r2")
       point(x3, y3, "(x3,y3)")
   end
end;