算額（その24）

算額（その24）

大きな円の中に，半径の異なる 4 種類の円が，図のように配置されている。各円の半径を求めよ。

岩手県東磐井郡藤沢町 藤勢寺 元治2年
http://www.wasan.jp/iwate/fujise.html

外側の円の半径を 2 とし，内部の円の中心座標と半径を小さい順に (0, y1, r1), (0, y2, r2), (x3, y3, r3), (1, 0, 1) と定める。

using SymPy

@syms r1::positive, r2::positive, r3::positive,
       y1::positive, y2::positive, x3::positive, y3::positive;

未知数は r1, r2, r3, x3, y1, y2, y3 の 7 個で，以下の 7 つの方程式を立てて解く。

eq1 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2;
eq2 = x3^2 + (y3 - y2)^2 - (r2 + r3)^2;
eq3 = 1 + y2^2 - (1 + r2)^2 |> expand
eq4 = (1 - x3)^2 + y3^2 - (1 + r3)^2;
eq5 = y1 - y2  - r1 - r2;
eq6 = y1 - 2 + r1;
eq7 = x3^2 + y3^2 - (2 - r3)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x3, y1, y2, y3));

解は以下のごとくである。

name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
   println("$(name[i]) = $(res[1][i]) = $(res[1][i].evalf())")
end

   r1 = 8/11 - 2*sqrt(5)/11 = 0.320714913181856
   r2 = 16*sqrt(5)/209 + 46/209 = 0.391277931291850
   r3 = 7/11 - sqrt(5)/11 = 0.433084729318201
   x3 = 1/11 + 3*sqrt(5)/11 = 0.700745812045397
   y1 = 2*sqrt(5)/11 + 14/11 = 1.67928508681814
   y2 = 68/209 + 60*sqrt(5)/209 = 0.967292242344437
   y3 = 2/11 + 6*sqrt(5)/11 = 1.40149162409079

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, x3, y1, y2, y3) = (8/11 - 2*sqrt(5)/11, 16*sqrt(5)/209 + 46/209, 7/11 - sqrt(5)/11, 1/11 + 3*sqrt(5)/11, 2*sqrt(5)/11 + 14/11, 68/209 + 60*sqrt(5)/209, 2/11 + 6*sqrt(5)/11)
   println("r1 = $r1,  r2 = $r2,  r3 = $r3")
   println("y1 = $y1,  y2 = $y2,  y3 = $y3")
   plot()
   circle(0, y1, r1, :green)
   circle(0, -y1, r1, :green)
   circle(0, y2, r2, :blue)
   circle(0, -y2, r2, :blue)
   circle( x3,  y3, r3, :magenta)
   circle(-x3,  y3, r3, :magenta)
   circle( x3, -y3, r3, :magenta)
   circle(-x3, -y3, r3, :magenta)
   circle( 1, 0, 1)
   circle(-1, 0, 1)
   if more
       point(0, y1, "(0,y1,r1)", :brown, :center)
       point(0, y2, "(0,y2,r2)", :blue, :center)
       point(x3, y3, "(x3,y3,r3)", :magenta, :center)
       point(1, 0, "(1,0,1)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.5)
   end
   circle(0, 0, 2, :black)
end;

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, x3, y1, y2, y3 = u
   return [
       x3^2 + (y1 - y3)^2 - (r1 + r3)^2,
       x3^2 + (y3 - y2)^2 - (r2 + r3)^2,
       1^2 + y2^2 - (1 +r2)^2,
       (1 - x3)^2 + y3^2 - (1 + r3)^2,
       (y1 - y2)  - (r1+ r2),
       y1 - 2 + r1,
       x3^2 + y3^2 - (2 - r3)^2
   ]
end;

iniv = [0.354644124017785, 0.33761490117551013, 0.4639474199979425, 0.6081577400061725, 1.3714562740209328, 0.8883769604436443, 1.4105324143047882]  # 初期値

res = nls(h, ini=iniv);

name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
   println("$(name[i]) = $(res[1][i])")
end

   r1 = 0.32071491318185646
   r2 = 0.3912779312918499
   r3 = 0.4330847293182009
   x3 = 0.7007458120453972
   y1 = 1.6792850868181435
   y2 = 0.9672922423444372
   y3 = 1.4014916240907944