算額（その824）

算額（その824）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に 2 本の斜線と，圭（二等辺三角形），甲円 3 個，乙円 2 個を入れる。乙円の直径が 3 寸のとき，甲円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 3r1)
乙円の半径と中心座標を r2, (x2, R - 2r1 +r2)
斜線と円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, x0, y0
eq1 = x2^2 + (R - 2r1 +r2)^2 - (R - r2)^2
eq2 = dist(x0, y0, 0, -R, x1, R - 3r1) - r1^2
eq3 = dist(x0, y0, 0, -R, 0, R - r1) - r1^2
eq4 = dist(x0, y0, 0, -R, x2, R - 2r1 + r2) - r2^2
eq5 = dist(sqrt(R^2 - (R - 2r1)^2), R - 2r1, 0, -R, x1, R - 3r1) - r1^2
eq6 = x0^2 + y0^2 - R^2;

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, r1, x1, x2, x0, y0) = u
   return [
       x2^2 - (R - r2)^2 + (R - 2*r1 + r2)^2,  # eq1
       -r1^2 + (-x0 + x0*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2) + x1)^2 + (R - 3*r1 - y0 - (-R - y0)*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq2
       -r1^2 + (x0*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2) - x0)^2 + (R - r1 - y0 - (-R - y0)*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq3
       -r2^2 + (-x0 + x0*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2) + x2)^2 + (R - 2*r1 + r2 - y0 - (-R - y0)*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq4
       -r1^2 + (-r1 - (-2*R + 2*r1)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2))^2 + (x1 + sqrt(R^2 - (R - 2*r1)^2)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2) - sqrt(R^2 - (R - 2*r1)^2))^2,  # eq5
       -R^2 + x0^2 + y0^2,  # eq6
   ]
end;

r2 = 3/2
iniv = BigFloat[7.9, 2.1, 3.5, 3.5, 2.3, 7.7]
res = nls(H, ini=iniv)

   ([8.0, 2.0, 3.4641016151377544, 3.4641016151377544, 2.262270442538942, 7.673469387755102], true)

乙円の直径が 3 寸のとき，甲円の直径は 4 寸である。

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

R = 8;  r1 = 2;  x1 = 3.4641;  x2 = 3.4641;  x0 = 2.26227;  y0 = 7.67347

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3//2
   (R, r1, x1, x2, x0, y0) = res[1]
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  x1 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", R, r1, x1, x2, x0, y0)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r1, r1)
   circle2(x2, R - 2r1 + r2, r2, :blue)
   circle2(x1, R - 3r1, r1)
   plot!([-x0, 0, x0], [y0, -R, y0], color=:magenta, lw=0.5)
   y00 = R - 2r1
   x00 = sqrt(R^2 - y00^2)
   plot!([-x00, 0, x00, -x00], [y00, -R, y00, y00], color=:magenta, lw=0.5)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R - r1, " 甲円:r1\n(0,R-r1)", :red, :center, :bottom, delta=delta/2)
       point(x1, R - 3r1, " 甲円:r1\n(x1,R-3r1)", :red, :center, delta=-delta/2)
       point(x2, R - 2r1 + r2, " 乙円:r2\n(x2,R-2r1+r2)", :blue, :center, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :magenta, :left, :bottom, delta=delta/2)
       point(sqrt(R^2 - (R - 2r1)^2), R - 2r1, "sqrt(R^2-(R-2r1)^2),R-2r1 ", :black, :right, delta=-delta/2)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
   end
end;