算額（その136）

算額（その136）

岐阜県郡上市 郡上八幡神社（小野天満宮）
https://toyo.repo.nii.ac.jp/index.php?action=pages_view_main&active_action=repository_action_common_download&item_id=2533&item_no=1&attribute_id=18&file_no=1&page_id=13&block_id=17

外円の中に，甲円，乙円，丙円，丁円が入っている。乙円，丁円の径がそれぞれ 9寸9分，2寸9分7厘のとき，外円の径はいかほどか。

乙円，丁円の半径をそれぞれ 990，297 とし，図のように記号を定め方程式を解く。

using SymPy

@syms r0, r1, x1::positive, y1::positive, r2::positive, r3, x3::positive, y3::negative, r3::positive, r4, x4, y4;
r3 = 990
r4 = 297
y1 = x1 / sqrt(Sym(3))
y3 = x3*sqrt(Sym(3))
y4 = x4*sqrt(Sym(3))
eq0 = r2 + 2r1 - r0 |> expand
eq1 = x1^2 + y1^2 - (r1 + r2)^2 |> expand
eq2 = x4^2 + y4^2 - (r2 + r4)^2 |> expand
eq3 = x3^2 + y3^2 - (r0 - r3)^2 |> expand
eq4 = x1^2 + y1^2 - (r0 - r1)^2 |> expand
eq5 = (x3 - x1)^2 + (y3 - y1)^2 - (r3 + r1)^2 |> expand
eq6 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2 |> expand;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r0, r1, x1, r2, x3, x4))

solve() では解が求まらないので，nlsolve() により数値解を求める。

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
println(eq6, ",")

   -r1^2 - 2*r1*r2 - r2^2 + 4*x1^2/3,
   -r2^2 - 594*r2 + 4*x4^2 - 88209,
   -r0^2 + 1980*r0 + 4*x3^2 - 980100,
   -r0^2 + 2*r0*r1 - r1^2 + 4*x1^2/3,
   -r1^2 - 1980*r1 + 4*x1^2/3 - 4*x1*x3 + 4*x3^2 - 980100,
   -r1^2 - 594*r1 + 4*x1^2/3 - 4*x1*x4 + 4*x4^2 - 88209,

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)
(r0, r1, x1, r2, x3, x4) = u
return [
       -r1^2 - 2*r1*r2 - r2^2 + 4*x1^2/3,
       -r2^2 - 594*r2 + 4*x4^2 - 88209,
       -r0^2 + 1980*r0 + 4*x3^2 - 980100,
       -r0^2 + 2*r0*r1 - r1^2 + 4*x1^2/3,
       -r1^2 - 1980*r1 + 4*x1^2/3 - 4*x1*x3 + 4*x3^2 - 980100,
       -r1^2 - 594*r1 + 4*x1^2/3 - 4*x1*x4 + 4*x4^2 - 88209,    
   ]
end;

iniv = [10000.0, 4000, 6000, 2000, 5000, 1500]
res = nls(H, ini=iniv)
println(res)

   ([10032.79204966505, 3808.477217382768, 5390.414765908734, 2415.8376148995144, 4521.396024832525, 1356.418807449757], true)

外円の半径 = 10032.79204966505
甲円の半径 = 3808.477217382768, x 座標 = 5390.414765908734
乙円の半径 = 2415.8376148995144
丙円の中心の x 座標 = 4521.396024832525
丁円の中心の x 座標 = 1356.418807449757

元の単位でいえば，外円の直径は 100寸3分3厘

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, x1, r2, x3, x4) = res[1]
   r3 = 990
   r4 = 297
   y1 = x1 / sqrt(3)
   y3 = x3 * sqrt(3)
   y4 = x4 * sqrt(3)
   plot()
   circle(0, 0, r0, :black)
   circle(x1, y1, r1, :orange)
   circle(x1, -y1, r1, :orange)
   circle(-x1, y1, r1, :orange)
   circle(-x1, -y1, r1, :orange)
   circle(0, r0 - r1, r1, :orange)
   circle(0, r1 - r0, r1, :orange)
   circle(0, 0, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(x3, -y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(-x3, -y3, r3, :green)
   circle(r0 - r3, 0, r3, :green)
   circle(r3 - r0, 0, r3, :green)
   circle(x4, y4, r4)
   circle(x4, -y4, r4)
   circle(-x4, y4, r4)
   circle(-x4, -y4, r4)
   circle(r2 + r4, 0, r4)
   circle(-r2 - r4, 0, r4)
   if more == true
       point(x1, y1, "甲:(x1,y1;r1)", :orange, :center)
       point(0, 0, "乙:(0,0;r2)", :blue, :center)
       point(x3, y3, " 丙:(x3,y3;r3)", :green)
       point(x4, y4, " 丁:(x4,y4;r4)", :red)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;