算額（その23）

算額（その23）

直角三角形に接する大円，小円の径がそれぞれ 3，2 である。直角三角形の外形を描け。

山形県天童市山元 若松観音堂
http://www.wasan.jp/yamagata/wakamatu2.html

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

using SymPy

@syms x::positive, x1::positive, y1::positive, x2::positive, y2::positive, x3::positive;

小円の中心座標を(x2, 2) とする。大円と小円が接していることから，

eq1 = 1 + (x2 - 3)^2 - 25
x2 = solve(eq1)[1]
x2 |> println
x2.evalf() |> println # x2

   3 + 2*sqrt(6)
   7.89897948556636

△ABC, △DEC が相似なので，(x3 - x2) / (x3 - 3) = 2//3

x3 = solve((x3 - x2)/ (x3 - 3) - 2//3, x3)[1]
x3 |>  println # x3

   3 + 6*sqrt(6)

大円と (x1, y1) で接し，(x3, 0) を通る直線について，
円周上に存在するので，

eq2 = (x1 - 3)^2 + (y1 - 3)^2 - 9;
eq3 = (x3 - x1)^2 + y1^2 - (x3 - 3)^2;

x1, y1 = solve([eq2, eq3], (x1, y1))[1] # x1, y1

   (12*sqrt(6)/25 + 3, 144/25)

△FGC，△IOC の相似により，

I = y1*x3/(x3-x1) |> simplify |> together |> println # intercept

   12*(sqrt(6) + 12)/23

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")
   x3 = 3 + 6sqrt(6)
   x2 = 3 + 2sqrt(6)
   (x1, y1) = (12*sqrt(6)/25 + 3, 144/25)
   println("x2 = $x2")
   println("x1 = $x1;  x2 = $x2")
   intercept = 12(sqrt(6) + 12)/23
   println("x3 = $x3;  intercept = $intercept")
   plot([0, x3, 0, 0], [0, 0, intercept, 0], color=:black, linewidth=0.5)
   circle(3, 3, 3, :green)
   circle(x2, 2, 2, :magenta)
   if more
       point(3, 3, "A:(3,3)", :green, :center)
       point(3, 0, "\nB:(3,0)", :green, :bottom)
       point(x3, 0, " C:(x3,0)", :black, :left, :bottom)
       point(x2, 2, "D:(x2,2)", :magenta, :center)
       point(x2, 0, "\nE:(x2,0)", :magenta, :bottom)
       point(x1, y1, " F:(x1,y1)", :green, :left, :bottom)
       point(x1, 0, "\n\nG:(x1,y1)", :green, :left, :bottom)
       point(0, 0, " O", :black, :left, :bottom)
       point(0, intercept, " I(0,intercept)", :black, :left, :bottom)
   end
end;

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