算額（その709）

算額（その709）

八六 加須市多聞寺 愛宕神社 明治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

キーワード：円8個，菱形，弦
#Julia, #SymPy, #算額, #和算

外円内に弦を境界として，上部に甲円 4 個，乙円 2 個をいれる。下部には菱形と，菱形に内接する楕円を容れる。
乙円の径を▢寸，楕円の短径を一寸五分としたとき，楕円の長径はいかほどか。

注：欠損した一文字は「答」から推測すると「一」であろう。なお，「答」にも「長径▢寸七分五厘有奇」と欠損文字があるが，こちらは「二」であろう。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 2r1); x1 = √3r1
乙円の半径と中心座標を r2, (x2, ya + r2); ya = R - 4r1 = 2r1 - R
楕円の長半径と短半径を a, b
楕円と菱形の接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::positive,
     a::positive, b::positive,
     x0::positive, y0::positive

x1 = sqrt(Sym(3))r1
yb = r1 - R
xb = sqrt(R^2 - yb^2)
eq1 = x1^2 + (R - 2r1)^2 - (R - r1)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (x1 - x2)^2 + (2r1 - r2)^2 - (r1 + r2)^2
eq4 = x0^2/a^2 + (y0 - r1 + R)^2/b^2 - 1
eq5 = -b^2*x0/(a^2*(y0 - r1 + R)) + r1/xb
eq6 = (y0 - yb)/(xb - x0) - r1/xb;

図の上半分は R, r1, x2 がわかれば描ける。

res = solve([eq1, eq2, eq3], (R, r1, x2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-2*sqrt(2)*r2 + 3*r2, r2*(3 - 2*sqrt(2))/3, -2*sqrt(6)*r2/3 + 2*sqrt(3)*r2/3)
    (2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)

2 組の解が得られるが，2 番目のものが適解である。

res[2][1] |> factor |> println
res[2][2] |> factor |> println
res[2][3] |> factor |> println

   r2*(2*sqrt(2) + 3)
   r2*(2*sqrt(2) + 3)/3
   2*r2*(sqrt(3) + sqrt(6))/3

r2 = 1/2
(2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)

   (2.914213562373095, 0.9714045207910317, 1.3938468501173515)

下半分も描くには SymPy は力不足なので，数値解を求める。

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

r2 = 1//2
b = 15//20
iniv = BigFloat[2.91, 0.97, 1.39, 0.88, -1.36, 1.38]
res = nls(H, ini=iniv)

   ([2.914213562373095, 0.9714045207910317, 1.3938468501173518, 0.8773124760905445, -1.363750587600455, 1.3804469258418701], true)

楕円の長径は 2.7608938516837402 （長半径は 1.3804469258418701）である。
「答」では「長径は▢寸七分5厘有奇」となっている。

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

   R = 2.91421;  r1 = 0.971405;  x2 = 1.39385;  ya = -0.971405;  x0 = 0.877312;  y0 = -1.36375;  a = 1.38045

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, b) = (1/2, 1.5/2)
   (R, r1, x2, x0, y0, a) = res[1]
   x1 = √3r1
   ya = R - 4r1  # = 2r1 - R
   xa = sqrt(R^2 - ya^2)
   yb = r1 - R
   xb = sqrt(R^2 - yb^2)
   println("長径 = $(2a) （長半径:a = $a）")
   @printf("r2 = %g;  b = %g\n", r2, b)
   @printf("R = %g;  r1 = %g;  x2 = %g;  ya = %g;  x0 = %g;  y0 = %g;  a = %g\n",
       R, r1, x2, ya, x0, y0, a)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, R - r1, r1)
   circle(0, R - 3r1, r1)
   circle(x1, R - 2r1, r1)
   circle(-x1, R - 2r1, r1)
   circle(0, r1 - R, r1, :gray90)
   circle(x2, ya + r2, r2, :orange) 
   circle(-x2, ya + r2, r2, :orange) 
   segment(-xa, ya, xa, ya)
   plot!([xb, 0, -xb, 0, xb], [yb, ya, yb, -R, yb], color=:green, lw=0.5)
   ellipse(0, r1 - R, a, b)
   if more
       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, :left, :vcenter)
       point(0, R - 3r1, "", :red)
       point(x1, R - 2r1, "(x1,R-2r1)", :red, :center, delta=-delta/2)
       point(x2, ya + r2, " 乙円:r2\n (x2,ya+r2)", :black, :left, :vcenter)
       point(0, ya, " ya=R-4r1", :black, :center, :bottom, delta=delta)
       segment(-xb, yb, xb, yb)
       segment(0, ya, 0, -R)
       point(0, r1 - R, "r1-R ", :black, :right, :bottom, delta=delta)
       point(a, r1 - R, "(a,r1-R) ", :red, :right, :bottom, delta=delta)
       point(0, r1 - R + b, " r1-R+b", :red, :left, delta=-delta)
   end
end;