算額（その2019）
(14) 京都府城陽市寺田 水度神社 明治18年(1885)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円1個，正方形5個
直角三角形の中に，木火土金水の 5 円を容れる。金円，水円の直径がそれぞれ 2.5 寸，1 寸のとき，土円の直径はいかほどか。
この問題において，直角三角形の中に入っているという条件は無用の飾りである。「問」，「答」，「術」のいずれも，直角三角形についてはなんの言及もない。実際，条件を満たすような円のパラメータが得られた後で，それらが内接する直角三角形は容易に求めることができる。問題の本質は，算額(その31)，算額（その443）と同じである。
ここでは，5 円を内包する直角三角形も同定しよう（SymPy の能力的に，数値解を求めるしかできないが）。
式の簡略化のために，図を左右反転させて考える。
直角三角形の直角を挟む二辺の短い方を 「鈎」，長い方を 「股」
木円の半径と中心座標を r1, (r1, r1)
火円の半径と中心座標を r2, (x2, r2)
土円の半径と中心座標を r3, (x3, r3)
金円の半径と中心座標を r4, (x4, r4)
水円の半径と中心座標を r5, (x5, r5)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, r2::positive, x2::positive,
r3::positive, x3::positive,
r4::positive, x4::positive,
r5::positive, x5::positive,
鈎::positive, 股::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq3 = (x4 - r1)^2 + (r1 - r4)^2 - (r1 + r4)^2 |> expand
eq4 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq5 = (x2 - x5)^2 + (r2 - r5)^2 - (r2 + r5)^2 |> expand
eq6 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2 |> expand
eq8 = dist2(0, 鈎, 股, 0, r1, r1, r1)
eq9 = dist2(0, 鈎, 股, 0, x2, r2, r2);
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (r1, r2, x2, r3, x3, x4, x5, 鈎, 股))
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, r2, x2, r3, x3, x4, x5, 鈎, 股) = u
return [
(-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2, # eq1
(-r1 + x3)^2 + (r1 - r3)^2 - (r1 + r3)^2, # eq2
(-r1 + x4)^2 + (r1 - r4)^2 - (r1 + r4)^2, # eq3
(r2 - r3)^2 - (r2 + r3)^2 + (x2 - x3)^2, # eq4
(r2 - r5)^2 - (r2 + r5)^2 + (x2 - x5)^2, # eq5
(r3 - r4)^2 - (r3 + r4)^2 + (x3 - x4)^2, # eq6
(r3 - r5)^2 - (r3 + r5)^2 + (-x3 + x5)^2, # eq7
股*鈎*(2*r1^2 - 2*r1*股 - 2*r1*鈎 + 股*鈎), # eq8
鈎*(-r2^2*鈎 + 2*r2*x2*股 - 2*r2*股^2 + x2^2*鈎 - 2*x2*股*鈎 + 股^2*鈎), # eq9
]
end;
r4 = 2.5/2
r5 = 2/2
iniv = BigFloat[14.5, 7.4, 35.1, 2.5, 26.5, 23.0, 29.7, 44, 56.5]
res = nls(H, ini=iniv)
([14.462681585949788, 7.363220593359079, 35.101665561800985, 2.5077640500378546, 26.50743030867362, 22.966410646176772, 29.67461457867614, 44.09240323596311, 56.50743030867362], true)
金円，水円の直径がそれぞれ 2.5 寸，2 寸のとき，土円の直径は 2*2.5077640500378546 = 5.015528100075709 である。
「術」はこう言っている。
「水円の径を金円の径で割り，1 を加えて，二乗したもので，水円の径の9倍を割る」
おみごと！！
水円径 = 2
金円径 = 2.5
9水円径 / ((sqrt(水円径/金円径) + 1)^2)
5.01552810007571
その他のパラメータは以下のとおりである。
r1 = 14.4627; r2 = 7.36322; x2 = 35.1017; r3 = 2.50776; x3 = 26.5074; x4 = 22.9664; x5 = 29.6746; 鈎 = 44.0924; 股 = 56.5074
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r4 = 2.5/2
r5 = 2/2
(r1, r2, x2, r3, x3, x4, x5, 鈎, 股) = res[1]
@printf("金円，水円の直径がそれぞれ %g，%g のとき，土円の直径は %g である。\n", 2r4, 2r5, 2r3)
@printf("r1 = %g; r2 = %g; x2 = %g; r3 = %g; x3 = %g; x4 = %g; x5 = %g; 鈎 = %g; 股 = %g\n", r1, r2, x2, r3, x3, x4, x5, 鈎, 股)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=1, lw=0.5)
circle(r1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :magenta)
circle(x4, r4, r4, :orange)
circle(x5, r5, r5, :brown)
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/2)
point(x2, r2, "火:r2,(x2,r2)", :blue, :center, delta=-delta/2)
point(x3, r3, " 土:r3,(x3,r3)", :magenta, :left, :bottom, delta=delta/2)
point(x4, r4, "金:r4,(x4,r4) ", :orange, :right, :bottom, delta=delta/2)
point(x5, r5, " 水:r5,(x5,r5) ", :brown, :left, :bottom, delta=delta/2)
point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
end
end;
draw(true)
