算額（その417）
北海道函館市住吉町 住吉大明神（住三吉神社） 文化2年(1805)
中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html
五芒星内に 5 個の等円を入れる。内側の正五角形の一辺の長さが 1 尺のとき，等円の直径を求めよ。
半径 1/(2*sind(36)) の円に内接する正五角形の一辺の長さが 1 である。
等円の半径と中心座標を r, (x, y) とし，以下の方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms x::positive, y::positive, r::positive, a::positive;
(cosd36, sind36, cosd18, sind18) = (cosd(Sym(36)), sind(Sym(36)), cosd(Sym(18)), sind(Sym(18)))
eq1 = distance(0, a, a*cosd18, a*sind18, x, y) - r^2
eq2 = distance(0, a, -a*cosd18, a*sind18, x, y) - r^2
eq3 = distance(a*sind36, -a*cosd36, a*cosd18, a*sind18, x, y) - r^2;
# res = solve([eq1, eq2, eq3], (x, y, r))
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 v, r.f_converged
end;
function H(u)
(r, x, y) = u
return [
-r^2 + (y - (-5*sqrt(2)*x - sqrt(10)*x + sqrt(sqrt(5) + 5)*(sqrt(5)*a + 3*a - sqrt(5)*y + 5*y))/(8*sqrt(sqrt(5) + 5)))^2 + (-a*sqrt(2*sqrt(5) + 10)/8 - sqrt(5)*x/8 + 5*x/8 + y*sqrt(2*sqrt(5) + 10)/8)^2, # eq1
-r^2 + (a*sqrt(2*sqrt(5) + 10)/8 - sqrt(5)*x/8 + 5*x/8 - y*sqrt(2*sqrt(5) + 10)/8)^2 + (-3*a/8 - sqrt(5)*a/8 - sqrt(2)*x*sqrt(sqrt(5) + 5)/8 + sqrt(5)*y/8 + 3*y/8)^2, # eq2
-r^2 + (x - (-5*sqrt(10)*a/2 + 5*sqrt(2)*a/2 - 5*x*sqrt(sqrt(5) + 5) - 2*x*sqrt(25 - 5*sqrt(5)) + 5*x*sqrt(5 - sqrt(5)) + 2*x*sqrt(5*sqrt(5) + 25) - 5*sqrt(10)*y + 10*sqrt(2)*y)/(-15*sqrt(sqrt(5) + 5) + sqrt(5)*sqrt(5 - sqrt(5)) + sqrt(5)*sqrt(sqrt(5) + 5) + 15*sqrt(5 - sqrt(5))))^2 + (y - (-15*a*sqrt(5 - sqrt(5))/4 - a*sqrt(5*sqrt(5) + 25)/4 - a*sqrt(25 - 5*sqrt(5))/4 + 15*a*sqrt(sqrt(5) + 5)/4 - 5*sqrt(10)*x + 10*sqrt(2)*x - 5*y*sqrt(sqrt(5) + 5) + 5*y*sqrt(5 - sqrt(5)))/(-15*sqrt(sqrt(5) + 5) + sqrt(5)*sqrt(5 - sqrt(5)) + sqrt(5)*sqrt(sqrt(5) + 5) + 15*sqrt(5 - sqrt(5))))^2, # eq3
]
end;
a = 1/2sind(36)
iniv = [big"0.36", 0.62, 0.85]
res = nls(H, ini=iniv);
println(res);
(BigFloat[0.3632712640026803262857421710703028576320274738964989772105963682762334090546295, 0.6180339887498946353591452335230846232288953659873432429140776623397096012159547, 0.8506508083520398708822131798566121390753734825678071178183687667072162785925422], true)
r = 0.363271; x = 0.618034; y = 0.850651; a = 0.850651
正5角形の一辺の長さを 1 としたとき，等円の直径 = 0.726543
等円の半径は 0.36327126400268034，直径は 0.7265425280053607 である。
なお，0.7265425280053607 = √(5 - 2√5) = tand(36) である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(cosd36, sind36, cosd18, sind18) =
((1 + √5)/4, √(10 - 2√5)/4, √(10 + 2√5)/4, (√5 - 1)/4)
a = 1/2sind(36)
ax = zeros(5)
ay = zeros(5)
r0 = a*cosd36 + a*sind(36)/sind(18)*cosd(18)
for i = 1:5
ax[i] = cosd(72i-18)r0
ay[i] = sind(72i-18)r0
end
(r, x, y) = res[1]
@printf("r = %g; x = %g; y = %g; a = %g\n", r, x, y, a)
@printf("正5角形の一辺の長さを 1 としたとき，等円の直径 = %g\n", 2r)
plot(a.*[0, -cosd18, -sind36, sind36, cosd18, 0], a.*[1, sind18, -cosd36, -cosd36, sind18, 1], color=:green, lw=1)
plot!(ax[[1, 3, 5, 2, 4, 1]], ay[[1, 3, 5, 2, 4, 1]], color=:blue, lw=0.5)
rotate(x, y, r, angle=72)
circle(0, 0, r0, :gray85)
circle(0, 0, a, :gray85)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r0, 0, " r0", :gray, delta=-delta)
point(x, y, "r,(x,y)", :red, :center, delta=-delta)
point(0, a, " a", :blue, :left, :vcenter)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
