算額（その398）

算額（その398）

東京都府中市 大國魂神社 明治18年(1885)
http://www.f.waseda.jp/takezawa/sousuukyo/paper/2011ronbun.pdf

半円内に甲円，乙円，丙円，丁円，戊円の 9 円をいれる。外円の直径が 10 寸のとき，丙円の直径を求めよ。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
以下の連立方程式を立て，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive, r5::positive, x5::positive, y5::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = x5^2 + (y5 - r1)^2 - (r1 + r5)^2
eq4 = (x4 - x2)^2 + (y4 - r2)^2 - (r2 + r4)^2
eq5 = (x2 - x5)^2 + (y5 - r2)^2 - (r2 + r5)^2
eq6 = (x4 - x3)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = x2^2 + r2^2 - (2r1 - r2)^2
eq10 = x3^2 + y3^2 - (2r1 - r3)^2
eq11 = x4^2 + y4^2 - (2r1 - r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r1, r3, x3, r4, x41, r42, x42, y42))

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=big"1e-14")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-14")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r2, x2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       x3^2 + (-r1 + y3)^2 - (r1 + r3)^2,  # eq2
       x5^2 + (-r1 + y5)^2 - (r1 + r5)^2,  # eq3
       (-r2 + y4)^2 - (r2 + r4)^2 + (-x2 + x4)^2,  # eq4
       (-r2 + y5)^2 - (r2 + r5)^2 + (x2 - x5)^2,  # eq5
       -(r3 + r4)^2 + (-x3 + x4)^2 + (y3 - y4)^2,  # eq6
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq7
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq8
       r2^2 + x2^2 - (2*r1 - r2)^2,  # eq9
       x3^2 + y3^2 - (2*r1 - r3)^2,  # eq10
       x4^2 + y4^2 - (2*r1 - r4)^2,  # eq11
   ]
end;

r1 = 5//2
iniv = [big"1.3", 3.45, 0.52, 2.76, 3.36, 0.34, 3.62, 2.93, 0.26, 2.84, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.25, 3.535533905932737622004221810526050308609409471587999437373189411634045676321378, 0.5000000000000000000914962217596418053047474197754641655637114736626653004349011, 2.828427124746190097802451876811358797663566934419370521520740609153604311115202, 3.49999999999999999972551133472107458408575774067360750330886560977611543609143, 0.4166666666666666664757008581209034940070734535154511572023238905891655189934959, 3.535533905932737622365050944667911778662710884062463557643917323522673436115187, 2.916666666666666666600917687737661690052746729593517424887121308994452716368222, 0.3333333333333333328134030564312976635069374013893610682603843011810001642030419, 2.828427124746190097243674263049636031780028552049905775368715116058220501682123, 2.666666666666666666169204262275820633241336297105830855213052541023599670593682], true)

r2 = 1.25;  x2 = 3.53553;  r3 = 0.5;  x3 = 2.82843;  y3 = 3.5;  r4 = 0.416667;  x4 = 3.53553;  y4 = 2.91667;  r5 = 0.333333;  x5 = 2.82843;  y5 = 2.66667
丙円の直径 = 1

外円の直径が 10 寸のとき，丙円の直径は 1 寸である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   @printf("r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", r2, x2, r3, x3, y3, r4, x4, y4, r5, x5, y5)
   @printf("丙円の直径 = %g\n", 2r3)
   plot()
   circle(0, 0, 2r1, beginangle=0, endangle=180, :black)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(x4, y4, r4, :orange)
   circle(x5, y5, r5, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1, "甲円  r1: ", :red, :right, :vcenter)
       point(x2, r2, " 乙円:r2,(x2,r2)", :blue, :center, delta=-delta)
       point(x3, y3, "丙円:r3,(x3,y3) ", :green, :right, :vcenter)
       point(x4, y4, "丁円:r4,(x4,y4) ", :orange, :right, :vcenter)
       point(x5, y5, "戊円:r5,(x5,y5) ", :magenta, :right, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;