算額（その1095）

算額（その1095）

五十七 岩手県花泉町 花泉天満宮 文政13年(1830)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.http://www.wasan.jp/yamamura/yamamura.html
キーワード：円9個，外円，楕円2個，正方形

正方形の中に内接する大円，2 個の楕円，8 個の等円を容れる。正方形の一辺の長さが 1 寸のとき，楕円の短径はいかほどか。

正方形の一辺の長さを 2a
等円の半径と中心座標を r, (a - r, a - r), (a - r - √2r, a - r - √2r)
楕円の長半径と短半径，中心座標を 0, 0, (a, b)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r::positive, x0::positive, y0::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (a - r -　√Sym(2)r - x0)/(a - r -　√Sym(2)r - y0)
eq2 = eq2 |> simplify |> numerator
eq3 = 2(a - r)^2 - (a + r)^2 |> expand
eq4 = (a - r -　√Sym(2)r - x0)^2 + (a - r -　√Sym(2)r - y0)^2 - r^2 |> expand;

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)
   (b, r, x0, y0) = u
   return [
       -1 + y0^2/b^2 + x0^2/a^2,  # eq1
       a^2*y0*(-a + r + sqrt(2)*r + x0) - b^2*x0*(-a + r + sqrt(2)*r + y0),  # eq2
       a^2 - 6*a*r + r^2,  # eq3
       2*a^2 - 4*sqrt(2)*a*r - 4*a*r - 2*a*x0 - 2*a*y0 + 5*r^2 + 4*sqrt(2)*r^2 + 2*r*x0 + 2*sqrt(2)*r*x0 + 2*r*y0 + 2*sqrt(2)*r*y0 + x0^2 + y0^2,  # eq4
   ]
end;

a = 1/2
iniv = BigFloat[0.25, 0.086, 0.267, 0.211]
res = nls(H, ini=iniv)

   ([0.24968850668279452, 0.08578643762690495, 0.2670793193663204, 0.2110827336905282], true)

正方形の一辺の長さが 1 寸のとき，楕円の短径は 0.49937701336558904 寸である。

「答」では 二分七厘五毛になっているが？

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1/2
   (b, r, x0, y0) = res[1]
   @printf("正方形の一辺の長さが %g のとき，楕円の短径は %g である。\n", 2a, 2b)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle(0, 0, a, :magenta)
   circle4(a - r, a - r, r)
   circle4(a - r - √2r, a - r - √2r, r)
   ellipse(0, 0, a, b, color=:green)
   ellipse(0, 0, b, a, color=:green)
   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(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta, deltax=-2delta)
       point(a - r, a - r, "等円:r\n(a-r,a-r)", :red, :center, :bottom, delta=delta/2)
       point(a - r - √2r, a - r - √2r, "等円:r\n(a-r-√2r,a-r2-√r)", :red, :center, :bottom, delta=delta/2)
       point(a - r - √2r, a - r - √2r)
       point(0, b, "b", :green, :center, :bottom, delta=delta/2)
       point(0, a, "a", :green, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end
end;