九十六 岩手県大船渡市立根町 気仙安養寺稲荷堂 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：累円，外円
#Julia #SymPy #算額 #和算 #数学
全円の中に甲円と心円，および乙円から始まる累円（丙円，丁円，戊円，己円...）が入っている。甲円の直径は全円の半径に等しい。心円の直径が 8 寸 9 分のとき，己円の直径を求めよ。
全円の半径と中心座標を \(R,\ (0,\ 0);\ R = 2r_1\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_1)\)
心円の半径と中心座標を \(r_2,\ (0,\ -r_2)\)
乙円の半径と中心座標を \(r_3,\ (x_3,\ y_3);\ x_3 = r_3,\ y_3 < 0\)
丙円の半径と中心座標を \(r_4,\ (x_4,\ y_4)\)
とおく。
まず，甲円の半径 \(r_1\) と乙円のパラメータ \(r_3,\ x_3,\ y_3\) を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive,
r3::positive, y3::negative,
r2::positive
R = 2r1
eq1 = r3^2 + (r1 - y3)^2 - (r1 + r3)^2
eq2 = r3^2 + (r2 + y3)^2 - (r3 + r2)^2
eq3 = r3^2 + y3^2 - (R - r3)^2
res = solve([eq1, eq2, eq3], (r1, r3, y3))[1]
(7*r2, 56*r2/9, -14*r2/3)
累円は，「ある円に外接する，半径が小さい次の円」の連続である。
乙円と丙円の関係は全円と甲円 の半径により，以下の連立方程式で得られる漸化式で計算される。
@syms R::positive, r1::positive,
r3::positive, y3::negative,
r2::positive
R = 2r1
@syms r4, x4, y4, x3
#x3 = r3
#(r1, r2, y3) = (7*r2, 56*r2/9, -14*r2/3)
eq4 = x4^2 + (r1 - y4)^2 - (r1 + r4)^2
eq5 = x4^2 + y4^2 - (R - r4)^2
eq6 = (x4 - x3)^2 + (y4 - y3)^2 - (r3 + r4)^2
res = solve([eq4, eq5, eq6], (r4, x4, y4))[1]
( (8*r1^3 + 4*r1^2*r3 - 20*r1^2*y3 - 2*r1*r3^2 - 4*r1*r3*y3 + 10*r1*x3^2 + 14*r1*y3^2 - r3^3 + 3*r3^2*y3 + r3*x3^2 + r3*y3^2 - 3*x3^2*y3 - 2*sqrt(2)*x3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 3*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)), (8*r1^2*x3 - 4*r1*r3*x3 - 4*r1*x3*y3 + 2*sqrt(2)*r1*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 4*r3^2*x3 + sqrt(2)*r3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) + 4*x3^3 + 4*x3*y3^2 - 3*sqrt(2)*y3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2), 3*sqrt(2)*x3*sqrt(-(-4*r1^2 - 4*r1*r3 - r3^2 + x3^2 + y3^2)*(2*r1*r3 - 2*r1*y3 - r3^2 + x3^2 + y3^2))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2) + (-8*r1^3 + 4*r1^2*r3 + 12*r1^2*y3 + 10*r1*r3^2 - 12*r1*r3*y3 + 2*r1*x3^2 - 6*r1*y3^2 + 3*r3^3 - 9*r3^2*y3 - 3*r3*x3^2 - 3*r3*y3^2 + 9*x3^2*y3 + 9*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)))
これにより得られる式は，丙円と丁円，丁円と己円...にも同じように適用される。以下の関数は，ある円 \(r_3,\ x_3,\ y_3\) の次の円の 半径，中心の \(x,\ y\) 座標を返す。
これを次々に適用することにより，累円のパラメータを求めることができる。
nextcircle(r1, r3, x3, y3) = ( (8*r1^3 + 4*r1^2*r3 - 20*r1^2*y3 - 2*r1*r3^2 - 4*r1*r3*y3 + 10*r1*x3^2 + 14*r1*y3^2 - r3^3 + 3*r3^2*y3 + r3*x3^2 + r3*y3^2 - 3*x3^2*y3 - 2*sqrt(2)*x3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 3*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)), (8*r1^2*x3 - 4*r1*r3*x3 - 4*r1*x3*y3 + 2*sqrt(2)*r1*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 4*r3^2*x3 + sqrt(2)*r3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) + 4*x3^3 + 4*x3*y3^2 - 3*sqrt(2)*y3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2), 3*sqrt(2)*x3*sqrt(-(-4*r1^2 - 4*r1*r3 - r3^2 + x3^2 + y3^2)*(2*r1*r3 - 2*r1*y3 - r3^2 + x3^2 + y3^2))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2) + (-8*r1^3 + 4*r1^2*r3 + 12*r1^2*y3 + 10*r1*r3^2 - 12*r1*r3*y3 + 2*r1*x3^2 - 6*r1*y3^2 + 3*r3^3 - 9*r3^2*y3 - 3*r3*x3^2 - 3*r3*y3^2 + 9*x3^2*y3 + 9*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)));
甲円の初期パラメータ \( (r_3,\ x_3,\ y_3) = (27.68888888888889,\ 27.68888888888889,\ -20.76666666666667)\) をはじめとして，累円のパラメータを実際に求めてみる。
nextcircle(31.15, 14.658823529411768, 43.976470588235316, 18.3235294117647)
(7.55151515151515, 37.75757575757576, 39.645454545454534)
nextcircle(31.15, 7.55151515151515, 37.75757575757576, 39.645454545454534)
(4.371929824561412, 30.603508771929828, 49.18421052631578)
プログラムで書いてみると以下のようになる。
r2 = 8.9/2
(r1, r3, y3) = (7*r2, 56*r2/9, -14*r2/3)
R = 2r1
(r, x, y) = (r3, r3, y3)
for i = 1:10
println( (r, x, y))
(r, x, y) = nextcircle(31.15, r, x, y)
end
(27.68888888888889, 27.68888888888889, -20.76666666666667)
(14.658823529411762, 43.976470588235294, 18.32352941176471)
(7.551515151515154, 37.75757575757576, 39.64545454545455)
(4.371929824561407, 30.60350877192981, 49.18421052631579)
(2.7999999999999976, 25.199999999999978, 53.899999999999984)
(1.931782945736425, 21.24961240310076, 56.5046511627907)
(1.4079096045197808, 18.302824858757038, 58.07627118644067)
(1.0695278969957025, 16.042918454935606, 59.09141630901287)
(0.8390572390572382, 14.263973063973097, 59.782828282828284)
(0.6753387533875355, 12.83143631436317, 60.2739837398374)
己円の半径は 2.7999999999999976 である（直径は 5.6 寸）。
なお，累円の半径は分数で表すことができ，分子は 1246 で，分母は \(g(n) = 20n^2 - 20n + 45\) である。
using Printf
g(n) = 20n^2 - 20n + 45;
for n = 1:10
@printf("n = %2d, 分母 = %4d, 半径 = %.15g\n", n, g(n), 1246/g(n))
end
\(n = 1,\ 分母 = 45,\ 半径 = 27.6888888888889\)
\(n = 2,\ 分母 = 85,\ 半径 = 14.6588235294118\)
\(n = 3,\ 分母 = 165,\ 半径 = 7.55151515151515\)
\(n = 4,\ 分母 = 285,\ 半径 = 4.3719298245614\)
\(n = 5,\ 分母 = 445,\ 半径 = 2.8\)
\(n = 6,\ 分母 = 645,\ 半径 = 1.93178294573643\)
\(n = 7,\ 分母 = 885,\ 半径 = 1.40790960451977\)
\(n = 8,\ 分母 = 1165,\ 半径 = 1.06952789699571\)
\(n = 9,\ 分母 = 1485,\ 半径 = 0.839057239057239\)
\(n = 10,\ 分母 = 1845,\ 半径 = 0.675338753387534\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
names = Char["甲乙丙丁戊己庚辛壬癸"...]
r2 = 8.9/2
(r1, r3, y3) = (7*r2, 56*r2/9, -14*r2/3)
R = 2r1
@printf("心円の直径が %g のとき，己円の直径は %g である。\n", 2r2, 2r3)
plot(showaxis=false)
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
circlef(0, 0, R, :orange)
circlef(0, r1, r1, :green)
circlef(0, -r2, r2, :purple)
(rnow, xnow, ynow) = (r3, r3, y3)
for i = 1:9
@printf("%s: R = %.15g; r1 = %.15g; r = %.15g; x = %.15g; y = %.15g\n", names[i+1], R, r1, rnow, xnow, ynow)
circlef(xnow, ynow, rnow, i)
circlef(-xnow, ynow, rnow, i)
point(xnow, ynow, names[i+1], :white, :center, :vcenter, mark=false)
(rnow, xnow, ynow) = nextcircle(31.15, rnow, xnow, ynow)
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(0, r1, "甲", :white, :center, :vcenter, mark=false)
point(0, -r2, "心", :white, :center, :vcenter, mark=false)
end
end;
心円の直径が 8.9 のとき，己円の直径は 55.3778 である。
乙: R = 62.3; r1 = 31.15; r = 27.6888888888889; x = 27.6888888888889; y = -20.7666666666667
丙: R = 62.3; r1 = 31.15; r = 14.6588235294118; x = 43.9764705882353; y = 18.3235294117647
丁: R = 62.3; r1 = 31.15; r = 7.55151515151515; x = 37.7575757575758; y = 39.6454545454545
戊: R = 62.3; r1 = 31.15; r = 4.37192982456141; x = 30.6035087719298; y = 49.1842105263158
己: R = 62.3; r1 = 31.15; r = 2.8; x = 25.2; y = 53.9
庚: R = 62.3; r1 = 31.15; r = 1.93178294573643; x = 21.2496124031008; y = 56.5046511627907
辛: R = 62.3; r1 = 31.15; r = 1.40790960451978; x = 18.302824858757; y = 58.0762711864407
壬: R = 62.3; r1 = 31.15; r = 1.0695278969957; x = 16.0429184549356; y = 59.0914163090129
癸: R = 62.3; r1 = 31.15; r = 0.839057239057238; x = 14.2639730639731; y = 59.7828282828283
デカルトの円定理を用いる方法
円の直径（半径）だけを求めるならば，デカルトの円定理を使って，累円の半径を求めてゆくことができる。
半径 \(r_1, prev, nxt\) の 3 円が内接する外円の半径を \(R\) とすると，以下の式が成り立つ。
using SymPy
@syms R, r1, prev, nxt
nxtr = -1/(1/r1 + 1/prev + 1/nxt - 2sqrt(1/(r1*prev) + 1/(prev*nxt) + 1/(nxt*r1))) - R
これを解いて \(nxt\) を求める式を得る。
result = solve(nxtr, nxt)[1]
result |> println
(R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2)
これを関数にする。この関数は，円の半径を入力すれば，その円に外接する次の円の半径を返す。
res2(prev, R, r1) = (R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2)
res2 (generic function with 1 method)
\(R, r_1, r_3\) が求まった段階で，次々に関数を適用してゆく。
r2 = 8.9/2
r1 = 7*r2
R = 2r1
r3 = r2 * 56/9 # 乙円
(r2, R, r1, r3)
(4.45, 62.300000000000004, 31.150000000000002, 27.68888888888889)
nxt = res2(r3, R, r1) # 丙円
14.65882352941177
nxt = res2(nxt, R, r1) # 丁円
7.551515151515152
nxt = res2(nxt, R, r1) # 戊円
4.371929824561404
nxt = res2(nxt, R, r1) # 己円
2.8000000000000003
nxt = res2(nxt, R, r1) # 庚円
1.9317829457364344
以下のアイコンをクリックして応援してください
