https://sangaku0418.hatenablog.com/entry/2024/07/11/211230

二〇埼玉県武州不動岡村不動堂文政元年(1818)

埼玉県立図書館：埼玉県史料集第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円6個，外円，正五角形，斜線
#Julia #SymPy #算額 #和算 #数学

外円の中に等円 5 個と正五角形を容れる。外円の直径が 15 寸，等円の直径が 3 寸のとき，正五角形の一辺の長さはいかほどか。

外円の半径と中心座標を \(R, (0, 0)\)
等円の半径と中心座標を \(r, (x, y)\)
正五角形の外接円の半径を \(a\)
とおき，以下の連立方程式を解く。

正五角形の一辺の長さは \(2a\ \sin(\pi/5)\) である。

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

using SymPy

@syms R::positive, a::positive, r::positive, x::positive, y::positive,
x1::positive, y1::positive,
x2::positive, y2::positive,
x0::positive, y0::positive
(x1, y1) = (a*cosd(Sym(18)), a*sind(Sym(18)))
(x2, y2) = (a*cosd(Sym(54)), -a*sind(Sym(54)))
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist(x1, y1, 0, a, x, y) - r^2
eq3 = dist(x1, y1, x2, y2, x, y) - r^2;

function H(u)
(a, x, y) = u
return [
   x^2 + y^2 - (R - r)^2, # eq1
   sqrt(5)*a^2/8 + 3*a^2/8 - a*x*sqrt(2*sqrt(5) + 10)/4 - 3*a*y/4 - sqrt(5)*a*y/4 - r^2 - sqrt(5)*x^2/8 + 5*x^2/8 + x*y*sqrt(2*sqrt(5) + 10)/4 + sqrt(5)*y^2/8 + 3*y^2/8, # eq2
   sqrt(5)*a^2/8 + 3*a^2/8 - 9*a*x*sqrt(10*sqrt(5) + 50)/16 - 19*a*x*sqrt(2*sqrt(5) + 10)/16 + 13*a*x*sqrt(10 - 2*sqrt(5))/8 + 3*a*x*sqrt(50 - 10*sqrt(5))/4 + a*y/2 - r^2 + sqrt(5)*x^2/8 + 5*x^2/8 - 11*x*y*sqrt(10*sqrt(5) + 50)/8 - 23*x*y*sqrt(2*sqrt(5) + 10)/8 + 37*x*y*sqrt(10 - 2*sqrt(5))/8 + 17*x*y*sqrt(50 - 10*sqrt(5))/8 - sqrt(5)*y^2/8 + 3*y^2/8, # eq3
]
end;

R = 15/2
r = 3/2
iniv = BigFloat[5, 4, 4]
res = nls(H, ini=iniv)

([5.430242979853682, 4.375871302877152, 4.105088347484888], true)

正五角形の外接円の半径は 5.430242979853682 なので，正五角形の一辺の長さは 5.430242979853682*sind(36)*2 = 6.3836334798454555 寸である。

res[1][1]*sind(36)*2

6.3836334798454555

図を描くために必要な座標値を計算する。

@syms R
eq11 = x0^2 + y0^2 - R^2
eq12 = y0 - y2 - (y1 - y2)/(x1 - x2)*(x0 - x2)
res2 = solve([eq11, eq12], (x0, y0))[1]

(-a*sqrt(50 - 10*sqrt(5))/40 + a*sqrt(10*sqrt(5) + 50)/40 + a*sqrt(10 - 2*sqrt(5))/8 + a*sqrt(2*sqrt(5) + 10)/8 + (a/4 - sqrt(10)*sqrt(-4*sqrt(5)*R^2 + 20*R^2 - 5*a^2 - sqrt(5)*a^2)/(2*(-sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 10)))*(-sqrt(10*sqrt(5) + 50)/10 + sqrt(50 - 10*sqrt(5))/10), -a/4 + sqrt(-40*sqrt(5)*R^2 + 200*R^2 - 50*a^2 - 10*sqrt(5)*a^2)/16 + sqrt(-8*sqrt(5)*R^2 + 40*R^2 - 10*a^2 - 2*sqrt(5)*a^2)/16)

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

function draw(R, r, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, x, y) = [5.430242979853682, 4.375871302877152, 4.105088347484888]
(x1, y1) = (a*cosd(18), a*sind(18))
(x2, y2) = (a*cosd(54), -a*sind(54))
(x0, y0) = (-a*sqrt(50 - 10*sqrt(5))/40 + a*sqrt(10*sqrt(5) + 50)/40 + a*sqrt(10 - 2*sqrt(5))/8 + a*sqrt(2*sqrt(5) + 10)/8 + (a/4 - sqrt(10)*sqrt(-4*sqrt(5)*R^2 + 20*R^2 - 5*a^2 - sqrt(5)*a^2)/(2*(-sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 10)))*(-sqrt(10*sqrt(5) + 50)/10 + sqrt(50 - 10*sqrt(5))/10), -a/4 + sqrt(-40*sqrt(5)*R^2 + 200*R^2 - 50*a^2 - 10*sqrt(5)*a^2)/16 + sqrt(-8*sqrt(5)*R^2 + 40*R^2 - 10*a^2 - 2*sqrt(5)*a^2)/16)
θ1 = -54:72:306
θ2 = atand(y0, x0):72:432
plot()
circle(0, 0, R)
rotate(x, y, r, angle=72, :blue)
for i = 1:5
segment(a*cosd(θ1[i]), a*sind(θ1[i]), R*cosd(θ2[i]), R*sind(θ2[i]), :green)
end
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(x1, y1, " (x1,y1)", :green, :left, :vcenter)
point(x2, y2, "(x2,y2) ", :green, :right, :bottom, delta=delta/2)
point(0, a, " a", :green, :left, :bottom, delta=delta/2)
point(0, R, " R", :red, :left, :bottom, delta=delta/2)
point(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
end
end;