算額（その97）

算額（その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