算額(その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;
※コメント投稿者のブログIDはブログ作成者のみに通知されます