算額(その97)
東京都台東区下谷池端 東淵寺 明治2年(1869)
http://www.wasan.jp/tokyo/tobuti.html
正方形の中に半円,大円,小円が図のように配置されている。大円の径が二寸四分八厘のとき,小円の径はいくつか。
大円の半径を 248 として,図のように記号を定めて方程式を解く。
using SymPy
@vars R, r, x, y (positive=true, real=true);
@syms R::positive, r::positive, x::positive, y::positive;
R =248
eq1 = x^2 + (y - r)^2 - 4r^2
eq2 = x^2 + y^2 - R^2
eq3 = x^2 + (R + 2r - y)^2 - (R + r)^2;
(r, x, y) = solve([eq1, eq2, eq3], (r, x, y))[1]
(-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108, 248*sqrt(6)*sqrt(-(109 + 27*sqrt(17))^(5/3)*(909*sqrt(17) + 3755) - 2*(109 + 27*sqrt(17))^(4/3)*(4627 + 1125*sqrt(17)) + 32*(109 + 27*sqrt(17))*(2943*sqrt(17) + 12137))/(9*(109 + 27*sqrt(17))^(3/2)), -1984/(81*(109/729 + sqrt(17)/27)^(1/3)) + 248/9 + 248*(109/729 + sqrt(17)/27)^(1/3))
小円の半径は,100.00743003769321 である。元の単位で径に直せばほぼ 1 寸ということで,算額の答えにあっている。
-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108
100.00743003769321
using Plots
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
(r, x, y) = (-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108, 248*sqrt(6)*sqrt(-(109 + 27*sqrt(17))^(5/3)*(909*sqrt(17) + 3755) - 2*(109 + 27*sqrt(17))^(4/3)*(4627 + 1125*sqrt(17)) + 32*(109 + 27*sqrt(17))*(2943*sqrt(17) + 12137))/(9*(109 + 27*sqrt(17))^(3/2)), -1984/(81*(109/729 + sqrt(17)/27)^(1/3)) + 248/9 + 248*(109/729 + sqrt(17)/27)^(1/3))
R = 248
println("r = $r; x = $x; y = $y")
plot()
circle(0,2r + R, R, :green)
circle(0, r, r)
circle(x, y, r)
circle(-x, y, r)
circle(0, 0, R + r, :blue, beginangle=0, endangle=180)
segment(-R-r, 0, R+r, 0, :blue)
l = r + R
plot!([-l, l, l, -l, -l], [0, 0, 2l, 2l, 0], color=:black, lw=0.5)
if more
point(0, 2r + R, "2r+R ", :green, :right)
point(0, r, "r ", :red, :right)
point(x, y, "(x,y)", :red, :top)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;
r = 100.00743003769321; x = 191.57807116330426; y = 157.48600778909824