https://sangaku0418.hatenablog.com/entry/2024/10/31/002422

七十三群馬県安中市下後閑威徳神社嘉永3年(1850)

群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，楕円，正方形
#Julia #SymPy #算額 #和算 #数学

第六問
「群馬の算額」によれば，「十問あるが文章は不明」とのことであるが，以下のようなものでもあろう。

正方形の中に楕円と 5 個の等円を容れる。正方形の一辺の長さが与えられたとき，等円の直径を求めよ。

または，等円の直径が与えられたとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さを \(2a\)
等円の半径と中心座標を \(r, (0, 0), (a - r, a - r)\)
楕円の長半径と短半径，中心座標を \(a, b, (0, 0)\)
楕円と等円の接点の座標を \( (x_0, y_0)\)
とおき，以下の連立方程式を解く。

しかし，SymPy の能力的に解析解を求めることができないので，\(a = 1/2\) とおき，数値解を求める（得られる図形は相似なので，一般性を損なわない）。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms a, b, r, x0, y0
b = r
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0*(a - r - y0) + (a - r - x0)*(a^2*y0)
eq3 = (x0 - (a - r))^2 + (y0 - (a - r))^2 - r^2;

function H(u)
(r, x0, y0) = u
return [
-1 + y0^2/r^2 + x0^2/a^2, # eq1
a^2*y0*(a - r - x0) - r^2*x0*(a - r - y0), # eq2
-r^2 + (-a + r + x0)^2 + (-a + r + y0)^2, # eq3
]
end;

a = 1/2
iniv = BigFloat[0.18, 0.28, 0.15]
res = nls(H, ini=iniv)

([0.1785016026524881, 0.27969014973090045, 0.14796196722191884], true)

正方形の一辺の長さが 1 のとき，等円の直径は 0.357003205304976 である。
（等円の直径が 1 のとき，正方形の一辺の長さは 2.8010952986983217 である。）

描画関数プログラムのソースを見る

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 1/2
(r, x0, y0) = [0.1785016026524881, 0.27969014973090045, 0.14796196722191884]
b = r
@printf("正方形の一辺の長さが %g のとき，等円の直径は %.15g である。\n", 2a, 2r)
@printf("r = %g; a = %g; b = %g; x0 = %g; y0 = %g\n", r, a, b, x0, y0)
plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
circle4(r - a, a - r, r)
circle(0, 0, r)
ellipse(0, 0, a, r, color=: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(a, a, "(a,a)", :blue, :right, :bottom, delta=delta/2)
point(a - r, a - r, "(a-r,a-r)", :red, :center, delta=-delta/2)
point(0, b, " b=r", :red, :center, :bottom, delta=delta/2)
end
end;