算額（その418）

算額（その418）

福島県田村市船引町北移東鳥堂 東鳥堂（東鳥神社）観音堂 文久2年(1862)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

団扇の中に弦の中点を通る2線分がある。3 個の等円は外円に接し，弦と2線分にも接している。弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，等円の直径を求めよ。

等円の半径と中心座標を r, (x, a + r) とし，以下の方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r::positive, x::positive,
     a::positive, b::positive, 矢::positive, 弦::positive;

矢 = r0 - a
弦 = 2sqrt(r0^2 - a^2)
eq1 = 弦 + 矢 - 16
eq2 = 弦*矢 - 48
eq3 = distance(0, a, b, sqrt(r0^2 - b^2), x, a + r) - r^2
eq4 = distance(0, a, b, sqrt(r0^2 - b^2), 0, r0 - r) - r^2
eq5 = x^2 + (a + r)^2 - (r0 - r)^2;
# solve([eq1, eq2, eq3, eq4, eq5], (r, r0, x, a, b))

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, r0, x, a, b) = u
   return [
       -a + r0 + 2*sqrt(-a^2 + r0^2) - 16,  # eq1
       2*(-a + r0)*sqrt(-a^2 + r0^2) - 48,  # eq2
       -r^2 + (-b*(-a*r + b*x + r*sqrt(-b^2 + r0^2))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) + x)^2 + (a + r - (a*(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) - (a - sqrt(-b^2 + r0^2))*(-a*r + b*x + r*sqrt(-b^2 + r0^2)))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq3
       b^2*(a^2 + a*r - a*r0 - a*sqrt(-b^2 + r0^2) - r*sqrt(-b^2 + r0^2) + r0*sqrt(-b^2 + r0^2))^2/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2)^2 - r^2 + (-r + r0 - (-a^2*r + a^2*r0 + a*b^2 + 2*a*r*sqrt(-b^2 + r0^2) - 2*a*r0*sqrt(-b^2 + r0^2) + b^2*r - b^2*r0 - r*r0^2 + r0^3)/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq4
       x^2 + (a + r)^2 - (-r + r0)^2,  # eq5
   ]
end;

iniv = [big"1.6", 5, 3, 2.5, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.499999999999999999999999999999999999999999999999999999999999999999999999999948, 6.499999999999999999999999999999999999999999999999999999999999999999999999999862, 3.000000000000000000000000000000000000000000000000000000000000000000000000000138, 2.499999999999999999999999999999999999999999999999999999999999999999999999999862, 2.594733192202055198398672253319262240463466167190260192229005823351113326367367], true)

r = 1.5;  r0 = 6.5;  x = 3;  a = 2.5;  b = 2.59473
矢 = 4;  弦 = 12
等円の直径 = 3

弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，弦は 12 寸，矢は 4 寸，等円の直径は 3 寸である。

算額に描かれている図は，「弦は 12 寸，矢は 4 寸」に対応していない。実際に図を描いてみると，とても「団扇」とは言えないものである。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r0, x, a, b) = res[1]
   矢 = r0 - a
   弦 = 2sqrt(r0^2 - a^2)
   @printf("r = %g;  r0 = %g;  x = %g;  a = %g;  b = %g\n", r, r0, x, a, b)
   @printf("矢 = %g;  弦 = %g\n", 矢, 弦)
   @printf("等円の直径 = %g\n", 2r)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r, r, :blue)
   circle(x, a + r, r, :blue)
   circle(-x, a + r, r, :blue)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a)
   segment(0, a, b, sqrt(r0^2 - b^2))
   segment(0, a, -b, sqrt(r0^2 - b^2))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x, a + r, "r,(x,a+r)", :blue, :center, delta=-delta)
       point(0, a, " a", :black)
       point(b, √(r0^2-b^2), "(b,√(r0^2-b^2))", :red, :left, :bottom, delta=delta)
       point(√(r0^2-a^2), a, "(√(r0^2-a^2),a)", :red, :right, :top, delta=-delta)
       point(0, r0 - r, " r0-r", :blue, :left, :vcenter)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;