算額（その475）

算額（その475）

宮城県石巻市雄勝町 葉山神社 明治17年(1884)

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

外円内に同じ大きさの菱形 2 個と円が入っている。円の直径が 1 寸のとき外円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
菱形の長い方の対角線と短い方の対角線の長さを 2b, 2a
円の半径と中心座標を r1, (x1, y1)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms r0::positive, a::positive, b::positive,
     r1::positive, x1::positive, y1::positive;

eq1 = 2a + 2b - 2r0
eq2 = (r0 - a)^2 + b^2 - r0^2
eq3 = x1^2 + y1^2 - (r0 - r1)^2
eq4 = distance(0, r0 - 2a, b, r0 - a, x1, y1) - r1^2
eq5 = distance(0, r0 - 2a, a, b - r0, x1, y1) - r1^2;

a, b は r0 が分かれば決まる。しかし，r0, x1, y1 は SymPy では有限の時間内に計算できないようなので，a, b も含めて数値解を求める。

solve([eq1, eq2], (a, b))

   1-element Vector{Tuple{Sym, Sym}}:
    (-sqrt(2)*r0/2 + r0, sqrt(2)*r0/2)

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)
   (r0, a, b, x1, y1) = u
   return [
       2*a + 2*b - 2*r0,  # eq1
       b^2 - r0^2 + (-a + r0)^2,  # eq2
       x1^2 + y1^2 - (r0 - r1)^2,  # eq3
       -r1^2 + (y1 - (a^2*y1 - 2*a*b^2 + a*b*x1 + b^2*r0)/(a^2 + b^2))^2 + (-b*(2*a^2 - a*r0 + a*y1 + b*x1)/(a^2 + b^2) + x1)^2,  # eq4
       -r1^2 + (y1 - (-a*(2*a^2 - a*r0 - 2*a*x1 + a*y1 - b*x1 + 2*r0*x1) + y1*(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2))/(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2))^2 + (-a*(4*a^2 + 2*a*b - 6*a*r0 + a*x1 + 2*a*y1 - b*r0 + b*y1 + 2*r0^2 - 2*r0*y1)/(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2) + x1)^2,  # eq5
   ]
end;

r1 = 1//2
iniv = [big"85", 25, 60, 48, 34] ./ 37
iniv = [big"2.297", 0.676, 1.622, 1.297, 0.919]
res = nls(H, ini=iniv);

   外円の直径 = 2.38017;  r0 = 1.19008;    a = 0.348568;  b = 0.841517;  x1= 0.653281;  y1= 0.222351

外円の直径は 2寸3分8毛有奇である。

using Plots

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