算額（その643）

算額（その643）

長野県軽井沢町峠 熊野神社 安政4年(1857)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円の中に 4 つの円弧と 2 つの等円が入っている。等円は互いに外接し円弧とも外接している。
大矢，小矢がそれぞれ 49寸，36.75寸のとき，等円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
等円の半径と中心座標を r1, (r1, y1)
上部の円弧を構成する円の半径と中心座標を r2, (r2, r0)
下部の円弧を構成する円の半径と中心座標を r3, (r3, -r0)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms r0::positive, r1::positive, y1::positive, r2::positive, r3::positive

eq1 = (r2 - r1)^2 + (r0 - y1)^2 - (r1 + r2)^2
eq2 = r2^2 + r0^2 - (r0 + r2 - 3675//100)^2
eq3 = r3^2 + r0^2 - (r0 + r3 - 49)^2
eq4 = r0 - sqrt(r2*r3)
eq5 = 1/sqrt(r2) + 1/sqrt(r3) - 1/sqrt(r1);

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)
   (r0, r1, y1, r2, r3) = u
   return [
       (r0 - y1)^2 + (-r1 + r2)^2 - (r1 + r2)^2,  # eq1
       r0^2 + r2^2 - (r0 + r2 - 147/4)^2,  # eq2
       r0^2 + r3^2 - (r0 + r3 - 49)^2,  # eq3
       r0 - sqrt(r2)*sqrt(r3),  # eq4
       1/sqrt(r3) + 1/sqrt(r2) - 1/sqrt(r1),  # eq5
   ]
end;

iniv = BigFloat[73.5, 20, 10, 54, 97]
res = nls(H, ini=iniv)

   (BigFloat[73.49999999999999999999999999999999999999999999999999999999999999999999999999889, 17.99999999999999999999999999999999999999999999999999999999999999999999511910947, 10.50000000000000000000000000000000000000000000000000000000000000000000678816348, 55.12500000000000000000000000000000000000000000000000000000000000000000000000055, 97.99999999999999999999999999999999999999999999999999999999999999999999999999779], true)

等円の半径は 18 寸である（直径 36寸 ）。
その他のパラメータは以下の通りである。

