算額（その599）

算額（その599）

 

七六 北埼玉郡騎西町中種足 雷神社 明治5年(1872)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

東京都中央区 福徳神社 令和2年(2020)
http://www.wasan.jp/tokyo/hukutoku1.html

正三角形を 2 本の線で 4 つの部分に分け，その中に大小の円をそれぞれ 2 個入れる。大円の直径が 1 のとき，小円の直径はいくらか。

正三角形の一辺の長さを 2a とする。
線分と斜辺の交点座標を (±x, (a - x)√3) とする。
大円の半径と中心座標を r1, (0, y1), (0, r1)
小円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, x::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive

eq1 = distance(0, sqrt(Sym(3))a, a, 0, 0, y1) - r1^2
eq2 = distance(0, sqrt(Sym(3))a, a, 0, x2, y2) - r2^2
eq3 = distance(-x, (a-x)sqrt(Sym(3)), a, 0, 0, y1) - r1^2
eq4 = distance(-x, (a-x)sqrt(Sym(3)), a, 0, x2, y2) - r2^2
eq5 = distance(x, (a-x)sqrt(Sym(3)), -a, 0, 0, r1) - r1^2
eq6 = distance(x, (a-x)sqrt(Sym(3)), -a, 0, x2, y2) - r2^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6])

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, x, y1, r2, x2, y2) = u
   return [
3*a^2/4 - sqrt(3)*a*y1/2 - r1^2 + y1^2/4,  # eq1
-r2^2 + (-3*a + 3*x2 + sqrt(3)*y2)^2/16 + (-sqrt(3)*a + sqrt(3)*x2 + y2)^2/16,  # eq2
(3*a^4/4 - 3*a^3*x/2 - sqrt(3)*a^3*y1/2 - a^2*r1^2 + 3*a^2*x^2/4 + a^2*y1^2/4 + a*r1^2*x + sqrt(3)*a*x^2*y1/2 + a*x*y1^2/2 - r1^2*x^2 + x^2*y1^2/4)/(a^2 - a*x + x^2),  # eq3
-r2^2 + (a + x)^2*(-sqrt(3)*a^2 + sqrt(3)*a*x + sqrt(3)*a*x2 + a*y2 - sqrt(3)*x*x2 + x*y2)^2/(16*(a^2 - a*x + x^2)^2) + (-3*a^3 + 6*a^2*x + 3*a^2*x2 + sqrt(3)*a^2*y2 - 3*a*x^2 - 6*a*x*x2 + 3*x^2*x2 - sqrt(3)*x^2*y2)^2/(16*(a^2 - a*x + x^2)^2),  # eq4
(3*a^4 - 2*sqrt(3)*a^3*r1 - 6*a^3*x - 3*a^2*r1^2 + 3*a^2*x^2 + 6*a*r1^2*x + 2*sqrt(3)*a*r1*x^2 - 3*r1^2*x^2)/(4*(a^2 - a*x + x^2)),  # eq5
-r2^2 + (-sqrt(3)*a^3 - sqrt(3)*a^2*x2 + a^2*y2 + sqrt(3)*a*x^2 + 2*a*x*y2 + sqrt(3)*x^2*x2 + x^2*y2)^2/(16*(a^2 - a*x + x^2)^2) + (3*a^3 - 6*a^2*x + 3*a^2*x2 - sqrt(3)*a^2*y2 + 3*a*x^2 - 6*a*x*x2 + 3*x^2*x2 + sqrt(3)*x^2*y2)^2/(16*(a^2 - a*x + x^2)^2),  # eq6
   ]
end;

r1 = 1/2
iniv = BigFloat[1.6, 0.66, 1.74, 0.322, 0.56, 1.1] .* (2r1)
res = nls(H, ini=iniv)

   (BigFloat[1.573132184970985912579820468331052169654464272601005512809483077205430728737341, 0.6610365985079724914150555424853615260166663839137217848002503187911426334707547, 1.724744871391588873493812097618873490781894353477567544930381227134074172251389, 0.3224744871391588398986467839454922576529915481447599591870862696864391756755701, 0.5585421958697395716072717462463142074739205600259983862524414059477805245308117, 1.112372435695794479712354981538013592707116186365025213663076513709077532563648], true)

大円の直径が 1 のとき，小円の直径は 0.644949 である。

その他のパラメータは以下の通り。
小円の直径 = 0.644949;  r1 = 0.5;  a = 1.57313;  x = 0.661037;  y1 = 1.72474;  r2 = 0.322474;  x2 = 0.558542;  y2 = 1.11237

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   (a, x, y1, r2, x2, y2) = res[1]
   @printf("小円の直径 = %g;  r1 = %g;  a = %g;  x = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", 2r2, r1, a, x, y1, r2, x2, y2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:magenta, lw=0.5)
   circle(0, y1, r1)
   circle(0, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   segment(-a, 0, x, (a - x)√3)
   segment(a, 0, -x, (a - x)√3)
   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, y1, " 大円:r1\n (0,y1)", :red, :left, :vcenter)
       point(0, r1, " 大円:r1\n (0,r1)", :red, :left, :vcenter)
       point(x2, y2, "小円:r2\n(x2,y2)", :blue, :center, delta=-delta)
       point(a, 0, "a", :magenta, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :magenta, :left, :vcenter)
       point(x, √3(a-x), " (x,√3(a-x))", :magenta, :left, :vcenter)
   end
end;