算額（その22）

算額（その22）

バージョンアップ　https://blog.goo.ne.jp/r-de-r/e/b049f677e607ba6a92b04bb805ae3b0f

一辺の長さが 2 の正三角形に 3 種類の半径を持つ円が図のように接している。それぞれの円の半径を求めよ。

岩手県遠野市附馬牛町 早池峰神社 弘化3年6月(1846)
http://www.wasan.jp/iwate/hayatine.html

東京都 住吉神社 嘉永8年(1851) 小円 2 個のない簡易版
http://www.wasan.jp/tokyo/sumiyosi.html

ほかに必要とする座標は図に示すとおりである。

using SymPy

@syms r1::positive, r2::positive, r3::positive, x2::positive, y3::positive;

以下の 5 式を立て，solve() で解こうとするも，一向に計算が終了しない。

eq1 = r3 / (sqrt(Sym(3)) - y3) - sin(PI/6);
eq2 = r2/(1 - x2) - tan(PI/6);
eq3 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2;
eq4 = x2^2 + (y3 - r2)^2 - (r2 + r3)^2;
eq5 = (x2 - r1)^2 + (r2 - r1)^2 - (r1 + r2)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r3, x2, y3))

そこで，nlsolve() で数値解を求めることにする。

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
      r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params...), [ini], ftol=1e-14)
      v = r.zero[1]
   else
      r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
      v = r.zero
   end
   return v, r.f_converged
end;

function h(u)
   r1, r2, r3, x2, y3 = u
   return [
       r3 / (sqrt(3) - y3) - sin(pi/6),
       r2　/　(1 - x2) - tan(pi/6),
       r1^2 + (y3 - r1)^2 - (r1 + r3)^2,
       x2^2 + (y3 - r2)^2 - (r2 + r3)^2,
       (x2 - r1)^2 + (r2 - r1)^2 - (r1 + r2)^2
   ]
end;

iniv = [0.20, 0.3, 0.5, 0.5, 0.8]  # 初期値

nls(h, ini=iniv)

   ([0.15054780115874736, 0.2613764437119185, 0.4830337288182606, 0.5472827195892904, 0.7659833499323558], true)

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, r3, x2, y3) = [0.1505478011587473, 0.2613764437119185, 0.4830337288182606, 0.5472827195892904, 0.7659833499323557]
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2;  y3 = $y3")
   plot([-1, 1, 0, -1], [0, 0, sqrt(3), 0], color=:black, linewidth=0.25)
   circle(0, sqrt(3)-2r3, r3)
   circle(r1, r1, r1, :blue)
   circle(-r1, r1, r1, :blue)
   circle(1-sqrt(3)*r2, r2, r2, :green)
   circle(sqrt(3)*r2-1, r2, r2, :green)
   if more
       point(r1, r1, "(r1,r1)", :blue, :center)
       point(x2, r2, "(x2,r2)", :green, :center)
       point(0, y3, " y3", :red, :left)
       point(0, y3+r3, "\ny3+r3", :red, :left)
       point(0, sqrt(3), " √3", :red, :left)
       vline!([0], linewidth=0.5)
   end
end;