算額（その481）

算額（その481）

宮城県丸森町小斎日向 鹿島神社 明治13年

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

長方形内に 2 本の斜線で隔てられた領域に甲，乙，丙，丁の 4 円が入っている。甲円の直径が 12 寸のとき，乙円の直径を最大にしたい。乙円，丙円，丁円の直径を求めよ。

長方形の長辺の長さを a
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (x3, 2r1 - r3)
丁円の半径と中心座標を r4, (r4, y4)
とする。

乙円が最大になるのは 2 本の斜線が直交するときである。

このとき，甲円が内接する直角三角形（の半分）と，乙円，丙円，丁円の 3 個の直角三角形は相似なので，相似比は 12:6:4:3 になるのは明らかである。

長方形の長辺の長さ a，丙円と丁円の中心座標を確定するためには，以下の連立方程式を nlsolve() で解き，数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, d::positive,
     r1::positive, r2::positive, r3::positive, x3::positive,
     r4::positive, y4::positive;

eq1 = distance(0, 2r1, b, 0, a - r1, r1) - r1^2 |> simplify
eq2 = distance(0, 2r1, b, 0, r2, r2) - r2^2 |> simplify
eq3 = distance(0, 2r1, b, 0, x3, 2r1 - r3) - r3^2 |> simplify
eq4 = distance(0, 2r1, b, 0, r4, y4) - r4^2 |> simplify
eq5 = distance(0, c, d, 2r1, a - r1, r1) - r1^2 |> simplify
eq6 = distance(0, c, d, 2r1, r2, r2) - r2^2 |> simplify
eq7 = distance(0, c, d, 2r1, x3, 2r1 - r3) - r3^2 |> simplify
eq8 = distance(0, c, d, 2r1, r4, y4) - r4^2 |> simplify
eq9 = b/2r1 - (2r1 - c)/d;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (a, b, c, d, r2, r3, x3, r4, y4))

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, b, c, d, r2, r3, x3, r4, y4) = u
   return [
       4*r1^2*(a^2 - a*b - 2*a*r1 + b*r1)/(b^2 + 4*r1^2),  # eq1
       4*b*r1*(b*r1 - b*r2 - 2*r1*r2 + r2^2)/(b^2 + 4*r1^2),  # eq2
       4*r1*(-b*r3*x3 - r1*r3^2 + r1*x3^2)/(b^2 + 4*r1^2),  # eq3
       b*(4*b*r1^2 - 4*b*r1*y4 - b*r4^2 + b*y4^2 - 8*r1^2*r4 + 4*r1*r4*y4)/(b^2 + 4*r1^2),  # eq4
       (a^2*c^2 - 4*a^2*c*r1 + 4*a^2*r1^2 - 2*a*c^2*d - 2*a*c^2*r1 + 6*a*c*d*r1 + 8*a*c*r1^2 - 4*a*d*r1^2 - 8*a*r1^3 + c^2*d^2 + 2*c^2*d*r1 - 2*c*d^2*r1 - 6*c*d*r1^2 + 4*d*r1^3)/(c^2 - 4*c*r1 + d^2 + 4*r1^2),  # eq5
       d*(c^2*d - 2*c^2*r2 - 2*c*d*r2 + 4*c*r1*r2 + 2*c*r2^2 - 4*r1*r2^2)/(c^2 - 4*c*r1 + d^2 + 4*r1^2),  # eq6
       (c^2*d^2 - 2*c^2*d*x3 - c^2*r3^2 + c^2*x3^2 - 4*c*d^2*r1 + 2*c*d^2*r3 + 8*c*d*r1*x3 - 2*c*d*r3*x3 + 4*c*r1*r3^2 - 4*c*r1*x3^2 + 4*d^2*r1^2 - 4*d^2*r1*r3 - 8*d*r1^2*x3 + 4*d*r1*r3*x3 - 4*r1^2*r3^2 + 4*r1^2*x3^2)/(c^2 - 4*c*r1 + d^2 + 4*r1^2),  # eq7
       d*(c^2*d - 2*c^2*r4 - 2*c*d*y4 + 4*c*r1*r4 + 2*c*r4*y4 - d*r4^2 + d*y4^2 - 4*r1*r4*y4)/(c^2 - 4*c*r1 + d^2 + 4*r1^2),  # eq8
       b/(2*r1) - (-c + 2*r1)/d,  # eq9
   ]
end;

r1 = 12/2
iniv = [big"18.2", 9.1, 4.4, 10.3, 3.1, 2.1, 4.1, 1.6, 7.4]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[17.99999999999999999999999999999958787838359120591652042567639170665603992771524, 8.99999999999999999999999999999984155332251522567635341448728327805462644595505, 4.500000000000000000000000000000024260203115491587793400726971558949746337902191, 9.99999999999999999999999999999961117790943389634414624752326251631979937803581, 2.999999999999999999999999999999520308595440307614021124672062062280275651040769, 1.999999999999999999999999999998867711776892728650702737150651030496089805787837, 3.999999999999999999999999999999609211568160356725082605146703705988921640964954, 1.499999999999999999999999999999653325970663365188170394965090885917020592201141, 7.499999999999999999999999999999917260822554019948727938965336654899285660744732], true)

   a = 18;  b =　9;  c =　4.5;  d =　10;  r1 = 6;  r2 =　3;  r3 =　2;  x3 =　4;  r4 = 1.5;  y4 = 7.5
   乙円の直径 = 6;  丙円の直径 = 4;  丁円の直径 = 3

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 12/2
   (a, b, c, d, r2, r3, x3, r4, y4)  = res[1]
   @printf("a = %g;  b =　%g;  c =　%g;  d =　%g;  r1 = %g;  r2 =　%g;  r3 =　%g;  x3 =　%g;  r4 = %g;  y4 = %g\n",  a, b, c, d, r1, r2, r3, x3, r4, y4)
   @printf("乙円の直径 = %g;  丙円の直径 = %g;  丁円の直径 = %g\n", 2r2, 2r3, 2r4)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:blue, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, r2, r2, :green)
   circle(x3, 2r1 - r3, r3, :magenta)
   circle(r4, y4, r4, :orange)
   segment(0, c, d, 2r1)
   segment(0, 2r1, b, 0)
   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(a, 0, "a", :blue, :center, :top, delta=-1.5delta)
       point(b, 0, "b", :blue, :center, :top, delta=-1.5delta)
       point(0, c, "c ", :blue, :right, :vcenter)
       point(d, 2r1, "(d,2r1)", :blue, :center, :bottom, delta=delta)
       point(0, 2r1, "2r1 ", :blue, :right, :vcenter)
       point(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, :top, delta=-delta)
       point(r2, r2, "乙円:r2,(r2,r2)", :green, :center, :top, delta=-delta)
       point(x3, 2r1 - r3, "丙円:x3\n(x3,2r1-r3)", :magenta, :center, :bottom, delta=delta)
       point(r4, y4, "丁円:r4\n(r4,y4)", :orange, :center, :vcenter)
       plot!(xlims=(-1.5, 19), ylims=(-1, 13))
   else
       plot!(showaxis=false)
   end
end;