算額（その104）

算額（その104）

愛媛県四国中央市川之江町 大岡八幡神社 文化元年(1804)4月
http://www.wasan.jp/ehime/kawanoe2.html

外円の中に甲円，乙円が1個ずつ，丙円，丁円，戊円が2つずつ互いに接して入っている。
甲円，乙円の直径をそれぞれ66寸，44寸としたとき，戊円の直径はいくつか。

甲円，乙円の半径をそれぞれ 33，22 とする。丙円，丁円，戊円の中心座標と半径を (x1, y1, r1), (x2, y2, r2), (x3, y3, r3) として，以下の9本の方程式を解く。

条件式と未知数の数が等しいので，解けるはずだが solve() では解けない。

using SymPy

@syms x1::positive, y1::positive, r1::positive, x2::positive, y2::positive, r2::positive, x3::positive, y3::positive, r3::positive;
@syms x1, y1, r1, x2, y2, r2, x3, y3, r3
eq1 = x3^2 + (y3 - 33)^2 - (33 + r3)^2 |> expand
eq2 = x1^2 + (y1 - 33)^2 - (33 + r1)^2 |> expand
eq3 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq4 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2 |> expand
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq6 = x2^2 + (y2 - 88)^2 - (22 + r2)^2 |> expand
eq7 = x3^2 + (88 - y3)^2 - (22 + r3)^2 |> expand
eq8 = x2^2 + (y2 - 55)^2 - (55 - r2)^2 |> expand
eq9 = x1^2 + (y1 - 55)^2 - (55 - r1)^2 |> expand
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9]);

今までの経験で，nlsolve() で解くことにする。

eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println

   -r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3
   -r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1
   -r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2
   -r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2
   -r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2
   -r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260
   -r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260
   -r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2
   -r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1

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)
(x1, y1, r1, x2, y2, r2, x3, y3, r3) = u
return [
       -r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3,
       -r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1,
       -r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2,
       -r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2,
       -r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2,
       -r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260,
       -r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260,
       -r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2,
       -r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1
     ];
end;

iniv = [40.0, 57, 14,
       33,    82, 12,
       20,    68, 8];

res = nls(H, ini=iniv)

   ([40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996], false)

収束していないが，ftol=1e-14 が厳しすぎるのである。図に描いてみると十分な精度である。
戊円の直径は 7.5*2 = 15 である。元の単位で言えば 15 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 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")
   (x1, y1, r1, x2, y2, r2, x3, y3, r3) = [40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996]
   plot()
   circle(0, 55, 55)
   circle(0, 33, 33)
   circle(0, 88, 22)
   circle(x1, y1, r1, :green)
   circle(-x1, y1, r1, :green)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
    if more
       a = @sprintf("丙円 (%.1f, %.1f, %.1f)\n", x1, y1, r1)
       b = @sprintf("丁円 (%.1f, %.1f, %.1f)\n", x2, y2, r2)
       c = @sprintf("戊円 (%.1f, %.1f, %.1f)", x3, y3, r3)
       point(-20, 22, a * b * c, :black, mark=false)
       point(0, 33, " 甲 (0,33,33)", :red)
       point(0, 88, " 乙 (0,88,22)", :red)
       point(x1, y1, "    丙\n    (x1,y1;r1)", :green, :top)
       point(x2, y2, "    丁\n    (x2,y2;r2)", :blue, :top)
       point(x3, y3, "    戊\n    (x3,y3;r3)", :magenta, :top)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;