https://sangaku0418.hatenablog.com/entry/2024/03/14/214509

福島県白河市南湖南湖神社昭和58年(1983)

http://www.wasan.jp/fukusima/nanko.html
キーワード：円7個，台形
#Julia #SymPy #算額 #和算 #数学

等脚台形の中に大円が内接しており，台形との隙間に 6 個の小円が入っている。小円の直径が 1 寸のとき，大円の直径はいかほどか。

「算額（その731）」から内部の楕円と内円を抜き取った簡易版である。

大円の中心を原点に置く。
大円の半径と中心座標を \(R, (0, 0)\)
台形の右上と右下の頂点の座標を \( (x_b, R), (x_a, -R)\)
小円の半径と中心座標を \(r_2, (x_1, R - r_2), (x_2, y_2), (x_3, r_2)\)
とおき以下の連立方程式を解く。

include("julia-source.txt"); # 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;

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; x_a = 3.96863; x_b = 3.08671; x_1 = 2.64575; x_2 = 3.30719\)
\(y_2 = -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;