算額（その2404） - 算額あれこれ

東都浅草浅草寺境内稲荷社文化7年(1810)

内藤忠辰校：賽祠神算
Kano Knowledge Portal
https://kanokp.lib.u-tokyo.ac.jp/item/tohoku-10010000022299
キーワード：円13個，外円
#Julia #SymPy #算額 #和算 #数学

外円の中に，甲，乙，丙，丁，戊，己，庚，合わせて 12 個の円を容れる。甲円の直径が 6 寸，乙円の直径が 3 寸のとき，庚円の直径はいかほどか。
注：甲円と乙円の大きさが等しい場合には「算額（その0350）」になる。

外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0, r_1 - R)\)
乙円の半径と中心座標を \(r_2,\ (0,\ R - r_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
丁円の半径と中心座標を \(r_4,\ (x_4,\ y_4)\)
戊円の半径と中心座標を \(r_5,\ (x_5,\ y_5)\)
己円の半径と中心座標を \(r_6,\ (x_6,\ y_6)\)
庚円の半径と中心座標を \(r_7,\ (x_7,\ y_7)\)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt"); # julia-source.txt ソース
function driver(r1, r2)
function H(u)
function parameters()
R = r1 + r2
eq1 = x3^2 + (r1 - R - y3)^2 - (r1 + r3)^2
eq2 = x7^2 + (r1 - R - y7)^2 - (r1 + r7)^2
eq3 = x6^2 + (R - r2 - y6)^2 - (r2 + r6)^2
eq4 = x7^2 + (R - r2 - y7)^2 - (r2 + r7)^2
eq5 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq6 = (x3 - x7)^2 + (y3 - y7)^2 - (r3 + r7)^2
eq7 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq8 = (x4 - x7)^2 + (y4 - y7)^2 - (r4 + r7)^2
eq9 = (x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2
eq10 = (x5 - x7)^2 + (y5 - y7)^2 - (r5 + r7)^2
eq11 = (x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2
eq12 = x3^2 + y3^2 - (R - r3)^2
eq13 = x4^2 + y4^2 - (R - r4)^2
eq14 = x5^2 + y5^2 - (R - r5)^2
eq15 = x6^2 + y6^2 - (R - r6)^2
return [eq1, eq2, eq3, eq4, eq5,
eq6, eq7, eq8, eq9, eq10,
eq11, eq12, eq13, eq14, eq15]
end;
(r3, x3, y3, r4, x4, y4, r5, x5, y5,
r6, x6, y6, r7, x7, y7) = u
return parameters()
end;
iniv = BigFloat[20, 70, -35, 12, 83, 15, 12, 75, 35,
16, 60, 50, 20, 50, 20]
iniv = BigFloat[0.39, 1.3, -0.41, 0.22, 1.5, 0.16, 0.21, 1.5, 0.58, 0.32, 1.1, 0.97, 0.40, 0.93, 0.27].*r1
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
(r1, r2) = (6/2, 3/2)
(r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) = driver(r1, r2)
(r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) |> println

(1.0, 3.4641016151377544, 0.5, 0.5, 3.4641016151377544, 2.0, 0.42857142857142855, 2.969229955832361, 2.7857142857142856, 0.6, 2.0784609690826525, 3.3, 0.9, 2.0784609690826525, 1.8)

# 庚円の直径: 2r7
2r7

1.8

甲円の直径が 6 寸，乙円の直径が 3 寸のとき，庚円の直径は 1.8 寸である。

術は以下のようになっており，同じ結果が得られる。

@syms 甲, 乙
子 = 甲 + 乙
丑 = 甲*乙*3
A = 4*子^2 - 丑
庚 = 子/A*丑

\(\displaystyle \frac{3 乙甲 \left(乙 + 甲\right)}{- 3 乙甲 + 4 \left(乙 + 甲\right)^{2}}\)

庚(甲 => 6, 乙 => 3)

\(\displaystyle \frac{9}{5}\)

庚(甲 => 6, 乙 => 3).evalf()

\(1.8\)

描画関数プログラムのソースを見る

function draw(r1, r2, more=true)
pyplot(size=(600, 600), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAexMincho")
(r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) = driver(r1, r2)
R = r1 + r2
plot()
circle(0, 0, R, :green)
circle(0, r1 - R, r1)
circle(0, R - r2, r2, :blue)
circle2(x3, y3, r3, :magenta)
circle2(x4, y4, r4, :purple)
circle2(x5, y5, r5, :black)
circle2(x6, y6, r6, :orange)
circle2(x7, y7, r7, :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(0, 0, "外円:R,(0,0)", :green, :center, delta=-delta)
point(0, r1 - R, "甲円:r1,(0,r1-R)", :red, :center, delta=-delta)
point(0, R - r2, "乙円:r2,(0,R-r2)", :blue, :center, delta=-delta)
point(x3, y3, "丙円:r3,(x3,y3)", :magenta, :center, delta=-delta)
point(x4, y4, "丁円:r4,(x4,y4)", :purple, :left, :bottom, deltax=2delta, delta=4delta)
point(x5, y5, " 戊円:r5,(x5,y5)", :black, :left, :bottom, delta=r5 + delta/2)
point(x6, y6, " 己円:r6,(x6,y6)", :orange, :left, :bottom, delta=5delta)
point(x7, y7, "庚円:r7,(x7,y7)", :tomato, :center, delta=-delta)
adjust_xlims!(0, right=5delta)
end
end;
(r1, r2) = (6/2, 3/2)
draw(r1, r2, true)