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

安島直円：不朽算法（上巻）

東北大学総合知デジタルアーカイブ
https://touda.tohoku.ac.jp/portal/item/10010000022977
キーワード：円6個
#Julia #SymPy #算額 #和算 #数学

第三十　大円，中円，小円が互いに接しており，その中にできる隙間に，互いに外接する甲円，乙円，丙円が入っている。大円，中円，小円の直径がそれぞれ 6 寸，5 寸，4 寸のとき，甲円，乙円，丙円の直径はいかほどか。

大円の半径と中心座標を \(r_1,\ (0,\ 0)\)
中円の半径と中心座標を \(r_2,\ (x_2,\ y_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)\)
とおき，以下の連立方程式の数値解を求める。

答では甲円径は 60/119 寸，乙円径は 30/61 寸，丙円径は 15/31 寸となっているが，乙円と丙円の直径が入れ替わっている。

すなわち，甲円径は 60/119 寸であるが，乙円径は 15/31 寸，丙円径は 30/61 寸である。

include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1, r2, x2, y2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6
eq1 = x4^2 + y4^2 - (r1 + r4)^2
eq2 = x6^2 + y6^2 - (r1 + r6)^2
eq3 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq4 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq5 = (r1 + r3 - x5)^2 + y5^2 - (r3 + r5)^2
eq6 = (r1 + r3 - x6)^2 + + y6^2 - (r3 + r6)^2
eq7 = x2^2 + y2^2 - (r1 + r2)^2
eq8 = (r1 + r3 - x2)^2 + y2^2 - (r2 + r3)^2
eq9 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq10 = (x4 - x6)^2 + (y4 - y6)^2 - (r4 + r6)^2
eq11 = (x5 - x6)^2 + (y5 -y6)^2 - (r5 + r6)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (x2, y2, r4, x4, y4, r5, x5, y5, r6, x6, y6))

println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
println(eq9, ", # eq9")
println(eq10, ", # eq10")
println(eq11, ", # eq11")

x4^2 + y4^2 - (r1 + r4)^2, # eq1
x6^2 + y6^2 - (r1 + r6)^2, # eq2
-(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2, # eq3
-(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2, # eq4
y5^2 - (r3 + r5)^2 + (r1 + r3 - x5)^2, # eq5
y6^2 - (r3 + r6)^2 + (r1 + r3 - x6)^2, # eq6
x2^2 + y2^2 - (r1 + r2)^2, # eq7
y2^2 - (r2 + r3)^2 + (r1 + r3 - x2)^2, # eq8
-(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2, # eq9
-(r4 + r6)^2 + (x4 - x6)^2 + (y4 - y6)^2, # eq10
-(r5 + r6)^2 + (x5 - x6)^2 + (y5 - y6)^2, # eq11

function H(u)
(x2, y2, r4, x4, y4, r5, x5, y5, r6, x6, y6) = u
return [
x4^2 + y4^2 - (r1 + r4)^2, # eq1
x6^2 + y6^2 - (r1 + r6)^2, # eq2
-(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2, # eq3
-(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2, # eq4
y5^2 - (r3 + r5)^2 + (r1 + r3 - x5)^2, # eq5
y6^2 - (r3 + r6)^2 + (r1 + r3 - x6)^2, # eq6
x2^2 + y2^2 - (r1 + r2)^2, # eq7
y2^2 - (r2 + r3)^2 + (r1 + r3 - x2)^2, # eq8
-(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2, # eq9
-(r4 + r6)^2 + (x4 - x6)^2 + (y4 - y6)^2, # eq10
-(r5 + r6)^2 + (x5 - x6)^2 + (y5 - y6)^2, # eq11
]
end;
r1 = 6/2
r2 = 5/2
r3 = 4/2
iniv = BigFloat[3.5, 4.23,
0.5, 2.9, 1.6,
0.5, 3.3, 1.5,
0.5, 3.0, 1.1,
]
res = nls(H, ini=iniv)

([3.5, 4.242640687119285, 0.25210084033613445, 2.8487394957983194, 1.5687074809516686, 0.24193548387096775, 3.338709677419355, 1.505453147042327, 0.2459016393442623, 3.0491803278688523, 1.1128237867853863], true)

# 甲径
res[1][3]*2, 60/119

(0.5042016806722689, 0.5042016806722689)

# 乙径
res[1][6]*2, 15/31

(0.4838709677419355, 0.4838709677419355)

# 丙径
res[1][9]*2, 30/61

(0.4918032786885246, 0.4918032786885246)

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

function draw(r1, r2, r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x2, y2, r4, x4, y4, r5, x5, y5, r6, x6, y6) = res[1]
@printf("大円径 = %g; 中円径 = %g; 小円径 = %g; 甲円径 = %.15g; 乙円径 = %.15g; 丙円径 = %.15g\n",
2r1, 2r2, 2r3, 2r4, 2r5, 2r6)
plot()
circle(0, 0, r1)
circle(x2, y2, r2, :blue)
circle(r1 + r3, 0, r3, :magenta)
circle(x4, y4, r4, :orange)
circle(x5, y5, r5, :green)
circle(x6, y6, r6, :black)
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, "大:r1,(0,0)", :red, :center, delta=-delta)
point(x2, y2, "中:r2,(x2,y2)", :blue, :center, delta=-delta)
point(r1 + r3, 0, "小:r3,(r1+r3,0)", :magenta, :center, delta=-delta)
point(x4, y4, "甲:r4,(x4,y4)", :orange, :center, :bottom, delta=3delta, deltax=4delta)
point(x5, y5, "乙:r5,(x5,y5)", :green, :left, delta=-delta/2, deltax=2delta)
point(x6, y6, "丙:r6,(x6,y6)", :black, :right, :vcenter, deltax=-2.5delta)
end
end
draw(6/2, 5/2, 4/2, true)