算額（その656）

算額（その656）

群馬の和算 111(1/2) 倉賀野神社 慶応3年
http://takasakiwasan.web.fc2.com/gunnsann/g111_1.html

正三角形内に最も大きな円弧を置き，その円弧と2辺に接する甲円を置く。甲円と円弧と底辺に接する乙円を置き，同様に順次丙円，丁円...を置く。
これらの直径を求めよ。

正三角形の一辺の長さを a とおき，甲円，乙円...の半径と中心座標を (r1, x1), (r2, x2)...とする。

甲円のパラメータ (r1, x1) が決まったあと，乙円のパラメータは (r1, x1) で記述できる。その後も，丙円，丁円と同じ関係式でパラメータが決まる。

まず最初に，正三角形内の最大弧を決定する。弧を形成する円は正三角形の 2 つの辺に頂点で接するものである。その中心座標は (a, a/√3) であることは簡単な計算で分る。この中心座標を (x, y) とおき，甲円の r1, x1 を以下の連立方程式を解いて求める。

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

using SymPy
@syms R::positive, a::positive,
     r1::positive, x1::positive,
     r2::positive, x2::positive
a = 8
R = a/sqrt(Sym(3))
(x, y) = (a, R)
eq1 = (x - x1)^2 + (y - r1)^2 - (R + r1)^2
eq2 = r1/x1 - y/x
res = solve([eq1, eq2], (r1, x1))

   2-element Vector{Tuple{Sym, Sym}}:
    (8*sqrt(3)/9, 8/3)
    (8*sqrt(3), 24)

2 組の解が得られるが，2 番目のものが適解である。
すなわち，r1 = 8*sqrt(3)/9, x1 = 8/3 である。

次に，乙円のパラメータ r2, x2 は，r1, x1 を未知数のままにして（r2, x2 が r1, x1 を含む式で表される）以下の連立方程式で求めることができる。

eq3 = (x - x2)^2 + (y - r2)^2 - (R + r2)^2
eq4 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2;

簡単な二元連立方程式であるが，SymPy の能力的に solve([eq3, eq4], (r2, x2)) は有限の時間内には解けないようなので，手動で解く。

まず，eq4 から x2 を求める。x2 は r2 が決まれば求めることができる。

# x2 を求める式
solve(eq4, x2) |> println

   Sym[-2*sqrt(r1)*sqrt(r2) + x1, 2*sqrt(r1)*sqrt(r2) + x1]

2 組の解が得られるが，2 番目のものが適解である。 すなわち，x2 = 2*sqrt(r1)*sqrt(r2) + x1 である。この式は算額において有名な式で，つまる所は，直線の上に載っている，互いに接している円の中心間の水平距離は 2 つの円の積の平方根の 2 倍である，つまり，x2 - x1 = 2√(r1*r2) にほかならない。

得られた x2 を eq3 に代入すると以下の方程式が得られる。

eq3 = (x - (2sqrt(r1)*sqrt(r2) + x1))^2 + (y - r2)^2 - (R + r2)^2
eq3 |> println

   (-r2 + 8*sqrt(3)/3)^2 - (r2 + 8*sqrt(3)/3)^2 + (-2*sqrt(r1)*sqrt(r2) - x1 + 8)^2

この方程式から r2 を求める。

res = solve(eq3, r2)
res |> println

   Sym[(x1 - 8)^2*(4*sqrt(2)*3^(3/4)*sqrt(r1)*(-3*r1^2 + 16*sqrt(3)*r1 - 64) + (3*r1 + 8*sqrt(3))*(3*r1^2 - 16*sqrt(3)*r1 + 64))/(4*(3*r1^2 - 16*sqrt(3)*r1 + 64)^2), (x1 - 8)^2*(4*sqrt(2)*3^(3/4)*sqrt(r1) + 3*r1 + 8*sqrt(3))/(4*(3*r1^2 - 16*sqrt(3)*r1 + 64))]

2 組の解が得られるが，最初のものが適解である。 簡約化すると以下のようになる。

