算額（その38）

算額（その38）

新潟県糸魚川市 天津神社 寛政12年
http://www.wasan.jp/niigata/tensin.html

図のように，外円の中に 5 種類の径（木，火，土，金，水）を持つ円が収まっている。土円，金円の径はそれぞれ 9分9釐（厘），一寸9分9釐（厘）である。水円の径を求めよ。

図のように記号を定め，方程式を解く。水円の径は r3 である。（プログラムでは，すべて直径ではなく半径を使っている。）

using SymPy
@syms　r1::positive, r2::positive, y2::positive, r0::positive, y3::positive, r3::positive, r4::positive, x4::positive, y4::positive, r5::positive, x5::positive, y5::positive;

r1 = 99
r2 = 198
y3 = r3
x5 = r5
x4 = r4
eq1 = x4^2+(2r0-r1-y4)^2-(r1+r4)^2
eq2 = x4^2+(y4-y2)^2-(r2+r4)^2
eq3 = x5^2+(y2-y5)^2-(r2+r5)^2
eq4 = x5^2+(y5-y3)^2-(r3+r5)^2
eq5 = x4^2+(y4-r0)^2-(r0-r4)^2b
eq6 = x5^2+(r0-y5)^2-(r0-r5)^2
eq7 = (x5 - x4)^2 + (y4 - y5)^2 - (r4 + r5)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (y2, r0, r3, r4, y4, r5, y5))

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)
   y2, r0, r3, r4, y4, r5, y5 = u
   return [
       r4^2 - (r4 + 99)^2 + (2*r0 - y4 - 99)^2
       r4^2 - (r4 + 198)^2 + (-y2 + y4)^2
       r5^2 - (r5 + 198)^2 + (y2 - y5)^2
       r5^2 + (-r3 + y5)^2 - (r3 + r5)^2
       r4^2 + (-r0 + y4)^2 - (r0 - r4)^2
       r5^2 - (r0 - r5)^2 + (r0 - y5)^2
       (-r4 + r5)^2 - (r4 + r5)^2 + (y4 - y5)^2
 ];
end;

iniv = [2260., 1521, 580, 326, 2670, 752, 1679];

res = nls(h, ini=iniv)

y2, r0, r3, r4, y4, r5, y5 = res[1]

   7-element Vector{Float64}:
    2260.2616848711004
    1521.243310737389
     580.0010348927411
     326.36779179808775
    2670.68298737871
     752.3774162147455
    1679.6194053038619

求める水円の径は 580.0010348927411 である。

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 = 99
   r2 = 198
   y2, r0, r3, r4, y4, r5, y5 = [2260.2616848711004, 1521.243310737389, 580.0010348927411, 326.36779179808775, 2670.68298737871, 752.3774162147455, 1679.6194053038619]
   y3 = r3
   x4 = r4
   x5 = r5
   println("r0 = $r0,  r1 = $r1,  r2 = $r2,  r3 = $r3,  r4 = $r4,  r5 = $r5")
   println("y2 = $y2,  y3 = $y3,  x4 = $x4,  y4 = $y4,  x5 = $x5,  y5 = $y5")
   plot()
   circle(0, r0, r0, :green)
   circle(0, 2r0-r1, r1)
   circle(0, y2, r2, :blue)
   circle(0, y3, r3, :brown)
   circle(x4, y4, r4, :magenta)
   circle(-x4, y4, r4, :magenta)
   circle(x5, y5, r5, :black)
   circle(-x5, y5, r5, :black)
   if more
       point(0, r0, "r0 ", :green, :right)
       point(0, 2r0-r1, "r1 ", :red, :right)
       point(0, y2, "y2 ", :blue, :right)
       point(0, y3, "r3 ", :brown, :right)
       point(x4, y4, "(x4,y4)", :magenta, :center)
       point(x5, y5, "(x5,y5)", :black, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;