算額（その615）

算額（その615）

福島県楢葉町北田 北田神社 明治27年(1894)
ならはの絵馬―村人の祈り― 018/036page
http://is2.sss.fukushima-u.ac.jp/fks-db/txt/10088.002/html/00018.html

福島県楢葉町北田には，北田天満宮と大山祇神社があるが北田神社はないようだ。

二円が交わってできる三日月型の中に仁，禮，義，智の 4 円が入っている。仁円と禮円の直径の和と智円の直径の積は 3 寸，智円と禮円の直径の和と仁円の直径の積は 4 寸，禮円の面積は 2 歩である。義円の面積はいかほどか。

三日月形を形成する2円の半径を r0，中心間の距離を x とする。
仁円の半径，中心座標を r1, (x1, y1)
禮円の半径，中心座標を r2, (x2, y2)
義円の半径，中心座標を r3, (x3, y3)
智円の半径，中心座標を r4, (x4, y4)
仁円と禮円の直径の和と智円の直径の積を「仁義和乗智」
智円と禮円の直径の和と仁円の直径の積を「智禮和乗仁」
禮円の面積を「禮円積」
とおき，以下の連立方程1式の数値解を nlsolve() で求める。

なお，算額では「仁義和乗智」= 3，「智禮和乗仁」= 4 の場合の解を要求しているが，算額の図は「仁義和乗智」= 4，「智禮和乗仁」= 3 の場合のもののようである。上に示した図も後者の設定での解である。

実際の解は以下の通り。

include("julia-source.txt");

using SymPy

@syms r0, x, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4,
     仁義和乗智, 智禮和乗仁, 禮円積
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = x3^2 + y3^2 - (r0 - r3)^2
eq4 = x4^2 + y4^2 - (r0 - r4)^2
eq5 = (x1 - x)^2 + y1^2 - (r0 + r1)^2
eq6 = (x2 - x)^2 + y2^2 - (r0 + r2)^2
eq7 = (x3 - x)^2 + y3^2 - (r0 + r3)^2
eq8 = (x4 - x)^2 + y4^2 - (r0 + r4)^2
eq9 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq10 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq11 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq12 = 2(r1 + r3)*2r4 - 仁義和乗智
eq13 = 2(r4 + r2)*2r1 - 智禮和乗仁
eq14 = (2r2)^2 * 0.79 - 禮円積;

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, x, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4) = u
   return [
       x1^2 + y1^2 - (r0 - r1)^2,  # eq1
       x2^2 + y2^2 - (r0 - r2)^2,  # eq2
       x3^2 + y3^2 - (r0 - r3)^2,  # eq3
       x4^2 + y4^2 - (r0 - r4)^2,  # eq4
       y1^2 - (r0 + r1)^2 + (-x + x1)^2,  # eq5
       y2^2 - (r0 + r2)^2 + (-x + x2)^2,  # eq6
       y3^2 - (r0 + r3)^2 + (-x + x3)^2,  # eq7
       y4^2 - (r0 + r4)^2 + (-x + x4)^2,  # eq8
       -(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2,  # eq9
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq10
       -(r3 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq11
       2*r4*(2*r1 + 2*r3) - 仁義和乗智,  # eq12
       2*r1*(2*r2 + 2*r4) - 智禮和乗仁,  # eq13
       3.16*r2^2 - 禮円積,  # eq14
   ]
end;
(仁義和乗智, 智禮和乗仁, 禮円積) = (3, 4, 2)
iniv = BigFloat[3.68, -1.6, 0.76, 2.7, 1.11, 0.8, 2.85, -0.43, 0.69, 2.36, -1.84, 0.52, 1.57, -2.75]
res = nls(H, ini=iniv)
println([round(Float64(x), digits=10) for x in res[1]], " 収束? ", res[2])

   [3.6827310728, -1.6028875866, 0.7619151516, 2.6996483947, 1.1149278859, 0.7955572842, 2.8542380589, -0.4348535258, 0.6889728183, 2.3644695589, -1.8362658282, 0.5169248182, 1.5738882089, -2.7468536816] 収束? true

義円の面積は 3.16*r3^2 ≒ 1.5 である。

その他のパラメータは以下の通り。
r0 = 3.68273;  x = -1.60289;
r1 = 0.761915;  x1 = 2.69965;  y1 = 1.11493
r2 = 0.795557;  x2 = 2.85424;  y2 = -0.434854
r3 = 0.688973;  x3 = 2.36447;  y3 = -1.83627
r4 = 0.516925;  x4 = 1.57389;  y4 = -2.74685

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, x, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4) = res[1]
   @printf("義円の面積 = %g;  r0 = %g;  x = %g;
r1 = %g;  x1 = %g;  y1 = %g
r2 = %g;  x2 = %g;  y2 = %g
r3 = %g;  x3 = %g;  y3 = %g
r4 = %g;  x4 = %g;  y4 = %g\n",
       3.16*r3^2, r0, x, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4)
   plot()
   circle(0, 0, r0, :black)
   circle(-x, 0, r0, :black)
   circle(-x1, y1, r1, :black)
   circle(-x2, y2, r2, :red)
   circle(-x3, y3, r3, :orange)
   circle(-x4, y4, r4, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(-x1, y1, "仁:r1,(x1,y1)", :black, :center, :vcenter, mark=false)
       point(-x2, y2, "禮:r2,(x2,y2)", :red, :center, :vcenter, mark=false)
       point(-x3, y3, "義:r3,(x3,y3)", :orange, :center, :vcenter, mark=false)
       point(-x4, y4, "智:r4,(x4,y4)", :green, :center, :vcenter, mark=false)
       plot!(xlims=(-4, 1))
   end
end;