算額（その405）

算額（その405）

長野県北佐久郡軽井沢町峠 熊野神社 明治5年(1872)
山口正義：「やまぶき3」
https://yamabukiwasan.sakura.ne.jp/ymbk42.pdf

外円の中に，元，利，貞，亨，無名の円が入っている。外円の直径が 16 寸のとき，亨円の直径を求めよ。

外円の中心を原点に置く。
外円の半径と中心座標を r0, (0, 0)
元円の半径と中心座標を r1, (0, r0 - r1) と (0, r0 - 3r1)
利円の半径と中心座標を r2, (x2, y2)
貞円の半径と中心座標を r3, (x3, y3)
亨円の半径と中心座標を r4, (0, r4 - r0)
無名円の半径と中心座標を r5, (x5, y5)
として連立方程式を立て nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

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

r0 = 16//2
eq1 = x2^2 + y2^2 - (r0 - r2)^2
eq2 = x3^2 + y3^2 - (r0 - r3)^2
eq3 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq4 = x5^2 + (r0 - 3r1 - y5)^2 - (r1 + r5)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq7 = x3^2 + (r4 - r0 - y3)^2 - (r3 + r4)^2
eq8 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq9 = x5^2 + (r4 - r0 - y5)^2 - (r4 + r5)^2
eq10 = 4r1 + 2r4 - 2r0
eq11 = r0 - 2r1 - y2
# res = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11))

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

function H(u)
   (r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5) = u
   return [
x2^2 + y2^2 - (8 - r2)^2,  # eq1
x3^2 + y3^2 - (8 - r3)^2,  # eq2
x2^2 - (r1 + r2)^2 + (-r1 - y2 + 8)^2,  # eq3
x5^2 - (r1 + r5)^2 + (-3*r1 - y5 + 8)^2,  # eq4
-(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq5
-(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2,  # eq6
x3^2 - (r3 + r4)^2 + (r4 - y3 - 8)^2,  # eq7
-(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq8
x5^2 - (r4 + r5)^2 + (r4 - y5 - 8)^2,  # eq9
4*r1 + 2*r4 - 16,  # eq10
-2*r1 - y2 + 8,  # eq11
      ]
end;
r0 = 16//2
iniv = [big"1.9", 2.49, 4.3, 3.37, 2.29, 5.61, -1.28, 3.9, 1.1, 2.35, 0.34]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.750000000000000000000000000000000000000000001683866081422877969875418056451415, 2.243589743589743589743589743589743589743589745370092093352019975896094213312947, 3.589743589743589743589743589743589743589743594636276696631029402293443881749929, 4.499999999999999999999999999999999999999999996632267837154244060249163887074371, 2.348993288590604026845637583892617449664429532022065401438839115069164780141786, 5.637583892617449664429530201342281879194630878160421974850153126915301979093254, 0.3892617449664429530201342281879194630872483133777124317852969056382365325617382, 4.499999999999999999999999999999999999999999996632267837154244060249163887074371, 1.017441860465116279069767441860465116279069768415921597062855408025815041846539, 2.44186046511627906976744186046511627906976744452530995904551667578153181660868, 1.447674418604651162790697674418604651162790691173766512477690728560896338607614], true)

r1 = 1.75;  r2 = 2.24359;  x2 = 3.58974;  y2 = 4.5;  r3 = 2.34899;  x3 = 5.63758;  y3 = 0.389262;  r4 = 4.5;  r5 = 1.01744;  x5 = 2.44186;  y5 =1.44767

亨円の直径は 9 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 16//2
   (r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5) = res[1]
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  r5 = %g;  x5 = %g;  y5 =%g\n", r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5)
   @printf("亨円の直径 = %g\n", 2r4)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r0 - 3r1, r1, :blue)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
   circle(0, r4 - r0, r4, :orange)
   circle(x5, y5, r5, :brown)
   circle(-x5, y5, r5, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " r0-r1", :blue)
       point(0.1r1, r0 - r1, " 元", :blue, :left, :bottom, delta=delta, mark=false)
       point(0, r0 - 3r1, " r0-3r1", :blue)
       point(x2, y2, " 利:r2,(x2,y2)", :green, :center, delta=-delta)
       point(x3, y3, " 貞:r3,(x3,y3)", :magenta, :center, :bottom, delta=delta/2)
       point(0, r4 - r0, " r4-r0  亨", :orange)
       point(x5, y5, " 無:r5,(x5,y5)", :brown, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;