算額（その250）

算額（その250）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（234)
長野県北佐久郡軽井沢町峠 熊野神社 明治5年(1872)

外円の弦の上に 9 個の円が入っている。等円の径が 1 寸のとき，外円の径はいかほどか。

弦と y 軸の交点の座標を (0, y) とする。 
外円の半径，中心座標を r0, (0, 0) とする。
等円の半径，中心座標を r1, (x1, y + r1) とする。
右の下円の半径，中心座標を r2, (x2, y + r2) とする。
中央の下円の半径，中心座標を r2, (0, y + r2) とする。
上円の半径，中心座標を r3, (0, r0 - r3) とする。
上円と下円に挟まれる3個の円のうち，
   右側の円の半径，中心座標を r4, (x4, r4) とする。
   中央の円の半径，中心座標を r5, (0, r0 - 2r3 - r5) とする。

以下の関数 H(u) 中の連立方程式を解く。

r1 は変数のままで，nlsolve() を使うときに定義する。

include("julia-source.txt");

using Printf
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, x1, r2, x2, r3, r4, x4, y4, r5, y) = u
   return [
       x2^2 - (r0 - r2)^2 + (r2 + y)^2,  # eq1
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq2
       x1^2 - (r0 - r1)^2 + (r1 + y)^2,  # eq3
       -(r1 + r4)^2 + (x1 - x4)^2 + (r1 + y - y4)^2,  # eq4
       x1^2 - (r1 + r3)^2 + (-r0 + r1 + r3 + y)^2,  # eq5
       x4^2 - (r4 + r5)^2 + (r0 - 2*r3 - r5 - y4)^2,  # eq6
       x1^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq7
       x4^2 - (r2 + r4)^2 + (-r2 - y + y4)^2,  # eq8
       x4^2 - (r3 + r4)^2 + (r0 - r3 - y4)^2,  # eq9
       -r0 + 2*r2 + 2*r3 + 2*r5 + y,  # eq10
   ]
end;

r1 = 0.5
iniv = [big"2.1", 0.69, 0.24, 1.39, 0.29, 0.07, 0.12, 1.36, 0.05, 0.83]
res = nls(H, ini=iniv);

names = ["r0", "x1", "r2", "x2", "r3", "r4", "x4", "y4", "r5", "y"]
for i in 1:length(names)
   @printf("%2s = %.7f\n", names[i], res[1][i])
end
println("収束した? ",res[2])

   r0 = 2.0000000
   x1 = 0.6966214
   r2 = 0.2426407
   x2 = 1.3932428
   r3 = 0.2928932
   r4 = 0.0736862
   x4 = 0.1239249
   y4 = 1.3621096
   r5 = 0.0502525
    y = 0.8284271
   収束した? true

外円の半径は r0 = 2.0，元の単位で直径は 4 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, x1, r2, x2, r3, r4, x4, y4, r5, y) = res[1]
   @printf("r0 = %.5f\n", r0)
   plot()
   x = sqrt(r0^2 - y^2)
   segment(x, y, -x, y)
   circle(0, 0, r0, :black)
   circle(x1, y + r1, r1)
   circle(-x1, y + r1, r1)
   circle(x2, y + r2, r2, :blue)
   circle(-x2, y + r2, r2, :blue)
   circle(0, y + r2, r2, :blue)
   circle(0, r0 - r3, r3, :green)
   circle(0, y + 2r2 + r5, r5, :red)
   circle(x4, y4, r4, :black)
   circle(-x4, y4, r4, :black)
   if more == true
       point(x1, y + r1, "(x1,y+r1)", :red, :center, :top)
       point(x2, y + r2, "(x2,y+r2)", :blue, :center, :top)
       point(0, r0, "(0,r0)", :black, :left, :bottom)
       point(0, y, "(0,y)", :black, :left, :top)
       point(x, y, "(x,y)", :black, :left, :top)
       point(0, r0 - r3, "r0-r3", :green)
       point(0, y + r2, "y+r2", :green)
       plot!([0, -1.2r1], [y + 2r2 + r5, r0], line=(:arrow, :red, 0.5))
       point(-1.2r1, r0, "(0,y+2r+r5)", :red, :right, :bottom, mark=false)
       # segment(x4, y4, 1.5, r0, :black)
       plot!([x4, 1.2r1], [y4, r0], line=(:arrow, :black, 0.5))
       point(1.2r1, r0, "(x4,y4)", :black, :left, :bottom, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;