算額（その1114）

算額（その1114）

百二 岩手県大船渡市猪川町 気仙猪川村安養寺 安政年間か
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，楕円，正方形

正方形の中に楕円 1 個と等円を5個容れる。等円の直径が与えられたとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さを 2a
楕円の長半径と短半径，中心座標を a, b, (0, b)
等円の半径と中心座標を r, (0, 2a - r), (r, b + r), (r, b - r)
右上の等円と楕円の接点座標を (x0, y0)
とおき，連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positice, b::positive, x0::positive,
     y0::positive, r::positive
@syms a, b, x0, y0, r
b = a - r
eq1 = x0^2/a^2 + (y0 - b)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - b - r)) +(y0 - b - r)/x0
eq3 = (x0 - r)^2 + (y0 - b - r)^2 - r^2;

eq1, eq2 から x0，y0 を求める。

(ans_x0, ans_y0) = solve([eq1, eq2], (x0, y0))[2]

   (sqrt(2)*a*sqrt(a^2 - r*(2*a - r + sqrt(2*a^2 - 4*a*r + r^2)))/(2*(a - r)), a - r/2 + sqrt(2*a^2 - 4*a*r + r^2)/2)

eq3 に ans_x0, ans_y0 を代入し,a を求める。

eq13 = eq3(x0 => ans_x0, y0 => ans_y0) |> simplify;
res = solve(eq13, a)[5]

   r*CRootOf(4*x^5 - 16*x^4 + 8*x^3 + 12*x^2 - 11*x + 2, 2)

ans_a = res(r => 1/2).evalf()
#= a =# ans_a |> println

   1.57180973798553

正方形の一辺の長さは，等円の半径の 2*1.57180973798553 = 3.14361947597106 倍である。
等円の直径が 1 寸のとき，正方形の一辺の長さは 3.14361947597106 寸である。

「術」では，正方形の一辺の長さは等円の直径の √7 + 0.5 = 3.1457513110645907 倍であるとしているが，違いがある。

図形描画には不要であるが，接点座標は以下のようになる。

#= x0 =# ans_x0(a => ans_a, r => 1/2).evalf() |> println
#= y0 =# ans_y0(a => ans_a, r => 1/2).evalf() |> println

   0.682602589915843
   2.03727326411408

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