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

橘田彌曾八元克天明八年戊申2月(1788)

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：円7個，外円
#Julia #SymPy #算額 #和算 #数学

大円の中に，甲，乙，丙，丁，戊，己の 6 個の円が入っている。丁円，戊円，己円の直径がそれぞれ 872.3 寸，671 寸，572 寸のとき，大円の直径はいかほどか。

（必要ならば）己円の中心が \(y\) 軸上に来るように回転する。中心の \(x\) 座標が 0 になる。
大円の半径と中心座標を \(r_0, (0, 0)\)
甲円の半径と中心座標を \(r_1, (x_1, y_1)\)
乙円の半径と中心座標を \(r_2, (x_2, y_2)\)
丙円の半径と中心座標を \(r_3, (x_3, y_3)\)
丁円の半径と中心座標を \(r_4, (x_4, y_4); r_4 = 872.3/2\)
戊円の半径と中心座標を \(r_5, (x_5, y_5); r_5 = 671/2\)
己円の半径と中心座標を \(r_6, (x_6, y_6); r_6 = 572/2, x_6 = 0\)
として以下の連立方程式の数値解を求める。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms r0,
r1, x1, y1, r2, x2, y2, r3, x3, y3,
r4, x4, y4, r5, x5, y5, r6, x6, y6

eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = x3^2 + y3^2 - (r0 - r3)^2
eq4 = x4^2 + y4^2 - (r0 - r4)^2
eq5 = x5^2 + y5^2 - (r0 - r5)^2
eq6 = x6^2 + y6^2 - (r0 - r6)^2
eq7 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq8 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq9 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq10 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq11 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq12 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq13 = (x2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq14 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq15 = (x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2;

function H(u)
(r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, x4, y4, x5, y5, y6) = u
return [
x1^2 + y1^2 - (r0 - r1)^2, # eq1
x2^2 + y2^2 - (r0 - r2)^2, # eq2
x3^2 + y3^2 - (r0 - r3)^2, # eq3
x4^2 + y4^2 - (r0 - r4)^2, # eq4
x5^2 + y5^2 - (r0 - r5)^2, # eq5
x6^2 + y6^2 - (r0 - r6)^2, # eq6
-(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2, # eq7
-(r1 + r3)^2 + (x1 - x3)^2 + (y1 - y3)^2, # eq8
-(r1 + r4)^2 + (x1 - x4)^2 + (y1 - y4)^2, # eq9
-(r1 + r5)^2 + (x1 - x5)^2 + (y1 - y5)^2, # eq10
-(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2, # eq11
-(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2, # eq12
-(r2 + r6)^2 + (x2 - x6)^2 + (y2 - y6)^2, # eq13
-(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2, # eq14
-(r3 + r6)^2 + (x3 - x6)^2 + (y3 - y6)^2, # eq15
]
end;

(r4, r5, r6) = (872.3, 671, 572) ./2
x6 = 0
iniv = [big"95.0", 57, 12, 40, 41, -40, -43, 34, 36, -47, -72, 18, 75, -8, -79] .* (620.3/28)
res = nls(H, ini=iniv);
println([round(Float64(x), digits=6) for x in res[1]], " 収束:", res[2]);

[1586.0, 830.761905, 265.84381, 706.902857, 726.916667, -668.763333, -539.24, 623.071429, 672.917143, -688.777143, -994.422, 577.304, 1248.06, 78.08, -1300.0] 収束:true

\(r_0 = 1586\)
\(r_1 = 830.762; x_1 = 265.844; y_1 = 706.903\)
\(r_2 = 726.917; x_2 = -668.763; y_2 = -539.24\)
\(r_3 = 623.071; x_3 = 672.917; y_3 = -688.777\)
\(r_4 = 436.15; x_4 = -994.422; y_4 = 577.304\)
\(r_5 = 335.5; x_5 = 1248.06; y_5 = 78.08\)
\(r_6 = 286; x_6 = 0; y_6 = -1300\)

大円の半径は 1586 寸（直径は 3172 寸）である。

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

function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r4, r5, r6) = (872.3, 671, 572) ./2
x6 = 0
(r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, x4, y4, x5, y5, y6) = res[1]
@printf("r0 = %g; r1 = %g; x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; x3 = %g; y3 = %g; r4 = %g; x4 = %g; y4 = %g; r5 = %g; x5 = %g; y5 = %g; r6 = %g; x6 = %g; y6 = %g\n", r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6)
@printf("大円の半径 = %g; 直径 = %g\n", r0, 2r0)
plot()
circle(0, 0, r0, :gray)
circle(x1, y1, r1)
circle(x2, y2, r2, :blue)
circle(x3, y3, r3, :green)
circle(x4, y4, r4, :magenta)
circle(x5, y5, r5, :orange)
circle(x6, y6, r6, :brown)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(x1, y1, " 甲:r1,(x1,y1)", :red, :left, :vcenter)
point(x2, y2, "乙:r2,(x2,y2)", :blue, :center, :top, delta=-delta)
point(x3, y3, "丙:r3,(x3,y3)", :green, :center, :top, delta=-delta)
point(x4, y4, "丁:r4,(x4,y4)", :magenta, :center, :top, delta=-delta)
point(x5, y5, "戊:r5,(x5,y5)", :orange, :center, :bottom, delta=delta/2)
point(x6, y6, " 己:r6,(x6,y6)", :black, :left, delta=-delta/2)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;