算額（その1107）
二 岩手県花巻市大田 清水寺 嘉永三年(1850)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，二等辺三角形
二等辺三角形の中に，中円 4 個，小円 1 個を容れる。中円 3 個に外接する大円の直径を求めよ。
二等辺三角形の底辺の長さを 2a，高さを h
大円の半径と中心座標を r0, (x0, r1*(1 + √3))
中円の半径と中心座標を r1, (r1, r1), (0, r1*(1 + √3)), (0, r1*(3 + √3))
小円の半径と中心座標を r2, (0, r1*(4 + √3) + r2)
とおき，以下の連立方程式を解く。
include("julia-source.txt")
using SymPy
@syms a::positive, b::positive,
r1::positive, r2::positive, r0::positive, x0::positive
eq1 = dist2(0, b, a, 0, 0, (4 + √Sym(3))r1 + r2, r2)
eq2 = dist2(0, b, a, 0, 0, (3 + √Sym(3))r1, r1)
eq3 = dist2(0, b, a, 0, r1, r1, r1)
eq4 = x0^2 + r1^2 - (r0 + r1)^2
eq5 = (x0 - r1)^2 + r1^2*(1 + √Sym(3))^2 - (r0 + r1)^2;
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)
(a, b, r1, r0, x0) = u
return [
(a^2*b^2 - 8*a^2*b*r1 - 2*sqrt(3)*a^2*b*r1 - 2*a^2*b*r2 + 8*sqrt(3)*a^2*r1^2 + 19*a^2*r1^2 + 2*sqrt(3)*a^2*r1*r2 + 8*a^2*r1*r2 - b^2*r2^2)/(a^2 + b^2), # eq1
(a^2*b^2 - 6*a^2*b*r1 - 2*sqrt(3)*a^2*b*r1 + 6*sqrt(3)*a^2*r1^2 + 11*a^2*r1^2 - b^2*r1^2)/(a^2 + b^2), # eq2
a*b*(a*b - 2*a*r1 - 2*b*r1 + 2*r1^2), # eq3
r1^2 + x0^2 - (r0 + r1)^2, # eq4
r1^2*(1 + sqrt(3))^2 - (r0 + r1)^2 + (-r1 + x0)^2, # eq5
]
end;
r2 = 1/2
iniv = BigFloat[1.96, 7.3, 0.85, 2.4, 3.17]
res = nls(H, ini=iniv)
([1.9558948090462631, 7.299498801620881, 0.8491981862085499, 2.431851652578137, 3.1692507766256446], true)
大円の半径は，小円の半径の 2.431851652578137/0.5 = 4.863703305156274 倍である。
小円の直径が 1 寸のとき，大円の直径は 4.863703305156274 寸である。
「術」では，「置四十八開平方加八個開平方加一個乗小円径」とあるが，これを山村は √48 + 1 倍とケアレスミスをしている。
術は，√(√48 + 8) + 1 倍といっており，これは二重根号を外して 1 + √2 + √6 となるが，いずれにしろそれは 4.863703305156273 倍ということで，上の数値解と一致している。
sqrt(√Sym(48) + 8)+1 |>sympy.sqrtdenest |> println
1 + sqrt(2) + sqrt(6) |> println
1 + sqrt(2) + sqrt(6)
4.863703305156273
その他のパラメータは以下のとおりである。
r2 = 0.5; a = 1.95589; b = 7.2995; r1 = 0.849198; r0 = 2.43185; x0 = 3.16925
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 1/2
(a, b, r1, r0, x0) = res[1]
@printf("小円の直径が %g のとき，大円の直径は %g である。\n", 2r2, 2r0)
@printf("r2 = %g; a = %g; b = %g; r1 = %g; r0 = %g; x0 = %g\n", r2, a, b, r1, r0, x0)
plot([a, 0, -a, a], [0, b, 0, 0], color=:magenta, lw=0.5)
circle2(r1, r1, r1)
circle(0, r1*(1 + √3), r1, :red)
circle(0, r1*(3 + √3), r1, :red)
circle(0, r1*(4 + √3)+ r2, r2, :blue)
circle(x0, r1*(2 + √3), r0, :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\n(r1,r1)", :red, :center, delta=-delta/2)
point(0, r1*(1 + √3), "中円:r1\n(0,r1*(1+√3))", :red, :center, delta=-delta/2)
point(0, r1*(3 + √3), "中円:r1\n(0,r1*(3+√3))", :red, :center, delta=-delta/2)
point(0, r1*(4 + √3) + r2, "小円:r1\n(0,r1*(4+√3)+r2)", :red, :center, delta=-delta/2)
point(x0, r1*(2 + √3), "大円:r0,(x0,r1*(2+√3))", :green, :center, delta=-delta/2)
point(a, 0, "a", :magenta, :left, :bottom, delta=delta/2)
point(0, b, " b", :magenta, :left, :vcenter)
end
end;