r0 = 73.5;  r1 = 18;  y1 = 10.5;  r2 = 55.125;  r3 = 98

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, y1, r2, r3) = res[1]
   @printf("r0 = %g;  r1 = %g;  y1 = %g;  r2 = %g;  r3 = %g\n", r0, r1, y1, r2, r3)
   plot()
   circle(0, 0, r0, :blue)
   circle(r1, y1, r1)
   circle(-r1, y1, r1)
   circle(r2, r0, r2, :green)
   circle(-r2, r0, r2, :green)
   circle(r3, -r0, r3, :magenta)
   circle(-r3, -r0, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(r1, y1, "(r1,y1)",)
       point(r2, r0, "(r2,r0)")
       point(r3, -r0, "(r3,-r0) ", :magenta, :right, :bottom, delta=delta/2)
       point(r0, 0, " r0", :blue, :left, :bottom, delta=delta/2)
       fr3 = Float64(r3)
       fr0 = Float64(r0)
       rect(-fr3, -fr0, fr3, fr0)
       plot!(xlims=(-fr3 - 5, fr3 + 5), ylims=(-fr0 - 5, fr0 + 5))
       θ = atand(r0/r2)
       (px1, py1, px2, py2) = (r0*cosd(θ), r0*sind(θ), 2(r0 - r2)*cosd(θ), 2(r0 - r2)*sind(θ))
       segment(0, 0, r2, r0, :gray60)
       segment(px1, py1, px2, py2, lw=1)
       segment(-px1, py1, -px2, py2, lw=1)
       point((px1 + px2)/2, (py1 + py2)/2, " 小矢", mark=false)
       θ = atand(-r0/r3)
       (px1, py1, px2, py2) = (r0*cosd(θ), r0*sind(θ), -(r0 - r3)*cosd(θ), -(r0 - r3)*sind(θ))
       segment(0, 0, r3, -r0, :gray60)
       segment(px1, py1, px2, py2, lw=1)
       segment(-px1, py1, -px2, py2, lw=1)
       point((px1 + px2)/2, (py1 + py2)/2, " 大矢", mark=false, delta=2delta)
   end
end;

算額（その642）

算額（その642）

長野県飯田市八幡町 松尾八幡神社（鳩ヶ嶺八幡宮） 明治6年(1873)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円に二等辺三角形が内接し，金円，銀円が入っている。
外円の直径が 10 寸のとき，金円，銀円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
金円の半径と中心座標を r1, (0, R - r1)
銀円の半径と中心座標を r2, (r2, R - 2r1 -r2), (x2, y2)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms R::positive, r1::positive, r2::positive, 弦::positive, x::positive, x::negative, x2::positive, y2::negative

y = 2R - 2r1
x = sqrt(R^2 - (R - 2r1)^2)
弦 = sqrt(y^2 + x^2)
eq1 = y + sqrt(R^2 - (R - 2r1)^2) - 弦 - 2r2
eq2 = sqrt(R^2 - 弦^2/4) + 2r2 - R |> simplify  # -R + 2*r2 + sqrt(R*r1)
eq3 = y2/x2 * y/x + 1
eq4 = 2sqrt(x2^2 + (y2 + R)^2 - r2^2) - 弦;
res = solve([eq1, eq2, eq3, eq4], (r1, r2, x2, y2))

   2-element Vector{NTuple{4, Sym}}:
    (R*(4 - sqrt(7))/8, R*(4 - sqrt(8 - 2*sqrt(7)))/8, R*(6 - sqrt(-16*sqrt(2)*sqrt(4 - sqrt(7)) - 4*sqrt(14)*sqrt(4 - sqrt(7)) + 8*sqrt(7) + 41))/16, R*(-24 - sqrt(-112*sqrt(2)*sqrt(4 - sqrt(7)) - 28*sqrt(14)*sqrt(4 - sqrt(7)) + 56*sqrt(7) + 287) + 6*sqrt(7) + 4*sqrt(-16*sqrt(2)*sqrt(4 - sqrt(7)) - 4*sqrt(14)*sqrt(4 - sqrt(7)) + 8*sqrt(7) + 41))/48)
    (R*(4 - sqrt(7))/8, R*(4 - sqrt(8 - 2*sqrt(7)))/8, R*(sqrt(-16*sqrt(2)*sqrt(4 - sqrt(7)) - 4*sqrt(14)*sqrt(4 - sqrt(7)) + 8*sqrt(7) + 41) + 6)/16, R*(-24 - 4*sqrt(-16*sqrt(2)*sqrt(4 - sqrt(7)) - 4*sqrt(14)*sqrt(4 - sqrt(7)) + 8*sqrt(7) + 41) + sqrt(-112*sqrt(2)*sqrt(4 - sqrt(7)) - 28*sqrt(14)*sqrt(4 - sqrt(7)) + 56*sqrt(7) + 287) + 6*sqrt(7))/48)

2 組の解が得られるが，2 番目のものが適解である。また，それぞれのパラメータの式は簡約化できる。

res[2][1] |> println
res[2][2] |> sympy.sqrtdenest |> simplify |> println
res[2][3] |> sympy.sqrtdenest |>  simplify |> println
res[2][4] |> sympy.sqrtdenest |>  simplify |> println

   R*(4 - sqrt(7))/8
   R*(5 - sqrt(7))/8
   R*(5 + 2*sqrt(7))/16
   -R*(2 + sqrt(7))/16

金円，銀円の直径は外円の直径の (4 - √7)/8 倍と (5 - √7)/8 倍である。
外円の直径が 10 寸のとき，金円，銀円の直径は 1.6928108611692616寸と 2.942810861169262寸である。

10*(4 - √7)/8, 10*(5 - √7)/8

   (1.6928108611692616, 2.942810861169262)

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

   r1 = 0.846405;  r2 = 1.47141;  x2 = 3.21609;  y2 = -1.4518

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 10//2
   (r1, r2, x2, y2) = (R/16) .* (8 - 2√7, 10 - 2√7, 5 + 2√7, -2 - √7)
   y1 = R - 2r1
   x1 = sqrt(R^2 - y1^2)
   弦 = sqrt((2R - 2r1)^2 + (R^2 - y1^2))
   @printf("金円の直径 = %g;  銀円の直径 = %g\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r1, r2, x2, y2)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(r2, R - 2r1 - r2, r2, :green)
   circle(-r2, R - 2r1 - r2, r2, :green)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   plot!([0, x1, -x1, 0], [-R, y1, y1, -R], color=:black, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r1, " 金円:r1,(0,R-r1)", :blue, :left, :vcenter)
       point(0, R - 2r1, " R-2r1", :blue, :left, delta=-delta/2)
       point(r2, R - 2r1 - r2, "銀円:r2\n(r2,R-2r1-r2)", :green, :center, delta=-delta/2)
       point(x2, y2, "銀円:r2\n(x2,y2)", :green, :center, delta=-delta/2)
   end
end;