算額（その919）

算額（その919）

一一〇 久喜市菖蒲町小林 小林神社 大正5年(1916)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

全円内に楕円と小円をそれぞれ 4 個入れる。楕円の長径，短径がそれぞれ 3 寸，2 寸のとき，全円の直径はいかほどか。

図形を 45°時計方向に回転させ，楕円の中心が x, y 軸上になるようにして以下の連立方程式を解く。
全円の半径と中心座標を R, (0, 0)
楕円の短半径，長半径と中心座標を a, b, (R - a, 0)
小円の半径と中心座標を r, (x, y); y = x
隣り合う 2 個の楕円の接点座標を (x0, y0); y0 = x0
として，以下の連立方程式を解くと，全円の半径 R を求めることができる。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     a::positive, b::positive, x0, y0, x1, y1
y0 =  x0
eq1 = (x0 - R + a)^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*(x0 - R + a)/(a^2*y0) - 1
(R, x0) = solve([eq1, eq2], (R, x0))[2]

   (a + sqrt(a^2 + b^2), b^2/sqrt(a^2 + b^2))

全円の直径は 2(a + sqrt(a^2 + b^2)) である。
楕円の長径，短径がそれぞれ 3 寸，2 寸のとき，a = 1.5, b = 1 なので，6.60555127546399 寸である。

a = 1.5
b = 1
2(a + sqrt(a^2 + b^2))

   6.60555127546399

算額の「問」の答えとしてはここまでである。

---

以下では，図を描くために，小円の半径と中心座標の数値解を求める。

using SymPy
@syms R::positive, r::positive, x::positive,
     a::positive, b::positive, x0, y0, x1, y1
@syms R, r, x, a, b, x0, y0, x1, y1
R = a + sqrt(a^2 + b^2)
y = x
eq3 = (x1 - x)^2 + (y1 - y)^2 - r^2
eq4 = (x1 - R + a)^2/a^2 + y1^2/b^2 - 1
eq5 = -b^2*(x1 - R + a)/(a^2*y1) - (x - x1)/(y - y1)
eq6 = 2x^2 - (R - r)^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, x, x1, y1) = u
   return [
       -r^2 + (-x + x1)^2 + (-x + y1)^2,  # eq3
       -1 + y1^2/b^2 + (x1 - sqrt(a^2 + b^2))^2/a^2,  # eq4
       -(x - x1)/(x - y1) - b^2*(x1 - sqrt(a^2 + b^2))/(a^2*y1),  # eq5
       2*x^2 - (a - r + sqrt(a^2 + b^2))^2,  # eq6
   ]
end;

(a, b) = (3, 2) .// 2
iniv = BigFloat[0.78, 1.78, 2, 1]
res = nls(H, ini=iniv)

   ([0.7824436199506534, 1.7821438606147606, 1.7711417205363233, 0.999777596442297], true)

function draw(more=false)
   pyplot(size=(600, 600), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3/2
   b = 2/2
   R = a + sqrt(a^2 + b^2)
   r = 0.7824436199506534
   @printf("楕円の長径が %g 寸，短径が %g 寸のとき，外円の直径は %g，等円の直径は %g である。\n", 2a, 2b, 2R, 2r)
   @printf("a = %g;  b = %g;  R = %g;  r = %g\n", a, b, R, r)
   plot()
   circle(0, 0, R)
   circle42(0, R - r, r, :blue)
   for theta = 45:90:315
       ox = (R - a)cosd(theta)
       oy = (R - a)sind(theta)
       ellipse(ox, oy, a, b, φ=theta, color=:green)
   end
   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", :green, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
       point(R - r, 0, "等円:r,(R-r,0)", :green, :center, delta=-delta/2)
   end
end;