# r2 を求める式
r2 = res[1] |> simplify |> println

   (x1 - 8)^2*(-4*sqrt(2)*3^(3/4)*sqrt(r1) + 3*r1 + 8*sqrt(3))/(4*(3*r1^2 - 16*sqrt(3)*r1 + 64))

以上で求めた r2, x2 を求める式は r1, x1 だけを含む。また，この式は互いに接する円の間に成り立つので，たとえば丙円と丁円の間においても成り立つ。そこで，以下のような関数としてまとめておく。

function nextcircle(r1, x1)
   # # r2 = (x1 - 8)^2*(-4*sqrt(Sym(2))*3^Sym(3//4)*sqrt(r1) + 3*r1 + 8*sqrt(Sym(3)))/(4*(3*r1^2 - 16*sqrt(Sym(3))*r1 + 64))
   r2 = (x1 - 8)^2*(-4*sqrt(2)*3^(3/4)*sqrt(r1) + 3*r1 + 8*sqrt(3))/(4*(3*r1^2 - 16*sqrt(3)*r1 + 64))
   x2 = 2*sqrt(r1)*sqrt(r2) + x1
   return (r2, x2)
end;

この関数を使って，累円のパラメータを求めてみよう。

(r1, x1) = (8*sqrt(3)/9, 8/3)

   (1.5396007178390019, 2.6666666666666665)

nextcircle(1.5396007178390019, 2.6666666666666665)

   (0.6188021535170058, 4.618802153517006)

nextcircle(0.6188021535170058, 4.618802153517006)

   (0.3316150746190428, 5.524791385931975)

nextcircle(0.3316150746190428, 5.524791385931975)

   (0.20626738450566878, 6.04786451314966)

   甲円の直径 = 3.0792014;  中心の X 座標 = 2.6666667
   乙円の直径 = 1.2376043;  中心の X 座標 = 4.6188022
   丙円の直径 = 0.6632301;  中心の X 座標 = 5.5247914
   丁円の直径 = 0.4125348;  中心の X 座標 = 6.0478645
   戊円の直径 = 0.2811508;  中心の X 座標 = 6.3884294
   己円の直径 = 0.2038283;  中心の X 座標 = 6.6278172
   庚円の直径 = 0.1545148;  中心の X 座標 = 6.8052841
   辛円の直径 = 0.1211510;  中心の X 座標 = 6.9421037
   壬円の直径 = 0.0975328;  中心の X 座標 = 7.0508060
   癸円の直径 = 0.0802036;  中心の X 座標 = 7.1392508

術では以下のような解が述べられている。関係式は簡潔であるが，直径（半径）の情報しか得られない。上述した式は半径と中心座標を求めることができる解である。

極 = 8/sqrt(0.75)
甲 = 極/3; 甲 |> println
子 = 甲*極; 乙 = 子/(2√子 + 極 + 甲); 乙 |> println
丑 = 乙*極; 丙 = 丑/(2√丑 + 極 + 乙); 丙 |> println

   3.0792014356780046
   1.2376043070340124
   0.6632301492380859

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   names = Char["甲乙丙丁戊己庚辛壬癸"...]
   n = length(names)
   rs = Vector{Float64}(undef, n)
   xs = Vector{Float64}(undef, n)
   a = 8
   R = a/sqrt(Sym(3))
   (x, y) = (a, R)
   plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:blue, lw=0.5)
   circle(x, y, R, :green)
   (r1, x1) = (8*sqrt(3)/9, 8/3)
   circle(x1, r1, r1)
   for i = 1:n
       @printf("%s円の直径 = %.7f;  中心の X 座標 = %.7f\n", names[i], 2r1, x1)
       (rs[i], xs[i]) = (r1, x1)
       (r2, x2) = nextcircle(r1, x1)
       circle(x2, r2, r2)
       (r1, x1) = (r2, x2)
   end
   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)
       for i = 1:3
           point(xs[i], rs[i], names[i], :red, :center, :vcenter, mark=false)
       end
   end
end;