算額（その2042）

算額（その2042）

(2) 兵庫県三田市藍本大丸 酒垂神社 文化8年(1811)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円2個，菱形，扇

地紙長が 6.4 寸，菱形の横の対角線の長さが 14.664 寸のとき，等円の直径はいかほどか。

扇長（要から先端までの長さ）を R
地紙長（紙の貼られている部分の長さ）を k
菱形の横の対角線の長さを 2a
等円の半径と中心座標を r, (x, y)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive, y::positive,
     a::positive, k::positive
eq1 = x^2 + y^2 - ((R - k) + r)^2
eq2 = dist2(0, R - k, a, sqrt(R^2 -a^2), x, y, r)
eq3 = dist2(0, 0, a, sqrt(R^2 -a^2), x, y, r)
eq4 = sqrt(R^2 -a^2) - (R - k/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 Float64.(v), r.f_converged
end;

function H(u)
   (R, r, x, y) = u
   return [
x^2 + y^2 - (R - k + r)^2,  # eq1
4*R^2*a^2 - 8*R*a^2*k - 8*R*a^2*y + 4*R*a*k*x + 4*a^2*k^2 + 8*a^2*k*y - 4*a^2*r^2 + 4*a^2*y^2 - 4*a*k^2*x - 4*a*k*x*y - k^2*r^2 + k^2*x^2,  # eq2
-4*R^2*r^2 + 4*R^2*x^2 - 8*R*a*x*y + 4*R*k*r^2 - 4*R*k*x^2 - 4*a^2*r^2 + 4*a^2*y^2 + 4*a*k*x*y - k^2*r^2 + k^2*x^2,  # eq3
-R + k/2 + sqrt(R^2 - a^2),  # eq4

   ]
end;

a = 14.664/2
k = 6.4
iniv = BigFloat[10, 1, 2.6, 3.7]
res = nls(H, ini=iniv)

   ([9.999722499999999, 0.9429406670508707, 2.61703497874062, 3.7130737360479014], true)

地紙長が 6.4，菱横が 14.664 のとき，等円の直径は 1.88588 である。

その他のパラメータは以下のとおりである

   a = 7.332; k = 6.4; R = 9.99972; r = 0.942941; x = 2.61703; y = 3.71307

 

function draw(a, k, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r, x, y) = res[1]
   @printf("地紙長が %g，菱横が %g のとき，等円の直径は %g である。\n", k, 2a, 2r)
   @printf("a = %g;  k = %g;  R = %g;  r = %g;  x = %g;  y = %g\n", a, k, R, r, x, y)
   θ = atand(R - k/2, a)
   plot([-a, 0, a, 0, -a], [R - k/2, R, R - k/2, R - k, R - k/2], color=:green, lw=0.5)
   plot!([-a, 0, a], [R - k/2, 0, R - k/2], color=:red, lw=0.5)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ)
   circle(0, 0, R - k, :blue, beginangle=θ, endangle=180 - θ)
   circle2(x, y, r, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, "R", :red, :center, :bottom, delta=delta)
       point(0, R - k/2, "R-k/2", :green, :center, :bottom, delta=delta)
       point(a, R - k/2, "(a,R-k/2)  ", :green, :right, :vcenter)
       point(0, R - k, "R-k", :green, :center, :bottom, delta=delta)
       point(x, y, "r,(x,y)", :magenta, :center, delta=-delta/2)
       dimension_line(2delta, -2delta, a + 2delta, R - k/2 - 2delta, "  扇長", :red, :left)
       dimension_line(a*(R - k)/R + 4delta, (R - k/2)*(R - k)/R - 4delta, a + 4delta, R - k/2 - 4delta, "  地紙長", :green, :left)
   end
end;

draw(14.664/2, 6.4, true)