算額（その651）

算額（その651）

長野県飯山市静間大久保 静間観音堂（静観庵）弘化5年(1848)

中村信弥「改訂増補 長野県の算額」県内の算額（177）
http://www.wasan.jp/zoho/zoho.html

日堂薫，宮崎雄也，鷲森勇希，伊藤栄一：「飯山市静間観音堂の算額」の復元
https://sbc21.co.jp/gakkokagaku/2015/27.pdf

直径と2本の弦を隔てて 4 個の等円が外円に内接し，直径と弦に接している。等円の直径を外円の直径で表わせ。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r, x, y1, x0, y2
y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + y2^2 - (R - r)^2
eq2 = r^2 + y1^2 - (R - r)^2
eq3 = dist(0, R, x0, -y0, r, y1) - r^2
eq4 = dist(0, R, x0, -y0, x, y2) - r^2
eq5 = y2/x * (y0 + R)/x0 - 1;

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, y1, x0, y2) = u
   t = -R - sqrt(R^2 - x0^2)
   t2 = x0^2 + t^2
   return [
       x^2 + y2^2 - (R - r)^2,  # eq1
       r^2 + y1^2 - (R - r)^2,  # eq2
       -r^2 + (r - x0*(r*x0 + (-R + y1)*t)/t2)^2 + (-R + y1 - t*(r*x0 + (-R + y1)*t)/t2)^2,  # eq3
       -r^2 + (x - x0*(x*x0 + (-R + y2)*t)/t2)^2 + (-R + y2 - t*(x*x0 + (-R + y2)*t)/t2)^2,  # eq4
       -1 + y2*(R + sqrt(R^2 - x0^2))/(x*x0),  # eq5
   ]
end;
R = 1
iniv = BigFloat[8, 17, -18, 18, 8] ./26
iniv = BigFloat[0.313, 0.637, -0.612, 0.694, 0.257] .* R
res = nls(H, ini=iniv)

   (BigFloat[0.3129063194259740325268103642250497724124515879919851568449522877018154669364085, 0.6371784719628400853411527943465112615089968502655608047796456430352258943786332, -0.6117085589952554644133700611927907136554589335804161767166315804026984221703721, 0.6940076375173001867788131238285698449294652216735132602273056924441179784109038, 0.2571017711954972908126786596219933186205074882118679191481568858007998140638422], true)

等円の直径は外円の直径の 0.312906 倍である。

   r = 0.3129063194259740325268103642250497724124515879919851568449522877018154669364085
   R = 1;  r = 0.312906;  x = 0.637178;  y1 = -0.611709;  y2 = 0.257102;  x0 = 0.694008

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x, y1, x0, y2) = res[1]
   println("r = $r")
   @printf("R = %g;  r = %g;  x = %g;  y1 = %g;  y2 = %g;  x0 = %g\n", R, r, x, y1, y2, x0)
   y0 = -sqrt(R^2 - x0^2) 
   plot([x0, 0, -x0], [y0, R, y0], color=:black, lw=0.5)
   circle(0, 0, R)
   circle(x, y2, r, :blue)
   circle(-x, y2, r, :blue)
   circle(r, y1, r, :blue)
   circle(-r, y1, r, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(x, y2, "(x,y2)")
       point(r, y1, "(r,y1)")
       point(x0, y0, " (x0,y0)", :red)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   end
end;

『「飯山市静間観音堂の算額」の復元』には，記号を上に合わせると，r について以下の 3 次式が示されている。
R = 1 にして方程式を解くと，実数解 1 つが得られ，
-(1133/8 + 15*sqrt(114))^(1/3)/3 + 71/(12*(1133/8 + 15*sqrt(114))^(1/3)) + 5/3
評価すると上と同じく r = 0.312906319425974 が得られる。

@syms R, r
eq = 4r^3 - 20R*r^2 + 57R^2*r - 16R^3;

res2 = solve(eq(R=>1));
res2[3] |> println

   -(1133/8 + 15*sqrt(114))^(1/3)/3 + 71/(12*(1133/8 + 15*sqrt(114))^(1/3)) + 5/3

res2[3].evalf() |> println

   0.312906319425974