算額（その779）
福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html
等脚台形の中に大円が内接しており，台形との隙間に 6 個の小円が入っている。小円の直径が 1 寸のとき，大円の直径はいかほどか。
算額（その731）から内部の楕円と内円を抜き取った簡易版である。
https://blog.goo.ne.jp/r-de-r/e/14641c44eab0e2f8256131152a6cd68b
大円の中心を原点に置く。
大円の半径と中心座標を R, (0, 0)
台形の右上と右下の頂点の座標を (xb, R), (xa, -R)
小円の半径と中心座標を r2, (x1, R - r2), (x2, y2), (x3, r2)
とおき以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, xa::positive, xb::positive, r2::positive, x1::positive, x2::positive, y2::negative, x3::positive
@syms d
eq1 = dist(xa, -R, xb, R, x1, R - r2) - r2^2
eq2 = dist(xa, -R, xb, R, x2, y2) - r2^2
eq3 = dist(0, 0, xa, - R, x2, y2) - r2^2
eq1 = numerator(apart(eq1, d))
eq2 = numerator(apart(eq2, d))
eq3 = numerator(apart(eq3, d))
eq4 = x3^2 + (r2 - R)^2 - (R + r2)^2;
eq5 = x2^2 + y2^2 - (R + r2)^2
eq6 = xa*xb - R^2;
eq7 = x1^2 + (R - r2)^2 - (R + r2)^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)
(R, xa, xb, x1, x2, y2, x3) = u
return [
-r2^2 + (2*R - 2*R*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) - r2)^2 + (x1 - xa - (-xa + xb)*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2, # eq1
-r2^2 + (R - 2*R*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) + y2)^2 + (x2 - xa - (-xa + xb)*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2, # eq2
-r2^2 + (x2 - xa*(-R*y2 + x2*xa)/(R^2 + xa^2))^2 + (R*(-R*y2 + x2*xa)/(R^2 + xa^2) + y2)^2, # eq3
x3^2 + (-R + r2)^2 - (R + r2)^2, # eq4
x2^2 + y2^2 - (R + r2)^2, # eq5
-R^2 + xa*xb, # eq6
x1^2 + (R - r2)^2 - (R + r2)^2, # eq7
]
end;
r2 = 1/2
iniv = BigFloat[7.1, 8, 6.1, 5.3, 6.6, -4.5, 3]
res = nls(H, ini=iniv)
([3.5, 3.968626966596886, 3.086709862908689, 2.6457513110645907, 3.307189138830738, -2.25, 2.6457513110645907], true)
大円の直径は小円の直径の 7 倍である。
小円の直径が 1 寸のとき，大円の直径は 7 寸である。
その他のパラメータは以下のとおりである。
R = 3.5; xa = 3.96863; xb = 3.08671; x1 = 2.64575; x2 = 3.30719; y2 = -2.25; x3 = 2.64575
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 1/2
(R, xa, xb, x1, x2, y2, x3) = res[1]
@printf("R = %g; xa = %g; xb = %g; x1 = %g; x2 = %g; y2 = %g; x3 = %g\n", R, xa, xb, x1, x2, y2, x3)
plot([xa, xb, -xb, -xa, xa], [-R, R, R, -R, -R], color=:blue, lw=0.5)
circle(0, 0, R)
circle2(x3, r2 - R, r2, :green)
circle2(x1, R - r2, r2, :green)
circle2(x2, y2, r2, :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(0, R, "R", :blue, :center, :bottom, delta=delta/2)
point(xa, -R, "(xa,-R)", :black, :right, :bottom, delta=delta/2)
point(xb, R, "(xb,R)", :black, :right, :bottom, delta=delta/2)
point(x1, R - r2, "(x1,R-r2)", :red, :right, :bottom, delta=delta)
point(x2, y2, "小円:r2,(x2,y2)", :red, :right, :vcenter)
point(x3, r2 - R, " (x3,r2-R)", :red, :left, :vcenter)
end
end;
