算額（その404）

算額（その404）

法道寺善の算変法「観術」
永井信一：早数教分科会レポート「和算と反転法について」
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/hanten.pdf

外円の中に 5 個の円が入っている。東円，西円，南円の直径がそれぞれ 3 寸，1 寸，2 寸のとき，外円，北円，心円の直径を求めよ。

外円の中心を原点に置き，西円の中心が x 軸上に来るように円を配置する。
東円の半径と中心座標を r1, (x1, y1)
西円の半径と中心座標を r2, (r0 - r2, 0)
南円の半径と中心座標を r3, (x3, y3)
北円の半径と中心座標を r4, (x4, y4)
心円の半径と中心座標を r5, (x5, y5)
外円の半径と中心座標を r0, (0, 0)
として連立方程式を立て nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::negative, y1::negative,
     r2::positive, r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::negative,
     r5::positive, x5::positive, y5::negative;

(r1, r2, r3) = (3, 1, 2) .// 2
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x3^2 + y3^2 - (r0 - r3)^2
eq3 = x4^2 + y4^2 - (r0 - r4)^2
eq4 = (x3 - x1)^2 + (y3 - y1)^2 - (r3 + r1)^2
eq5 = (x4 - x1)^2 + (y4 - y1)^2 - (r4 + r1)^2
eq6 = (x5 - x1)^2 + (y5 - y1)^2 - (r5 + r1)^2
eq7 = (r0 - r2 - x3)^2 + y3^2 - (r2 + r3)^2
eq8 = (r0 - r2 - x4)^2 + y4^2 - (r2 + r4)^2
eq9 = (r0 - r2 - x5)^2 + y5^2 - (r2 + r5)^2
eq10 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq11 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2;
# 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)
   (r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       x1^2 + y1^2 - (r0 - 3/2)^2,  # eq1
       x3^2 + y3^2 - (r0 - 1)^2,  # eq2
       x4^2 + y4^2 - (r0 - r4)^2,  # eq3
       (-x1 + x3)^2 + (-y1 + y3)^2 - 25/4,  # eq4
       -(r4 + 3/2)^2 + (-x1 + x4)^2 + (-y1 + y4)^2,  # eq5
       -(r5 + 3/2)^2 + (-x1 + x5)^2 + (-y1 + y5)^2,  # eq6
       y3^2 + (r0 - x3 - 1/2)^2 - 9/4,  # eq7
       y4^2 - (r4 + 1/2)^2 + (r0 - x4 - 1/2)^2,  # eq8
       y5^2 - (r5 + 1/2)^2 + (r0 - x5 - 1/2)^2,  # eq9
       -(r5 + 1)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq10
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq11
   ]
end;
iniv = [big"2.5", -0.9, -0.4, 0.9, 1.2, 0.6, 1.2, -1.3, 0.32, 1.0, -0.3]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[2.508041569829081421977001552141168300225929489745528932614456499113765038860902, -0.485951405752271286255384816114094927579869826307774264025836820627077860874849, -0.8831755418663212924359755945019721217093601316808355690943080396556492904135786, 1.010043911301963852566206096056080557623999906394429728983586117023084099186487, 1.119821715084321284656363526042374037999885992006369003023642241853344206655946, 0.6000000000000000000000000000000000000000012240791396860081807414352057058378827, 1.609242974712810880330524278490115654664770306212569723396079629646188139868433, -1.025163245797121287768208353426213271483677234175174482387089368453727261649156, 0.3262233880108996037951833703229864820441137963024028723844040481034973892347998, 1.20445443902395988951867252073114035193644842415907539462550103774599845159441, -0.1920750116506622924019613444332251808267586158530185578363832744817135812429724], true)

   r0 = 2.50804;  x1 = -0.485951;  y1 = -0.883176;  x3 = 1.01004;  y3 = 1.11982;  r4 = 0.6;  x4 = 1.60924;  y4 = -1.02516;  r5 = 0.326223;  x5 = 1.20445;  y5 = -0.192075
   北径 = 1.2;  外径 = 5.01608;  心径 = 0.652447

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (3, 1, 2) .// 2
   (r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   @printf("r0 = %g;  x1 = %g;  y1 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5)
   @printf("北径 = %g;  外径 = %g;  心径 = %g\n", 2r4, 2r0, 2r5)  
   plot()
   circle(0, 0, r0, :black)
   circle(x1, y1, r1)
   circle(r0 - r2, 0, r2, :blue)
   circle(x3, y3, r3, :orange)
   circle(x4, y4, r4, :magenta)
   circle(x5, y5, r5, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, " 東:r1,(x1,y1)", :red, :center, delta=-delta)
       point(r0 - r2, 0, " 西:r2\n(r0-r2,0)", :blue, :center, :bottom, delta=delta)
       point(x3, y3, " 南:r3,(x3,y3)", :orange, :center, delta=-delta)
       point(x4, y4, " 北:r4,(x4,y4)", :magenta, :center, delta=-delta)
       point(x5, y5, " 心:r5\n(x5,y5)", :brown, :center, :vcenter, delta=delta/2)
       point(0, r0, " r0", :black, :left, :top, delta=-delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;