算額（その740）
浅草観音堂の算額
会田安明：『他流諸国之表題集』
山口正義：やまぶき4，第84号
https://yamabukiwasan.sakura.ne.jp/ymbk84.pdf
直角三角形の中に，大円，中円，小円を入れる。
中円，小円の直径をそれぞれ 3.65767079 寸，1.95349271 寸としたとき，鈎，股，大円の直径を求めよ。
後述するが，1.92349271 の誤記である。
図形的には，算額（その453）と同じである。しかし，条件の付け方により，この問題は SymPy では解けず，NLsolve で数値解を求めるしかない。
直角三角形の直角を挟む二辺を「鈎」，「股」とする
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (r3, y3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1, r2, x2, r3, y3, 鈎, 股, 弦, d
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = r1*(股 - x2) - r2*(股 - r1)
eq4 = r1*(鈎 - y3) - r3*(鈎 - r1)
eq5 = 鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r1;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, 鈎, 股, x2, y3))
「問」で与えられた条件，中円直径 = 3.65767079，小円直径 = 1.95349271 では，「答」のように「きれいな数」は出ない。
r1, 鈎, 股 は 3.021927068553583, 8.099020725821061, 14.930767011037842 になる。
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 Float64.(v), r.f_converged
end;
function H(u)
(r1, 鈎, 股, x2, y3) = u
return [
(x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2,
(r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2,
r1*(股 - x2) - r2*(股 - r1),
r1*(鈎 - y3) - r3*(鈎 - r1),
鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r1
]
end;
r2 = 3.65767079/2
r3 = 1.95349271/2 # 誤記のまま。正しくは 1.92349271 である
iniv = BigFloat[6/2, 8, 15, 7.7, 6.4]
res = nls(H, ini=iniv)
([3.021927068553583, 8.099020725821061, 14.930767011037842, 7.723674481624338, 6.458004046320843], true)
「答」がきれいな数になるように，中円，小円の直径をあれほど細かく書いたのだろうに。
「答」の方から逆算すると，小円直径 = 1.92349271 の誤記であることがわかる。
訂正して，計算し直すと r1, 鈎, 股 は 3.000046427350913, 8.000208871784782, 14.999849358214375 となり，3, 8, 15 というきれいな数になる。
その他のパラメータは以下のとおりである。
r1 = 3.00005; 鈎 = 8.00021; 股 = 14.9998; x2 = 7.68474; y3=6.39727
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 3.65767079/2
r3 = 1.92349271/2
(r1, 鈎, 股, x2, y3) = res[1]
@printf("r1 = %g; 鈎 = %g; 股 = %g; x2 = %g; y3=%g\n", r1, 鈎, 股, x2, y3)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
circle(r1, r1, r1)
circle(x2, r2, r2, :blue)
circle(r3, y3, r3, :green)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta)
point(r3, y3, " 小円:r3,(r3,y3)", :green, :left, :vcenter)
point(股, 0, "股", :magenta, :left, :bottom, delta=delta)
point(0, 鈎, " 鈎", :magenta, :left, :bottom)
end
end;
