算額（その95）

算額（その95）

東京都台東区下谷池端 東淵寺 明治2年(1869)
http://www.wasan.jp/tokyo/tobuti.html

半円の中に等円 2 個と半円 1 個が入っている。外円の径を 2 寸 7 分 とすると等円の径はいくつか。

図のように記号を定め方程式を解く。

using SymPy

@syms R::positive, r::positive, x::positive, y::positive;

R = 27
eq1 = x^2 + (y - r)^2 - 4r^2 |> expand;
eq2 = x^2 + r^2 - (R - r)^2 |> expand;
eq3 = r^2 + y^2 - R^2 |> expand;
a = solve([eq1, eq2, eq3], (r, x, y));

方程式自体は簡単なものなのであるが，方程式の解 a をそのまま書き出すと，虚数単位 I が含まれたりする恐ろしい数式になる。それを .evalf() で数値になおすと以下の3組になる。
まだ I が残っているが 4.74210720465186e-23*I などで誤差の範囲で実数である。2 番目と 3 番目の解は不適切解である。
つまるところ r = 10.0929641600229, x = 13.5639203535985, y = 25.0426051852537 ということであろう。
よって等円の径は 1 寸 0 分 0 厘 9 毛である。

for i in 1:3
   println("r = $(a[i][1].evalf());  x = $(a[i][2].evalf());  y = $(a[i][3].evalf())")
end

   r = 10.0929641600229;  x = 13.5639203535985 + 4.74210720465186e-23*I;  y = 25.0426051852537 + 0.e-22*I
   r = -22.2346710141611 - 4.33680868994202e-19*I;  x = 43.928034724589 + 1.17139889977585e-22*I;  y = -15.3172910428713 - 0.e-21*I
   r = 17.5417068541382 + 4.33680868994202e-19*I;  x = 14.7733601500627*I;  y = -20.5253141423824 + 0.e-21*I

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; fontsize=10, mark=true)
  mark && scatter!([x], [y], color=color, markerstrokewidth=0)
  annotate!(x, y, text(string, fontsize, vertical, position, color))
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, r, x, y) = (27, 10.0929641600229, 13.5639203535985, 25.0426051852537)
   plot()
   circle(0, 0, R, :red, beginangle=0, endangle=180)
   circle(x, r, r, :blue)
   circle(-x, r, r, :blue)
   circle(0, y, r, :blue, beginangle=180, endangle=360)
   segment(-r, y, r, y, :blue)
   segment(-R, 0, R, 0, :red)
   if more
       point(0, y, "y ", :blue, :right)
       point(0, r, "r ", :blue, :right)
       point(x, r, "(x,r)", :blue, :top)
       point(x, 0, "  x", :blue, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;