算額（その1300）

算額（その1300）

百三 群馬県高崎市八幡町 八幡宮 安政7年(1860)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円3個，長方形

長方形の中に相似な楕円 3 個を容れる。長方形の短辺が与えられたとき，長辺はいかほどか。

長方形の長辺を c
楕円の長半径，短半径，中心座標を a, b, (b, a), (c - a, b)
とおき，以下の連立方程式を解く。
長方形の短辺は 2a = 4b である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     x0::positive, y0::positive;
eq1 = (x0 - b)^2/b^2 + (y0 - a)^2/a^2 - 1
eq2 = (x0 -c + a)^2/a^2 + (y0 - b)^2/b^2 - 1
eq3 = a^2*(x0 - b)/(b^2*(y0 - a)) - b^2*(x0 - c + a)/(a^2*(y0 - 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 Float64.(v), r.f_converged
end;

function H(u)
   (c, x0, y0) = u
   return [
-1 + (-b + x0)^2/b^2 + (-a + y0)^2/a^2,  # eq1
-1 + (-b + y0)^2/b^2 + (a - c + x0)^2/a^2,  # eq2
a^2*(-b + x0)/(b^2*(-a + y0)) - b^2*(a - c + x0)/(a^2*(-b + y0)),  # eq3
   ]
end;
a = 1/2
b = a/2
iniv = BigFloat[2.6475964903602565, 0.8542484563784245, 0.5046106449806166] .*a
res = nls(H, ini=iniv)

   ([1.4708869390890313, 0.4745824757657914, 0.28033924721145365], true)

長方形の長辺は，楕円の長半径の 2.9417738781780627 倍である。
楕円の長半径が 1/2（すなわち，長方形の短辺が 1）のとき，長辺は 1.47089 である。

楕円の長半径，短半径は 0.5, 0.25，接点の座標は (0.474582, 0.280339) である。

function draw(短辺, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 短辺/2
   b = a/2
   (c, x0, y0) = a .*[2.9417738781780627, 0.9491649515315828, 0.5606784944229073]
   @printf("長方形の短辺が %g のとき，長辺は %g である。\n楕円の長半径，短半径は %g, %g，接点の座標は (%g, %g) である。", 短辺, c, a, b, x0, y0)
   plot([0, c, c, 0, 0], [0, 0, 2a, 2a, 0], color=:green, lw=0.5)
   ellipse(b, a, b, a, color=:red)
   ellipse(c - a, b, a, b, color=:red)
   ellipse(c - a, 3b, a, b, color=:red)
   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(b, a, "(b,a)")
       point(c - a, b, "(c-a,b)")
       point(x0, y0, "(x0,y0)")
       point(c, 4b, "(c,2a)", :green, :right, :bottom, delta=delta/2)
   end

end;

draw(1, true)