算額（その436）

算額（その436）

橘田彌曾八元克 天明八年戊申2月(1788)

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の中に，甲円 3 個，乙円 2 個，丙円 1 個が入っている。甲円は弦に接している。
外円の直径が 3 寸 6 分 のとき，甲円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (2r1, a + r1)
乙円の半径と中心座標を r2, (r1, y2)
丙円の半径と中心座標を r3, (0, r0 - r3)
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, a::positive, r1::negative,
     r3::positive, r2::positive, y2::positive;
r0 = 36//20
eq1 = (2r1)^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = r1^2 + (y2 - a - r1)^2 - (r1 + r2)^2
eq3 = r1^2 + y2^2 - (r0 - r2)^2
eq4 = r1^2 + (r0 - r3 - y2)^2  - (r3 + r2)^2
eq5 = 2r3 + 2r1 + a - r0;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r3, r2, y2))

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)
   (a, r1, r3, r2, y2) = u
   return [
       4*r1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq1
       r1^2 - (r1 + r2)^2 + (-a - r1 + y2)^2,  # eq2
       r1^2 + y2^2 - (r0 - r2)^2,  # eq3
       r1^2 - (r2 + r3)^2 + (r0 - r3 - y2)^2,  # eq4
       a - r0 + 2*r1 + 2*r3,  # eq5
   ]
end;
r0 = 3.6/2
iniv = r0 .* [big"0.22", 0.24, 0.14, 0.12, 0.69]
res = nls(H, ini=iniv);
println([round(Float64(x), digits=6) for x in res[1]], " 収束:", res[2]);

   [0.330519, 0.500029, 0.234712, 0.282615, 1.43263] 収束:true

a = 0.330519;  r1 = 0.500029;  r3 = 0.234712;  r2 = 0.282615;  y2 = 1.43263

外円の直径が 3.6 のとき，甲円の直径 は 1.00006 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, r3, r2, y2) = res[1]
   @printf("a = %g;  r1 = %g;  r3 = %g;  r2 = %g;  y2 = %g\n", a, r1, r3, r2, y2)
   @printf("外円の直径 = %g;  甲円の直径 = %g\n", 2r0, 2r1)
   plot()
   circle(0, 0, r0, :gray)
   circle(0, a + r1, r1)
   circle(2r1, a + r1, r1)
   circle(-2r1, a + r1, r1)
   circle(0, r0 - r3, r3, :blue)
   circle(r1, y2, r2, :green)
   circle(-r1, y2, r2, :green)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, a + r1, " 甲円:r1,(0,a+r1)", :red, :center, :top, delta=-delta/2)
       point(2r1, a + r1, " (2r1,a+r1)", :red, :center, :bottom, delta=delta)
       point(0, r0 - r3, " 丙円:r3,(0,r0-r3)", :black, :center, :bottom, delta=delta)
       point(r1, y2, " 乙円:r2,(r1,y2)", :black, :left, :vcenter)
       point(0, a, " a", :magenta, :left, :top, delta=-delta/2)
       point(0, r0, " r0", :gray, :left, :bottom, delta=delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;