算額（その26 ）

算額（その26）

福島県田村郡三春町馬場 三春大神宮 明治 9 年
http://www.wasan.jp/fukusima/miharudaijingu.html

図のように，外円の中に 5 種類の半径を持つ円（木，火，土，金，水）が配置されている。「土」の径を 1 寸とすると，「金」の径は如何ほどか。

以下のように記号を定め，方程式を立てる。

外円の半径を 1 とする。「水」の中心の x 座標を x4，「金」の中心座標を (x5, r5) と置く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, r4::positive, r5::positive, x4::positive, x5::positive;

eq1 = x5^2 + (1 - r1 - r5)^2 - (r1 + r5)^2; # 木-金
eq2 = x5^2 + (r5 -(1 - 2r1 - r2))^2 - (r2 + r5)^2; # 火-金
eq3 = (x5 - x4)^2 + r5^2 - (r4 + r5)^2; # 水-金
eq4 = x4^2 + (r2 + 2r3)^2 - (r2 + r4)^2; # 火-水
eq5 = x4^2 + r3^2 - (r3 + r4)^2; # 土-水 
eq6 = 2(r1 + r2 + r3) - 1; # 半径の合計に関する制約
eq7 = x5^2 + r5^2 - (1 - r5)^2; # 金が外円に内接すること

しかし，なぜか solve() では解けない。

そこで，これも 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, r4, r5, x4, x5 = u
  return [
       x5^2 + (1 - r1 - r5)^2 - (r1 + r5)^2, # 木-金
       x5^2 + (r5 -(1 - 2r1 - r2))^2 - (r2 + r5)^2, # 火-金
       (x5 - x4)^2 + r5^2 - (r4 + r5)^2, # 水-金
       x4^2 + (r2 + 2r3)^2 - (r2 + r4)^2, # 火-水
       x4^2 + r3^2 - (r3 + r4)^2, # 土-水 
       2(r1 + r2 + r3) - 1, # 半径の合計に関する制約
       x5^2 + r5^2 - (1 - r5)^2 # 金が外円に内接すること
  ];
end;

iniv = [0.25, 0.2, 0.1, 0.2, 0.4, 0.2, 0.5]
res = nls(h, ini=iniv)[1];
r1, r2, r3, r4, r5, x4, x5 = res;

name = ["r1", "r2", "r3", "x3", "r4", "r5", "x4", "x5"]
for i = 1:7
  println("$(name[i]) = $(res[i])")
end

   r1 = 0.2795863507293534
   r2 = 0.1773223386410461
   r3 = 0.04309131062960051
   x3 = 0.14151590982250398
   r4 = 0.3602068246353233
   r5 = 0.1795075619334039
   x4 = 0.5287592559278309

r5/r3

   8.359152213576166

答え：「金」の径は，約 8.36 寸である。

算額には「9 寸」と描かれている。

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 draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2, r3, r4, r5, x4, x5 = [0.27958635072935356, 0.177322338641047, 0.04309131062959942, 0.14151590982250395, 0.3602068246353232, 0.179507561933404, 0.528759255927831]   
   println("r1 = $r1;  r2 = $r2;  r3 = $r3")
   println("r4 = $r4;  r5 = $r5")
   println("x4 = $x4;  x5 = $x5")
   plot()
   circle(0, 1-r1, r1, :green)
   circle(0, r1-1, r1, :green)
   circle(0, r2+2r3, r2, :magenta)
   circle(0, 2r1+r2-1, r2, :magenta)
   circle(0,  r3, r3, :red)
   circle(0, -r3, r3, :red)
   circle( x4, 0, r4, :black)
   circle(-x4, 0, r4, :black)
   circle( x5,  r5, r5, :brown)
   circle(-x5,  r5, r5, :brown)
   circle( x5, -r5, r5, :brown)
   circle(-x5, -r5, r5, :brown)
   if more
       point(0.15, 0.9, "木", :green, mark=false)
       point(0.05, 0.4, "火", :magenta, mark=false)
       point(0, -0.01, "土", :red, :center, mark=false)
       point(0.6, 0.6, "金", :brown, mark=false)
       point(0.2, 0.1, "水", :black, mark=false)
       point(0, 1-r1, "1-r1 ", :green, :right)
       point(0, r2+2r3, "r2+2r3 ", :magenta, :right)
       point(0, r3, "r3  ", :red, :right)
       point(x4, 0, "x4", :black)
       point(x5, r5, " (x5,r5)", :brown)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
   circle(0, 0, 1, :black,)
end;

算額（その25）

算額（その25）

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

図のように，外円にひし形の相対する 2 つの頂点が内接し，そのひし形の内部に大小の円が互いに外接している。大円の径が 4 寸のとき，小円の径はいかほどか。

以下のように記号を定め，方程式を立てる。

大円の半径を 1 とする。
x0 は，大円と小円の共通する接線（接点をそれぞれ (x1,y1), (x2, y2) とする）が x 軸と交わるときの x 座標である。

using SymPy

@syms x1::positive, y1::positive, x2::positive, y2::positive, r2::positive, x0::positive;

eq1 = (2 - y1) / (-x1) * y1 / (x1 - 1) + 1; # 接線と半径が直交する
eq2 = (x0-2 - r2) / (x0 - 1) - y2/y1; # 相似関係
eq3 = (x0-2 - r2) / (x0 - 1) - r2; # 相似関係
eq4 = (x1 - 1)^2 + y1^2 - 1; # (x1, y1) が大円の円周上にある
eq5 = y1/(x0-x1) - y2/(x0 - x2); # 相似関係
eq6 = y1 / (x0 - x1) - 2/x0; # 相似関係

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

   1-element Vector{NTuple{6, Sym}}:
    (1/4, 8/5, 4/5, 12/5, 1/5, 8/3)

name = ["r2", "x1", "y1", "x2", "y2", "x0"]
for (name0, res0) in zip(name, res[1])
   println("$name0 = $res0 = $(res0.evalf())")
end

   r2 = 1/4 = 0.250000000000000
   x1 = 8/5 = 1.60000000000000
   y1 = 4/5 = 0.800000000000000
   x2 = 12/5 = 2.40000000000000
   y2 = 1/5 = 0.200000000000000
   x0 = 8/3 = 2.66666666666667

答え：小円の半径は大円の半径の 1/4 である。

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")
   r2, x1, y1, x2, y2, x0 = (1/4, 8/5, 4/5, 12/5, 1/5, 8/3)    
   println("r2 = $r2")
   println("x1 = $x1;  y1 = $y1;  x2 = $x2;  y2 = $y2")
   println("x0 = $x0")
   plot()
   circle(1, 0, 1, :green)
   circle(-1, 0, 1, :green)
   circle(2 + r2, 0, r2, :blue)
   circle(-2 - r2, 0, r2, :blue)
   plot!([0, x0, 0, -x0, 0], [2, 0, -2, 0, 2], color=:black, linewidth=0.5)
   if more
       point(0, 2, "2 ", :black, :right, :bottom)
       point(x1, y1, "(x1,y1) ", :green, :right, :top)
       point(x2, y2, "(x2,y2)", :blue, :left, :bottom)
       point(2 + r2, 0, "2+r2", :blue, :center, :top)
       point(x0, 0, "  x0", :magenta, :center)
       point(1, 0, "1   ", :green, :top)
       point(2, 0, "2   ", :green, :top)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
   circle(0, 0, 2, :red)
end;