算額（その686）

算額（その686）

三七 児玉郡上里村勅使河原 丹生神社 弘化3年(1846)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

和算問題あれこれ 2 令和5年11月の問題-No.2（『神壁算法』第34問）
https://gunmawasan.web.fc2.com/k-n-mondai.html

長方形内に斜線を2本入れ，区分された領域に内接円を 4 個入れる。亨円，貞円の直径がそれぞれ 44 寸，33 寸のとき，利円の直径はいかほどか。

長方形の左下頂点を原点とし，長辺の長さを a とする（短辺の長さは元円の直径 2r1）
斜線と長辺の交点座標を (c, 0 ), (b, 0), (d, 2r1)
元円の半径と中心座標を r1, (a - r1, r1)
亨円の半径と中心座標を r2, (x2, 2r1 - r2)
利円の半径と中心座標を r3, (r3, r3)
貞円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
eq1〜eq7 は円の中心から斜線への距離に関する方程式，eq8 は元円と亨円の半径の相似比についての方程式である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a, b, c, d, r1, r2, x2, r3, r4, x4

#(r2, r4) = (44, 33).//2

eq1 = dist(c, 0, d, 2r1, a - r1, r1) - r1^2
eq2 = dist(c, 0, d, 2r1, x2, 2r1 - r2) - r2^2
eq3 = dist(c, 0, d, 2r1, r3, r3) - r3^2
eq4 = dist(c, 0, d, 2r1, x4, r4) - r4^2
eq5 = dist(b, 0, 0, 2r1, x2, 2r1 - r2) - r2^2
eq6 = dist(b, 0, 0, 2r1, r3, r3) - r3^2
eq7 = dist(b, 0, 0, 2r1, x4, r4) - r4^2
eq8 = r2/x2 - r1/(a - r1);

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)
   (a, b, c, d, r1, x2, r3, x4) = u
   return [
       -r1^2 + (-2*r1*(2*r1^2 + (-c + d)*(a - c - r1))/(4*r1^2 + (-c + d)^2) + r1)^2 + (a - c - r1 - (-c + d)*(2*r1^2 + (-c + d)*(a - c - r1))/(4*r1^2 + (-c + d)^2))^2,  # eq1
       -r2^2 + (-c + x2 - (-c + d)*(2*r1*(2*r1 - r2) + (-c + d)*(-c + x2))/(4*r1^2 + (-c + d)^2))^2 + (2*r1 - 2*r1*(2*r1*(2*r1 - r2) + (-c + d)*(-c + x2))/(4*r1^2 + (-c + d)^2) - r2)^2,  # eq2
       -r3^2 + (-2*r1*(2*r1*r3 + (-c + d)*(-c + r3))/(4*r1^2 + (-c + d)^2) + r3)^2 + (-c + r3 - (-c + d)*(2*r1*r3 + (-c + d)*(-c + r3))/(4*r1^2 + (-c + d)^2))^2,  # eq3
       -r4^2 + (-2*r1*(2*r1*r4 + (-c + d)*(-c + x4))/(4*r1^2 + (-c + d)^2) + r4)^2 + (-c + x4 - (-c + d)*(2*r1*r4 + (-c + d)*(-c + x4))/(4*r1^2 + (-c + d)^2))^2,  # eq4
       -r2^2 + (-b + b*(-b*(-b + x2) + 2*r1*(2*r1 - r2))/(b^2 + 4*r1^2) + x2)^2 + (2*r1 - 2*r1*(-b*(-b + x2) + 2*r1*(2*r1 - r2))/(b^2 + 4*r1^2) - r2)^2,  # eq5
       -r3^2 + (-2*r1*(-b*(-b + r3) + 2*r1*r3)/(b^2 + 4*r1^2) + r3)^2 + (-b + b*(-b*(-b + r3) + 2*r1*r3)/(b^2 + 4*r1^2) + r3)^2,  # eq6
       -r4^2 + (-2*r1*(-b*(-b + x4) + 2*r1*r4)/(b^2 + 4*r1^2) + r4)^2 + (-b + b*(-b*(-b + x4) + 2*r1*r4)/(b^2 + 4*r1^2) + x4)^2,  # eq7
       -r1/(a - r1) + r2/x2,  # eq8
   ]
end;

(r2, r4) = (44, 33)./2
iniv = BigFloat[205, 155, 40, 150, 42, 90, 30, 90]
res = nls(H, ini=iniv)

   (BigFloat[209.9999999999999999999999999999999999999999999999999999999999999998827654490999, 157.4999999999999999999999999999999999999999999999999999999999999999096530095862, 41.99999999999999999999999999999999999999999999999999999999999999999764231208424, 153.9999999999999999999999999999999999999999999999999999999999999998816882373083, 42.0000000000000000000000000000000000000000000000000000000000000000018357614881, 87.99999999999999999999999999999999999999999999999999999999999999994964199517249, 31.49999999999999999999999999999999999999999999999999999999999999999749742450348, 91.49999999999999999999999999999999999999999999999999999999999999994744891428587], true)

元円，利円の半径はそれぞれ 42, 31.5 である。
したがって，亨円，貞円の直径がそれぞれ 44 寸，33 寸のとき，利円の直径は 63 寸である。

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

   a = 210;  b = 157.5;  c = 42;  d = 154;  r1 = 42;  r2 = 22;  x2 = 88;  r3 = 31.5;  r4 = 16.5;  x4 = 91.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d, r1, x2, r3, x4) = [189, 113, 52, 100, 54, 62, 32, 76]
   (a, b, c, d, r1, x2, r3, x4) = res[1]
   (r2, r4) = (44, 33)./2
   println("利円の直径 = $(round(Int, Float64(2r3)))")
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  r4 = %g;  x4 = %g\n", a, b, c, d, r1, r2, x2, r3, r4, x4)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(x2, 2r1 - r2, r2, :blue)
   circle(r3, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   segment(b, 0, 0, 2r1)
   segment(c, 0, d, 2r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(a - r1, r1, "元円:r1,(a-r1,r1)", :red, :center, delta=-2delta)
       point(x2, 2r1 - r2, "亨円:r2\n(x2,2r1-r2", :blue, :center, :bottom, delta=2delta)
       point(r3, r3, "利円:r3,(r3,r3)", :green, :center, delta=-2delta)
       point(x4, r4, "貞円:r4,(x4,r4)", :black, :center, delta=-2delta)
       point(a, 0, "a", :black, :center, :top, delta=-delta)
       point(b, 0, "b", :black, :center, :top, delta=-delta)
       point(c, 0, "c", :black, :center, :top, delta=-delta)
       point(d, 2r1, "(d,2r1)", :black, :center, :bottom, delta=delta)
       point(0, 2r1, " 2r1", :black, :left, :bottom, delta=delta)
       plot!(ylims=(-10, 90))
   end
end;