算額（その20）

算額（その20）

4 種類の円と 1 つの半円が図のように配置されている。それぞれの円の半径を求めよ。

福島県田村郡三春町御木沢 厳島神社 明治 18 年
http://www.wasan.jp/fukusima/miharuitukusima1.html

半円の半径を 1 とする。そのほかに必要とする座標は図に示すとおりである。

using SymPy

@syms r1::positive, r2::positive, x1::positive, x2::positive, y2::positive;

以下の 5 式を立て，解く

half = 1//2;
eq1 =x1^2 + r1^2 - (1 - r1)^2;
eq2 = x2^2 + y2^2 - (1 + r2)^2;
eq3 = x2^2 +(1 + r1 - y2)^2 - (r1 + r2)^2;
eq4 = x2^2 + (y2 - half - r1)^2 - (half + r1 - r2)^2;
eq5 = x1^2 + half^2 - (half + 2r1)^2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x1, x2, y2))

   1-element Vector{NTuple{5, Sym}}:
    (-1/2 + sqrt(2)/2, 1/6, sqrt(2 - sqrt(2)), sqrt(4/9 - 2*sqrt(2)/9), 1/6 + 2*sqrt(2)/3)

name = ["r1", "r2", "x1", "x2", "y2"]
for i in 1:5
   println("$(name[i]) = $(simplify(res[1][i])) = $(res[1][i].evalf())")
end

   r1 = -1/2 + sqrt(2)/2 = 0.207106781186548
   r2 = 1/6 = 0.166666666666667
   x1 = sqrt(2 - sqrt(2)) = 0.765366864730180
   x2 = sqrt(4 - 2*sqrt(2))/3 = 0.360797400097465
   y2 = 1/6 + 2*sqrt(2)/3 = 1.10947570824873

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)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x1, x2, y2) = (-1/2 + sqrt(2)/2, 1/6, sqrt(2 - sqrt(2)), sqrt(4/9 - 2*sqrt(2)/9), 1/6 + 2*sqrt(2)/3)
   println("r1 = $r1,  r2 = $r2,  x1 = $x1,  x2 = $x2,  y2 = $y2")
   plot()
   circle(0, 0, 1, beginangle=0, endangle=180)
   circle(0, 1 + r1, r1, :blue)
   circle( x1, r1, r1, :blue)
   circle(-x1, r1, r1, :blue)
   circle(0, 0.5, 0.5, :brown)
   circle(0, r1 + 0.5, r1 + 0.5, :magenta)
   circle( x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(x1, r1, "(x1,r1)", :blue, :center)
       point(0, 0.5, "1/2 ", :brown, :right)
       point(0, 0.5+r1, "r1+1/2 ", :magenta, :right)
       point(0, 1, "1 ", :black, :right)
       point(0, 1+r1, "1+r1 ", :brown, :right)
       point(0, 1+2r1, "1+2r1 ", :black, :center, :bottom)
       point(x2, y2, "(x2,y2)", :green, :center)
       plot!([x2, x2], [y2, y2+r2],color=:green, linewidth=0.25)
       annotate!(x2, y2+r2/2, text("r2 ", 10, :green, :right))
       vline!([0], color=:black, linewidth=0.25)
   end
end;