算額（その417）

算額（その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;