算額（その707）

算額（その707）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

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

キーワード：円5個，外円，弦，斜線2本
#Julia, #SymPy, #算額, #和算

直径が 2.5 寸の外円内に弦と斜線を入れ，区画された領域に等円を 4 個入れる。等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
弦と外円の交点座標を (x1, y1); y1 = 2r - R, x1 = sqrt(R^2 - y1^2)
斜線と外円の交点座標を (x2, y2); x2 = sqrt(R^2 - y2^2)
等円の半径と中心座標を r, (0, r - R), (0, R - r), (x, y1 + r)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     x1::positive, y1::positive,
     x2::positive, y2::positive
y1 = 2r - R
x1 = sqrt(R^2 - y1^2)
x2 = sqrt(R^2 - y2^2)
eq1 = dist(0, y1, x2, y2, 0, R - r) - r^2
eq2 = dist(0, y1, x2, y2, x, y1 + r) - r^2
eq3 = x^2 + (y1 + r)^2 - (R - r)^2;
# res = solve([eq1, eq2, eq3], (r, x, y2))

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

function H(u)
   (r, x, y2) = u
   return [
       -r^2 + (2*R - 3*r)^2*(R^2 - y2^2)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2)^2 + (2*R - 3*r - (2*R - 3*r)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq1
       -r^2 + (r - (r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))*(R - 2*r + y2)/(R^2 - y2^2 + (R - 2*r + y2)^2))^2 + (x - sqrt(R^2 - y2^2)*(r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq2
       x^2 + (-R + 3*r)^2 - (R - r)^2,  # eq3
   ]
end;

R = 25//20
iniv = BigFloat[0.5, 0.8, 1.1]
res = nls(H, ini=iniv)

   ([0.4734152343522919, 0.7576940667754868, 1.0587960931595892], true)

外円の直径が 2.5 寸のとき，等円の直径は 2*0.4734152343522919 = 0.9468304687045838 寸である。

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

r = 0.473415;  x = 0.757694;  x1 = 1.21268;  y1 = -0.30317;  x2 = 0.664418;  y2 = 1.0588

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25//20
   (r, x, y2) = res[1]
   y1 = 2r - R
   x1 = sqrt(R^2 - y1^2)
   x2 = sqrt(R^2 - y2^2)
   @printf("等円の直径 = %g\n", 2r)
   @printf("r = %g;  x = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n", r, x, x1, y1, x2, y2)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, r - R, r)
   circle(x, y1 + r, r)
   circle(-x, y1 + r, r)
   circle(0, R - r, r)
   segment(-x1, y1, x1, y1, :green)
   segment(0, y1, x2, y2, :green)
   segment(0, y1, -x2, y2, :green)
   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(x, y1 + r, " (x,y1+r)", :red, :left, :vcenter)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(0, r - R, " r-R", :red, :left, :vcenter)
       point(0, y1, " y1", :green, :left, delta=-delta/2)
       point(x1, y1, "(x1,y1) ", :green, :right, delta=-delta/2)
       point(x2, y2, "(x2y2)", :green, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;