算額(その581)
長崎市 鎮西大社諏訪神社 明治20年(1887)
米光丁: 長崎県の和算の概説
http://hyonemitsu.web.fc2.com/Nagasakiwasan.pdf
問題 4. 外円の中に大小の円を入れる。外円と大円の直径が 32 寸,18 寸のとき,小円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y - r1)
小円の半径と中心座標を r2,(x2, y - r2), (0, y + r2)
弦と y 軸の交点座標を (0, y)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, y::positive, r1::positive, x1::negative,
r2::positive, x2::positive
eq1 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x1^2 + (y - r1)^2 - (R - r1)^2
eq3 = x2^2 + (y - r2)^2 - (R - r2)^2
eq4 = y + 2r2 - R
res = solve([eq1, eq2, eq3, eq4], (y, x1, r2, x2))
1-element Vector{NTuple{4, Sym}}:
((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))
小円の半径を表す式は長いが,問で与えられた数数値を代入すると,小円の半径は6,直径は 12 寸である。
res[1][3](R => 32//2, r1 => 18//2) |> println
6
その他のパラメータは以下の通りである。
R = 16; r1 = 9; x1 = -4.89898; r2 = 6; x2 = 9.79796; y = 4
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1) = (32, 18) .// 2
(y, x1, r2, x2) = ((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))
@printf("小円の直径 = %g; R = %g; r1 = %g; x1 = %g; r2 = %g; x2 = %g; y = %g\n",
2r2, R, r1, x1, r2, x2, y)
plot()
circle(0, 0, R, :magenta)
circle(x1, y - r1, r1)
circle(x2, y - r2, r2, :blue)
circle(0, y + r2, r2, :blue)
x = sqrt(R^2 - y^2)
segment(-x, y, x, y)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(x1, y - r1, "大円:r1,(x1,y-r1)", :red, :center, :top, delta=-delta/2)
point(x2, y - r2, "小円:r2,(x2,y-r2)", :blue, :center, :top, delta=-delta/2)
point(0, y + r2, " 小円:r2\n (x2,y+r2)", :blue, :left, :vcenter)
point(0, y, " y", :black, :left, :top, delta=-delta/2)
point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
end
end;