算額（その354）

算額（その354）

福島県楢葉町北田 北田神社 明治27年(1894)
ならはの絵馬―村人の祈り― 018/036page
http://is2.sss.fukushima-u.ac.jp/fks-db/txt/10088.002/html/00018.html
福島県楢葉町北田には，北田天満宮と大山祇神社があるが北田神社はないようだ。

外円内に春夏秋冬の6円を入れる。春円，夏円の直径をそれぞれ 3 寸，4 寸としたとき，秋円の直径を求めよ。

外円の半径，中心座標を r0, (0, 0)
春円の半径，中心座標を r1, (0, r0 - r1)
夏円の半径，中心座標を r2, (x2, y2)
秋円の半径，中心座標を r3, (r3, y3)
冬円の半径，中心座標を r4, (0, r0 - 2r1 - r4)
とおき，以下の連立方程式の数値解を nlsolve() で求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive,
     y2::positive, r3::positive, y3::negative, r4::positive;

(r1, r2) = (3, 4) .// 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (r0 - 2r1 - r4 - y2)^2 - (r2 + r4)^2
eq3 = x2^2 + y2^2 - (r0 - r2)^2
eq4 = (x2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = r3^2 + y3^2 - (r0 - r3)^2
eq6 = r3^2 + (r0 - 2r1 - r4 - y3)^2 - (r3 + r4)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6]);

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, x2, y2, r3, y3, r4) = u
   return [
       x2^2 + (r0 - y2 - 3/2)^2 - 49/4,  # eq1
       x2^2 - (r4 + 2)^2 + (r0 - r4 - y2 - 3)^2,  # eq2
       x2^2 + y2^2 - (r0 - 2)^2,  # eq3
       (-r3 + x2)^2 - (r3 + 2)^2 + (y2 - y3)^2,  # eq4
       r3^2 + y3^2 - (r0 - r3)^2,  # eq5
       r3^2 - (r3 + r4)^2 + (r0 - r4 - y3 - 3)^2,  # eq6
   ]
end;

iniv = [28.0, 7.2, 8.0, 6.0, -5.2, 2.8]
res = nls(H, ini=iniv);

   r0 = 8.69972; x2 = 3.23607; y2 = 5.86635; r3 = 4.34164; y3 = -0.378172; r4 = 1.77267

秋円の直径は 8.68328 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (3, 4) .// 2
   (r0, x2, y2, r3, y3, r4) = res[1]
   @printf("r0 = %g; x2 = %g; y2 = %g; r3 = %g; y3 = %g; r4 = %g\n", r0, x2, y2, r3, y3, r4)
   @printf("秋円直径 = %g 寸\n", 2r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :green)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   circle(r3, y3, r3, :gray)
   circle(-r3, y3, r3, :gray)
   circle(0, r0 - 2r1 - r4, r4, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r0 - r1, "春:r1\nr0-r1", :green, :left, :vcenter)
       point(x2, y2, "夏:r2,(x2,y2)", :red, :center, delta=-delta)
       point(0, r0 - 2r1 - r4, "冬:r4\nr0-2r1-r4", :orange, :left, :vcenter)
       point(r3, y3, "秋:r3,(r3,y3)", :gray, :center, delta=-delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;