算額（その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;
