長野県飯山市静間大久保 静間観音堂（静観庵）弘化5年(1848)
中村信弥「改訂増補 長野県の算額」県内の算額(P.177)
http://www.wasan.jp/zoho/zoho.html
日堂薫，宮崎雄也，鷲森勇希，伊藤栄一：「飯山市静間観音堂の算額」の復元
https://sbc21.co.jp/gakkokagaku/2015/27.pdf
キーワード：累円，正方形，斜線
#Julia #SymPy #算額 #和算 #数学
正方形の中に，四分円，半円（大円），天円，累円（甲円，乙円，丙円 ...）があり，天円，累円は隣り合う円と外接しさらに正方形の一辺と半円にも外接している。さらに，全円，中円，小円と斜線があり互いに接している（斜線と接しているのは全円と右側の小円のみ）。
累円の個数 n を，全円，中円，小円の直径と n 番目の累円の直径であらわせ。
注：半円と小円は外接していない。
正方形の一辺の長さを \(a\)，斜線と正方形の辺の交点座標を \( (0,\ b),\ (c,\ 0)\)
全円の半径と中心座標を \(r_1,\ (r_1,\ r_1)\)
中円の半径と中心座標を \(r_2,\ (r_2,\ r_2),\ (3r_2 - 2r_3,\ r_2),\ (5r_2 - 4r_3,\ r_2)\)
小円の半径と中心座標を \(r_3,\ (2r_2 - r_3,\ r_2),\ (4r_2 - 3r_3,\ r_2)\)
とおく。
順次以下の手順にしたがって，解を求める。
まず，以下の連立方程式を解き \(a,\ b,\ c\) を求める（\(r_1,\ r_2,\ r_3\) を含む式で表される）。
include("julia-source.txt"); # julia-source.txt ソース
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]
(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))
正方形の一辺の長さは \(4r_1(r_2 - r_3)/(r_1 - r_2)\) である。
次に，大円の半径と中心座標を \(r_4,\ (x_4,\ a - r_4)\) とおき，以下の連立方程式を解く。
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))[1] # 1 of 2
(a/4, a*(1 - sqrt(2)/2))
天円の半径は，正方形の一辺の長さの 1/4 である。
次に累円間の関係式を求める。
天円を初円 \( (x_4,\ a - r_4)\) とし，次の甲円 \( (x_5,\ a - r_5)\) を決める。
@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))[1] # 1 of 2
( (-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))
この解は，左側にある円に接する円のパラメータなので，これにより求められた円とその右側にある円とも同じ関係が成り立つ。
以下の関数は現在の累円から次の累円を求める関数である。
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\) 番目の累円の半径 \(r_n\) については
\(1/\sqrt{r_n} = n/\sqrt{大円の半径} + 1/\sqrt{天円の半径}\)
が成り立つ。
大円，天の半径は上述の通り，\(a/2,\ a/4\) なので
\(1/\sqrt{r_n} = n/\sqrt{a/2} + 1/\sqrt{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}\cdot (\sqrt{r_1(r_2 - r_3)/(r_1 - r_2)/(r_n)} - 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)
@printf("n = %d; n = %d; r = %.15g; x = %.15g\n", round(Int, n), n, r4, 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;
以下のアイコンをクリックして応援してください
