算額（その668）

算額（その668）

長野県飯山市静間大久保 静間観音堂（静観庵）弘化5年(1848)

中村信弥「改訂増補 長野県の算額」県内の算額（177）
http://www.wasan.jp/zoho/zoho.html

日堂薫，宮崎雄也，鷲森勇希，伊藤栄一：「飯山市静間観音堂の算額」の復元
https://sbc21.co.jp/gakkokagaku/2015/27.pdf

正方形の中に，四分円，半円（大円），天円，累円（甲円，乙円，丙円 ...）があり，天円，累円は隣り合う円と外接しさらに正方形の一辺と半円にも外接している。さらに，全円，中円，小円と斜線があり互いに接している（斜線と接しているのは全円と右側の小円のみ）。
累円の個数 n を，全円，中円，小円の直径と n 番目の累円の直径であらわせ。
注：半円と小円は外接していない。

正方形の一辺の長さを a，斜線と正方形の辺の交点座標を (0, b), (c, 0)
全円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (r2, r2), (3r2 - 2r3, r2), (5r2 - 4r3, r2)
小円の半径と中心座標を r3, (2r2 - r3, r2), (4r2 - 3r3, r2)
とおく。

順次以下の手順にしたがって，解を求める。

まず，以下の連立方程式を解き a, b, c を求める（r1, r2, r3 を含む式で表される）。

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

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, r3::positive,
     x4::positive, r4::positive
eq1 = r2/(a - (5r2 - 4r3)) - r1/(a - r1)
eq2 = r2/(a - c) - b/a
eq3 = a + b - sqrt(a^2 + b^2) - 2r1
res = solve([eq1, eq2, eq3], (a, b, c))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*r1*(r2 - r3)/(r1 - r2), r1*(r1 - 5*r2 + 4*r3)/(r1 - 3*r2 + 2*r3), 4*(r1^2*r2 - r1^2*r3 - 6*r1*r2^2 + 10*r1*r2*r3 - 4*r1*r3^2 + 3*r2^3 - 5*r2^2*r3 + 2*r2*r3^2)/(r1^2 - 6*r1*r2 + 4*r1*r3 + 5*r2^2 - 4*r2*r3))

正方形の一辺の長さは 4r1*(r2 - r3)/(r1 - r2) である。

次に，大円の半径と中心座標を r4, (x4, a - r4) とおき，以下の連立方程式を解く。

eq4 = (a - x4)^2 + r4^2 - (a - r4)^2
eq5 = (a - x4)^2 + (a/2 - r4)^2 - (a/2 + r4)^2
res2 = solve([eq4, eq5], (r4, x4))

   2-element Vector{Tuple{Sym, Sym}}:
    (a/4, a*(1 - sqrt(2)/2))
    (a/4, a*(sqrt(2)/2 + 1))

2 組の解が得られるが，最初のものが適解である。
天円の半径は，正方形の一辺の長さの 1/4 である。

次に累円間の関係式を求める。
天円を初円(x4, a - r4) とし，次の甲円 (x5, a - r5) を決める。

@syms r5, x5
eq6 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
eq7 = (a - x5)^2 + (a/2 - r5)^2 - (a/2 + r5)^2
solve([eq6, eq7], (r5, x5))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-a + x4)*(-2*sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) + a + 2*a*(-2*r4 + x4)/(-a + 2*r4) + x4)/(2*(-a + 2*r4)), sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) - a*(-2*r4 + x4)/(-a + 2*r4))
    ((-a + x4)*(2*sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) + a + 2*a*(-2*r4 + x4)/(-a + 2*r4) + x4)/(2*(-a + 2*r4)), -sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) - a*(-2*r4 + x4)/(-a + 2*r4))

2 組の解が得られるが，最初のものが適解である。
この解は，左側にある円に接する円のパラメータなので，これにより求められた円とその右側にある円とも同じ関係が成り立つ。

以下の関数は現在の累円から次の累円を求める関数である。

nextcircle(a, r4, x4) = ((x4 - a)*(-2sqrt(2a*r4)*(x4 - a)/(2r4 - a) + a + 2a*(x4 - 2r4)/(2r4 - a) + x4)/(2(2r4 - a)),
       sqrt(2a*r4)*(x4 - a)/(2r4 - a) - a*(x4 - 2*r4)/(2r4 - a))

ただし，このまま長々しい数式のままでの取り扱いは難しい。

直線の上に互いに接する3個の円の半径についてはデカルトの円定理より，以下のような関係式が成り立つ。

大円・天円・甲円の半径においては，1/√甲円の半径 = 1/√大円の半径 + 1/√天円の半径
大円・甲円・乙円の半径においては，1/√乙円の半径 = 1/√大円の半径 + 1/√甲円の半径
大円・乙円・丙円の半径においては，1/√丙円の半径 = 1/√大円の半径 + 1/√乙円の半径
大円・丙円・丁円の半径においては，1/√丁円の半径 = 1/√大円の半径 + 1/√丙円の半径
　：
下から順次代入していくと，甲円から数えて4番目の丁円の半径について
1/√丁円の半径 = 1/√大円の半径 + (1/√大円の半径 + 1/√乙円の半径)
   = 2/√大円の半径 + 1/√乙円の半径 = 2/√大円の半径 + (1/√大円の半径 + 1/√甲円の半径)
   = 3/√大円の半径 + 1/√甲円の半径 = 3/√大円の半径 + (1/√大円の半径 + 1/√天円の半径)
   = 4/√大円の半径 + 1/√天円の半径
