六十 群馬県佐波郡玉村町飯塚 光琳寺 天保10年(1837)
六十一 群馬県高崎市石原町 清水寺 天保10年(1837)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：累円，半円，台形
#Julia #SymPy #算額 #和算 #数学
等脚台形の中に半円，全円，次円，三円，四円，五円を容れる。上底（大頭）が 90 寸のとき，各円の直径はいかほどか。
上底，下底を \(2b, 2a\)
半円の半径と中心座標を \(r_0, (0, r_0); r_0 = 2r_1 = b/2\)
全円の半径と中心座標を \(r_1, (0, r_1)\)
次円の半径と中心座標を \(r_2, (x_{21}, r_2), (x_{22}, r_0 - r_2)\)
三円の半径と中心座標を \(r_3, (x_3, r_3)\)
四円の半径と中心座標を \(r_4, (x_4, r_4)\)
五円の半径と中心座標を \(r_5, (x_5, r_5)\)
とおき，以下の連立方程式を解く。
四円のパラメータまでは一度に求めることができるが，SymPy の性能上，五円のパラメータも一緒に求めることができないので，分けて解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, r0::positive, r1::positive,
r2::positive, x21::positive, x22::positive,
r3::positive, x3::positive,
r4::positive, x4::positive,
r5::positive, x5::positive;
r0 = 2r1
eq1 = x21^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = x22^2 + r2^2 - (r0 - r2)^2
eq3 = x22^2 + (r0 - r2 - r1)^2 - (r1 + r2)^2
eq4 = dist2(r0, r0, a, 0, x21, r2, r2)
eq5 = x3^2 + (r0 - r3)^2 - (r0 + r3)^2
eq6 = (x21 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq7 = x4^2 + (r0 - r4)^2 - (r0 + r4)^2
eq8 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r2, x21, x22, r3, x3, r4, x4, a))[4] # 4 of 6
(r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))
@syms r5::positive, x5::positive
(r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(Sym(2))*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(Sym(2))/2 + 2))
eq9 = x5^2 + (r0 - r5)^2 - (r0 + r5)^2
eq10 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
#res2 = solve([eq7, eq8], (r4, x4,))
res2 = solve([eq9, eq10], (r5, x5))[1] # 1 of 2
(2*r1/25, 4*r1/5)
次円（二円），三円，四円，五円の半径は 全円の半径の 1/2, 2/9, 1/8, 2/25 倍である。
\(n\) 番目の円の半径は，全円の半径の \(2/n^2\) である。
上底が 90 のとき，次円の直径 = 22.5，三円の直径 = 10，四円の直径 = 5.625，五円の直径 = 3.6 である。
全てのパラメータは以下のとおりである。
\(r_2 = 11.25; x_{21} = 45; x_{22} = 31.8198; r_3 = 5; x_3 = 30\)
\(r_4 = 2.8125; x_4 = 22.5; r_5 = 1.8; x_5 = 18; a = 60.9099\)
描画関数プログラムのソースを見る
function draw(b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r0 = b/2
r1 = r0/2
(r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))
(r5, x5) = (2*r1/25, 4*r1/5)
@printf("上底が %g のとき，次円の直径 = %g，三円の直径 = %g，四円の直径 = %g，五円の直径 = %g である。\n", b, 2r2, 2r3, 2r4, 2r5)
@printf("r2 = %g; x21 = %g; x22 = %g; r3 = %g; x3 = %g; r4 = %g; x4 = %g; r5 = %g; x5 = %g; a = %g\n", r2, x21, x22, r3, x3, r4, x4, r5, x5, a)
plot([a, r0, -r0, -a, a], [0, r0, r0, 0, 0], color=:magenta, lw=0.5)
circle(0, r0, r0, beginangle=180, endangle=360)
circle(0, r1, r1, :blue)
circle2(x21, r2, r2, :green)
circle2(x22, r0 - r2, r2, :green)
circle2(x3, r3, r3, :orange)
circle2(x4, r4, r4, :brown)
circle2(x5, r5, r5, :tomato)
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(b/2, b/2, "(b/2,b/2)", :magenta, :right, :bottom, delta=delta)
point(a, 0, "a", :magenta, :center, delta=-4delta)
point(0, r1, "全円:r1,(0,r1)", :blue, :center, delta=-delta)
point(x21, r2, "次円:r2\n(x21,r2)", :green, :center, :bottom, delta=delta)
point(x22, r0 - r2, "次円:r2\n(x22,r0-r2)", :green, :center, :bottom, delta=delta)
point(x3, r3, "三円:r3,(x3,r3)", :orange, :left, delta=-18delta)
point(x4, r4, "四円:r4,(x4,r4)", :brown, :left, delta=-22delta)
point(x5, r5, "五円:r5,(x5,r5)", :tomato, :left, delta=-28delta)
end
end;
draw(90, true)
以下のアイコンをクリックして応援してください
