算額（その1337）

算額（その1337）

新潟県小千谷市　小千谷二荒神社 天保4年(1833)
http://www.wasan.jp/niigata/niko.html

涌田和芳，外川一仁:小千谷二荒神社の紛失算額，長岡工業高等専門学校研究紀要，第 51 巻，p.35-40，2015.
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_51/51_35wakuta.pdf

キーワード：累円
#Julia, #SymPy, #算額, #和算

長方形の中に円弧，等円から始まる累円を入れる。長方形の長辺の長さが 28寸，等円の直径が 4 寸，累円の松園の直径が 1 寸のとき，（等円の次から末円までの）左右の円の個数はいくつか。

長方形の長辺の長さを 2a
円弧の半径と中心座標を R, (0, R)
等円の半径と中心座標を (x1, r1); x1 = a - r1
二円，三円の半径と中心座標を [r2, (x2, r2)], [r3, (x3, r3)], ...
末円の半径と中心座標を rn, (xn, rn)
とおく。

まず，円弧と等円，二円の関連について以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive,
     r1::positive, x1::posistive,
     r2::positive, x2::posistive;
x1 = a - r1
eq1 = x1^2 + (R - r1)^2 - (R + r1)^2
eq2 = x2^2 + (R - r2)^2 - (R + r2)^2
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (R, r2, x2))[2];

二円の半径は r1*(r1 - a)^2/(a + r1)^2 である。

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

   (-a + r1)^2/(4*r1)
   r1*(-a + r1)^2/(a + r1)^2
   (-a + r1)^2/(a + r1)

次に，R, r2, x2 が既知として，三円のパラメータ r3, x3 を求める。

@syms r3::positive, x3::posistive;
(R, r2, x2) = (x1^2/(4*r1), r1*x1^2/(2*r1 + x1)^2, x1^2/(2*r1 + x1))
eq4 = x3^2 + (R - r3)^2 - (R + r3)^2
eq5 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res2 = solve([eq4, eq5], (r3, x3))[2];

res2[1] |> simplify |> factor |> println
res2[2] |> simplify |> factor |> println

   r1*(-a + r1)^2/(a + 3*r1)^2
   (-a + r1)^2/(a + 3*r1)

三円の半径は r1*(r1 - a)^2/(a + 3*r1)^2 である。

さらに，r3, x3 も既知として，四円のパラメータ r4, x4 を求める。

@syms r4::positive, x4::posistive;
(R, r2, x2) = (x1^2/(4*r1), r1*x1^2/(2*r1 + x1)^2, x1^2/(2*r1 + x1))
(r3, x3) = (r1*(-a + r1)^2/(a + 3*r1)^2, (-a + r1)^2/(a + 3*r1))
eq6 = x4^2 + (R - r4)^2 - (R + r4)^2
eq7 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res3 = solve([eq6, eq7], (r4, x4))[2];

res3[1] |> simplify |> factor |> println
res3[2] |> simplify |> factor |> println

   r1*(-a + r1)^2/(a + 5*r1)^2
   (-a + r1)^2/(a + 5*r1)

ここまでやってくると，累円の半径は規則的な式で表されていることがわかる。
r1*(r1 - a)^2/(a + 1*r1)^2
r1*(r1 - a)^2/(a + 3*r1)^2
r1*(r1 - a)^2/(a + 5*r1)^2

なお，各段階での累円の半径は，r1, a のみで表されることがわかり，四円の半径は 0.5 であることもわかる。

res3[1](r1 => 4//2, a => 28//2) |> println

   1/2

もう少し考察を進めていく予定であったが，答えが出てしまった。
つまり，四円の直径がちょうど 1 寸であり，これが末円である。
累円は，二円，三円，四円の 3 個，左右合わせて 6 個というのが解である。

r2 = 1.125, x = 9, 直径 = 2.25
r3 = 0.72, x = 7.2, 直径 = 1.44
r4 = 0.5, x = 6, 直径 = 1
r5 = 0.367347, x = 5.14286, 直径 = 0.734694
r6 = 0.28125, x = 4.5, 直径 = 0.5625
r7 = 0.222222, x = 4, 直径 = 0.444444
r8 = 0.18, x = 3.6, 直径 = 0.36
r9 = 0.14876, x = 3.27273, 直径 = 0.297521
r10 = 0.125, x = 3, 直径 = 0.25
r11 = 0.106509, x = 2.76923, 直径 = 0.213018

以下では追加的に一般解を求める。等円の次の円から順に一番目，二番目，...，n 番目とすれば，n 番目の円の
半径は r1*(r1 - a)^2/(a + （2n - 1)*r1)^2
中心の x 座標は (r1 - a)^2/(a + (2n - 1)*r1)
である。

等円の直径が 4 寸，長方形の長辺が 28 寸のとき，「半径 - 0.5」の式を n で解くと n = 3 が得られる。

@syms n, rn
rn = r1*(r1 - a)^2/(a + (2n - 1)*r1)^2
rn = rn(r1 => 4//2, a => 28//2)
rn |> println

   288/(4*n + 12)^2

eq = rn - 1//2
eq |> println

   -1/2 + 288/(4*n + 12)^2

solve(eq, n)[2] |> println

   3

function draw(a, r1, more)
   pyplot(size=(700, 700), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   names = ["二", "三", "四"]
   (R, r2, x2) = ((-a + r1)^2/(4*r1), r1*(-a + r1)^2/(a + r1)^2, (a^2 - 2*a*r1 + r1^2)/(a + r1))
   x1 = a - r1
   θ = atand(-sqrt(R^2 - a^2), a)
   plot()
   rect(-a, 0, a, R - sqrt(R^2 - a^2), :green)
   circlef(0, R, R, :cornsilk2, beginangle=180 - θ, endangle=360 + θ)
   circle2f(x1, r1, r1)
   for n = 1:100
       i = 2n - 1
       r = r1*(r1 - a)^2/(a + i*r1)^2
       x = (r1 - a)^2/(a + i*r1)
       n <= 10 && @printf("r%d = %g, x = %g, 直径 = %g\n", n + 1, r, x, 2r)
       circle2f(x, r, r, i)
       n <= 3 && point(x, r, names[n], :white, :center, :vcenter, mark=false)
   end
   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(x1, r1, "等円:r1\n(x1,r1)", :white, :center, delta=-3delta)
       point(a, R - sqrt(R^2 - a^2), "(a,R-√(R^2-a^2))", :green, :right, :bottom, delta=delta/2)
   end  
end;

draw(28/2, 4/2, true)