甲円から数えて n 番目の累円の半径 rn については
1/√(rn) = n/√大円の半径 + 1/√天円の半径
が成り立つ。

大円，天の半径は上述の通り，a/2, a/4 なので
1/√(rn) = n/√(a/2) + 1/√(a/4)

@syms rn, n, r1, r2, r3
a = 4*r1*(r2 - r3)/(r1 - r2)
eq = 1/√rn - (n/√(a/2) + 1/√(a/4))

方程式を解いて n を求める。

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

   sqrt(2)*(-sqrt(rn) + sqrt(r1*(r2 - r3)/(r1 - r2)))/sqrt(rn)

SymPyではこれ以上簡約化できないが，以下のように簡約化できる。
sqrt(2)*(sqrt(r1*(r2 - r3)/(r1 - r2)/(rn)) - 1)

実行例を見れば，19 番目の累円の半径は 0.015874124835590503，でこのとき n は 19 となる。

(r1, r2, r3) = (4.3, 3, 1) ./ 2
rn = 0.015874124835590503
sqrt(2)*(sqrt(r1*(r2 - r3)/(r1 - r2)/(rn)) - 1)

   19.000000000000043

   r1 = 2.15;  r2 = 1.5;  r3 = 0.5;  a = 13.2308;  b = 5.33519;  c = 9.51091
   大円の半径 = 6.615384615384616
   天円の半径 = 3.307692307692308
   -19.619822485207095
   n = 0;  n = 0.0;  r = 3.307692307692308;  x = 3.875202587377986
   n = 1;  n = 1.0000000000000009;  r = 1.1350205593713572;  x = 7.750405174755974
   n = 2;  n = 1.9999999999999998;  r = 0.5675102796856791;  x = 9.355566643391244
   n = 3;  n = 2.9999999999999996;  r = 0.33950675324251667;  x = 10.233458601408486
   n = 4;  n = 3.999999999999998;  r = 0.22567545882547532;  x = 10.787058971033915
   n = 5;  n = 4.999999999999996;  r = 0.16079341811242404;  x = 11.168042584375158
          :
   n = 18;  n = 18.00000000000003;  r = 0.01755155072579467;  x = 12.549270122486426
   n = 19;  n = 19.000000000000043;  r = 0.015874124835590503;  x = 12.5826536817502

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (4.3, 3, 1) ./ 2
   (a, b, c) = (4*r1*(r2 - r3)/(r1 - r2), r1*(r1 - 5*r2 + 4*r3)/(r1 - 3*r2 + 2*r3), 4*(r1^2*r2 - r1^2*r3 - 6*r1*r2^2 + 10*r1*r2*r3 - 4*r1*r3^2 + 3*r2^3 - 5*r2^2*r3 + 2*r2*r3^2)/(r1^2 - 6*r1*r2 + 4*r1*r3 + 5*r2^2 - 4*r2*r3))
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  a = %g;  b = %g;  c = %g\n",
       r1, r2, r3, a, b, c)
   println("大円の半径 = $(a/2)")
   println("天円の半径 = $(a/4)")
   println((a - (5r2 - 4r3))^2 + (a - r2)^2 - (a + r2)^2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(3r2 - 2r3, r2, r2, :blue)
   circle(5r2 - 4r3, r2, r2, :blue)
   circle(2r2 - r3, r2, r3, :green)
   circle(4r2 - 3r3, r2, r3, :green)
   segment(0, b, a, 0, :orange)
   circle(a, a, a, :magenta, beginangle=180, endangle=270)
   circle(a, a/2, a/2, :blue, beginangle=90, endangle=270)
   (r4, x4) = (a/4, a*(1 - sqrt(2)/2))
   for i = 1:20
       n = sqrt(2)*(-sqrt(r4) + sqrt(r1*(r2 - r3)/(r1 - r2)))/sqrt(r4)
       println("n = $(round(Int, n));  n = $n;  r = $r4;  x = $(x4)")
       circle(x4, a - r4, r4)
       (r4, x4) = nextcircle(a, r4, x4)
   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)
       point(a, a/2, "(a,a/2)")
       point(r1, r1, "全円:r1,(r1,r1)", :black, :center, :bottom, delta=delta)
       point(5r2 - 4r3, r2, "中円:r2,(5r2-4r3,r2)", :black, delta=-delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom)
       point(c, 0, " c", :black, :left, :bottom, delta=delta/2)
       segment(0, r2, c, r2, :gray80)
       segment(c, 0, c, r2, :gray80)
   end
end;