『雑題』三十（大西佐兵衛）
愛媛和算研究会：『雑題』三十巻（大西佐兵衛著）を読むにあたって（その2），令和5年7月。
https://ehimewasan.com/wp-content/uploads/2024/09/72a03830eef2b8f8be223d0e0ff00b5e-1.pdf
キーワード：円7個
#Julia #SymPy #算額 #和算 #数学
雑題【11の4】日円，月円，火円，土円が互いに接している。その隙間に，水円，木円，金円を容れる。火円，水円，木円，金円の直径が 17.16 寸，3.9 寸，2.86 寸，4.29 寸のとき，土円の直径はいかほどか。
月円の中心を座標原点に置き，木円の中心が \(y\) 軸上になるように配置する。
日円の半径と中心座標を \(r_1,\ (x_1,\ y_1)\)
月円の半径と中心座標を \(r_2,\ (0,\ 0)\)
火円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
水円の半径と中心座標を \(r_4,\ (x_4,\ y_4)\)
木円の半径と中心座標を \(r_5,\ (0,\ r_2 + r_5)\)
金円の半径と中心座標を \(r_6,\ (x_6,\ y_6)\)
土円の半径と中心座標を \(r_7,\ (x_7,\ y_7)\)
とおき，以下の連立方程式の数値解を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1, x1, y1, r2, r3, x3, y3, r4, x4, y4,
r5, y5, r6, x6, y6, r7, x7, y7
y5 = r2 + r5
x5 = 0
eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + y3^2 - (r2 + r3)^2
eq3 = x7^2 + y7^2 - (r2 + r7)^2
eq4 = (x1 - x7)^2 + (y1 - y7)^2 - (r1 + r7)^2
eq5 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq6 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq7 = x4^2 + (y4 - y5)^2 - (r4 + r5)^2
eq8 = x6^2 + (y5 - y6)^2 - (r5 + r6)^2
eq9 = (x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2
eq10 = (x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2
eq11 = x1^2 + (y1 - y5)^2 - (r1 + r5)^2
eq12 = x4^2 + y4^2 - (r4 + r2)^2
eq13 = x6^2 + y6^2 - (r2 + r6)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13], (r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7))
function H(u)
(r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7) = u
return [
(x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2, # eq1
x3^2 + y3^2 - (r2 + r3)^2, # eq2
x7^2 + y7^2 - (r2 + r7)^2, # eq3
(x1 - x7)^2 + (y1 - y7)^2 - (r1 + r7)^2, # eq4
(x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2, # eq5
(x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2, # eq6
x4^2 + (y4 - (r2 + r5))^2 - (r4 + r5)^2, # eq7
x6^2 + (r2 + r5 - y6)^2 - (r5 + r6)^2, # eq8
(x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2, # eq9
(x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2, # eq10
x1^2 + (y1 - (r2 + r5))^2 - (r1 + r5)^2, # eq11
x4^2 + y4^2 - (r4 + r2)^2, # eq12
x6^2 + y6^2 - (r2 + r6)^2, # eq13
]
end;
(r3, r4, r5, r6) = (17.16, 3.9, 2.86, 4.29) ./ 2
iniv = BigFloat[14, 0.6, 22, 5, 13, 4, 3, 7, -4, 7, 13, -18, 3]
res = nls(H, ini=iniv)
([14.211137652235761, 0.6147384017903117, 22.463598083149616, 5.404545528578323, 13.496531066824165, 3.662125448852095, 3.3670139679392537, 6.538543941248187, -3.570333542168671, 6.651943782515173, 13.200000000000012, -18.404061703588194, 2.7238992520764294], true)
# r7: 土円の直径
2 * res[1][11]
26.400000000000023
土円の直径は 26.4 寸である。
術は以下のような式を示している。
(火径, 水径, 木径, 金径) = (17.16, 3.9, 2.86, 4.29)
天 = 火径*(水径 - 木径)
地 = 天*水径
人 = 天 - (金径 - 木径)*(火径 - 金径)*水径^2/金径^2
土 = 地/人
26.400000000000034
他の値を試しても数値解と同じ答えが得られる。
やはり，例示された図は与えられた数値（直径）を反映したものではないようだ。
描画関数プログラムのソースを見る
function draw(r3, r4, r5, r6, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7) =res[1]
# [14.211137652235761, 0.6147384017903117, 22.463598083149616, 5.404545528578323, 13.496531066824165, 3.662125448852095, 3.3670139679392537, 6.538543941248187, -3.570333542168671, 6.651943782515173, 13.200000000000012, -18.404061703588194, 2.7238992520764294]
x5 = 0
y5 = r2 + r5
@printf("火円 = %.15g, 水円 = %.15g, 木円 = %.15g, 金円 = %.15g のとき，土円 = %.15g\n", 2r3, 2r4, 2r5, 2r6, 2r7)
plot()
circle(x1, y1, r1)
circle(0, 0, r2, :blue)
circle(x3, y3, r3, :green)
circle(x3, y3, r3, :magenta)
circle(x4, y4, r4, :darkorange)
circle(x5, y5, r5, :purple)
circle(x6, y6, r6, :tomato)
circle(x7, y7, r7, :blueviolet)
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(x1, y1, "日円:r1,(x1,y1)", :red, :center, delta=-delta)
point(0, 0, "月円:r2,(0,0)", :blue, :center, delta=-delta)
point(x3, y3, "火円:r3,(x3,y3)", :magenta, :center, delta=-delta)
point(x4, y4, "水円:r4,(x4,y4)", :darkorange, :left, :vcenter, delta=delta, deltax=5delta)
point(x5, y5, "木円:r5,(x5,y5)", :purple, :center, delta=-4delta)
point(x6, y6, "金円:r6,(x6,y6)", :tomato, :right, :vcenter, deltax=-4delta)
point(x7, y7, "土円:r7,(x7,y7)", :blueviolet, :center, delta=-delta)
end
end;
draw(17.16/2, 3.9/2, 2.86/2, 4.29/2, true)
火円 = 17.16, 水円 = 3.9, 木円 = 2.86, 金円 = 4.29 のとき，土円 = 26.4
以下のアイコンをクリックして応援してください
