算額（その490）

算額（その490）

新潟県長岡市上岩井 根立寺 嘉永2年(1849)
http://www.wasan.jp/niigata/konryuji.html
涌田和芳，外川一仁: 三島根立寺の算額，長岡工業高等専門学校研究紀要，第53巻(2017)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_53/53_17wakuta.pdf

正方形内に 2 本の斜線で分割された領域に大円 1 個，小円 3 個が入っている。
大円の直径が 10 寸のとき，小円の直径はいかほどか。

正方形の左下隅を原点とし，正方形の一辺の長さを a とする
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (r2, a - r2), (r2, a - 3r2), (a - r2, a - r2)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive,
     x1::positive, r2::positive;

eq1 = distance(0, 0, b, a, r2, a - 3r2) - r2^2
eq2 = distance(0, 0, b, a, x1, r1) - r1^2
eq3 = distance(a, 0, b, a, a - r2, a - r2) - r2^2
eq4 = distance(a, 0, b, a, x1, r1) - r1^2;
# res = solve([eq1, eq2, eq3, eq4], (b, r1, x1, r2))

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)
   (b, r1, x1, r2) = u
   return [
       -r2^2 + (-b*(a^2 - 3*a*r2 + b*r2)/(a^2 + b^2) + r2)^2 + (a - a*(a^2 - 3*a*r2 + b*r2)/(a^2 + b^2) - 3*r2)^2,  # eq1
       -r1^2 + (-a*(a*r1 + b*x1)/(a^2 + b^2) + r1)^2 + (-b*(a*r1 + b*x1)/(a^2 + b^2) + x1)^2,  # eq2
       -r2^2 + (a - r2 - (a^3 - a^2*b + a*b^2 + a*b*r2 - b^2*r2)/(2*a^2 - 2*a*b + b^2))^2 + (-a*(a^2 - b*r2)/(2*a^2 - 2*a*b + b^2) + a - r2)^2,  # eq3
       -r1^2 + (x1 - (a*(a^2 - a*r1 - a*x1 + b*r1) + x1*(2*a^2 - 2*a*b + b^2))/(2*a^2 - 2*a*b + b^2))^2 + (-a*(a^2 - a*b + a*r1 - a*x1 + b*x1)/(2*a^2 - 2*a*b + b^2) + r1)^2,  # eq4
   ]
end;
a = 1
iniv = BigFloat[0.61, 0.3, 0.57, 0.14]
res = nls(H, ini=iniv)

   (BigFloat[0.6234898018587335305250048840042398106322575229856815465655969032398117579169351, 0.3079785283699041303721851029979308598027401532891747156054051562155987534336586, 0.5549581320873711914221948710064104810672818660788279185605428752567220611398606, 0.1539892641849520651860925514989654299013683741434972076685407572160588959569412], true)

   a = 1;  b = 0.62349;  r1 = 0.307979;  x1 = 0.554958;  r2 = 0.153989

小円の直径は大円の直径の 1/2 である。
大円の直径が 10 寸のとき，小円の直径は 5 寸である。

res[1][4]/res[1][2]  # r2/r1

   0.4999999999999999999999999999999999999999944720136849757617152696086889978777665

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (b, r1, x1, r2) = res[1]
   @printf("a = %g;  b = %g;  r1 = %g;  x1 = %g;  r2 = %g\n", a, b, r1, x1, r2)
   @printf("大円の直径が 10 寸のとき，小円の直径は %g 寸である\n", 10r2/r1)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :red)
   circle(r2, a - 3r2, r2, :green)
   circle(r2, a - r2, r2, :green)
   circle(a - r2, a - r2, r2, :green)
   plot!([0, b, a], [0, a, 0], color=:blue, lw=0.5)
   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(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(r2, a - r2, "小円:r2,(r2,a-r2)", :green, :center, delta=-delta)
       point(r2, a - 3r2, "小円:r2,(r2,a-3r2)", :green, :center, delta=-delta)
       point(a - r2, a - r2, "小円:r2,(a-r2,a-r2)", :green, :center, delta=-delta)
       point(b, a, "(b,a)", :black, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;